matlab penrose广义逆,Moore-Penrose广义逆：可解决MATLAB报错“矩阵接近奇异值，或者缩放错误。结果可能......

EDA365欢迎您登录！

您需要 登录 才可以下载或查看，没有帐号？注册

x

上一篇讲到：方程AX=b的解的讨论(特解、通解、零空间向量等概念)及其MATLAB实现，程序中用到的是mldivide或者A\b的方法(二者相同)来解方程。6 f& `3 z( p9 b

7 f+ m$ Y8 g  w9 G1 E! s- V* X

但实际上运行过程中我们会遇到：当AX=b线性方程组是一个病态方程组；或者A是奇异矩阵(即det(A)=0，不可逆)，没法求逆，用不了inv(A)方法只能用A\b，此时MATLAB会报错“矩阵接近奇异值，或者缩放错误。结果可能不准确”…网络上很多人问这个问题怎么解决，其实不是MATLAB的问题，而是MATLAB内置算法的鲁棒性问题，直接用A\b方法无法处理这个棘手的问题。如果没有矩阵论或数值计算方法基础的同学可能会一头雾水。本文借用Moore-Penrose广义逆来解决这个问题，帮助大家理解带奇异矩阵的病态方程组如何解决。9 n5 m% D/ x2 t) E; F% u5 K

1 q% }- r( g" j

首先我们先来看下mldivide, \在MATLAB中的含义：2 H, y% F9 S6 L3 ?6 A

, T5 Y; }7 o& Q8 b$ z* q  }

8.png (77.1 KB, 下载次数: 0)

2020-6-5 15:15 上传

" J1 l; y  @4 X) I

% B* f+ m# D/ W$ |& a: f

7.png (104.47 KB, 下载次数: 0)

2020-6-5 15:15 上传

2 ]. o8 d4 Y6 O' A  z! c0 a  X/ f

$ f/ d2 S0 X( ]1 a

! q6 F3 v. q; P8 M6 O) K$ u: {& Q, g. J  g! d

也就是说，A\b的方法是可以求包含奇异矩阵的方程组的，但是可能会出错。而且错误可能非常离谱。(这个例子告诫我们，不要以为MATLAB算出来的结果都是准确的，MATLAB也不过是调用一些算法进行运算，每个算法都可能存在一些缺陷，无法处理某些极端的情况)。

8 ?- H, o! s/ }3 P; I7 }; C# U1 B5 F* B

; s8 v) E) L5 [  ?& v9 M

这里就涉及到数值计算方法领域矩阵的性态的问题了。我们可以直观来感受一下：

4 G' ^7 W* Q  L9 y7 S

$ p- A8 v, G0 M1 c假设如下方程组：! V) @  g% h3 l( ~% D2 ?/ O; c

% D: Y" ^- w$ A5 v

6.png (19.8 KB, 下载次数: 0)

2020-6-5 15:15 上传

, z8 J2 d; N) h, p2 N, }7 o$ v

# ^& D3 P8 T2 u- Q/ h3 p! Z/ w其精确解是(1,1)。% ~. a& h9 `# ~% B

8 p& i9 O3 R2 X- F- k

若对左右边都做一些非常非常微小的变化：: V6 p+ K7 t% r: U% L  \

# {+ L7 u# `3 Y+ Y& r

5.png (20.17 KB, 下载次数: 0)

2020-6-5 15:15 上传

4 Q$ ~+ {5 \3 u

6 [( @0 V! Y1 r其精确解变为：(10,-2)。

l  i$ k! M8 h0 ]2 ^4 a1 B" v8 U

一个非常非常微小的扰动就让方程的解产生巨大的变动，我们称上述方程组是病态方程组，系数矩阵A是病态矩阵。

/ N& p( N3 Y  w9 w1 [* i! a( s# r

) d  ?" C; N2 Q0 R9 m' k4 c; R如果我们遇到不是方阵的矩阵，或者不能求逆的方阵，要想求解AX=b，避免奇异值导致MATLAB产生错误的情况，我们可以采用“伪逆”来帮助我们解决这个问题。

9 F: ~, r; c) S; g) R( Y0 S! |: \7 L3 L( k& _

广义逆矩阵：

% A. ]* X/ R- K: {6 w6 `" Z

$ N6 D/ }- i% _. F0 T# `# F对任意一个矩阵A，提出四个条件：

9 x! O+ {9 f0 b3 {, F

$ B, s* ^7 V6 O# |

4.png (111.24 KB, 下载次数: 0)

2020-6-5 15:15 上传

* |2 Z% e. |$ G" d& R$ O# ~5 c) E

; n7 [. J( f+ f( e1 T2 f7 v+ W

如果存在矩阵G满足上述的一部分或全部条件，G就可以称为A的广义逆矩阵。最常用的四种广义逆矩阵定义如下：0 @5 g" @3 R# z

0 v- q, n( N* `9 p0 Y7 @  K

3.png (114.92 KB, 下载次数: 0)

2020-6-5 15:15 上传

`& Z" L+ ]% @5 a. E

2 E4 ?# V9 i( _& F: t' L

MATLAB中自带的pinv方法，就可以求矩阵的M-P广义逆，即A+矩阵。

" x$ o7 ^, j; O0 L& s1 ^$ a$ I

! z; }1 W* c! r0 @  |大家可以查阅官方文档看具体应用实例。4 i$ W4 J) Z8 S) ~; U+ O8 i; @4 `

1 q& G7 F5 X, k7 N

2.png (122.89 KB, 下载次数: 0)

2020-6-5 15:15 上传

: |; R( s6 y, [* w: \+ l& H1 \7 {( ]5 g( v* H2 }

如果我们求出A+，就可以有另一种思路来解AX=b了：a' X4 h# F8 n0 l% O1 |6 H5 m

h& _& w* y! K

1.png (132.21 KB, 下载次数: 0)

2020-6-5 15:15 上传

2 P  ^( P" ]- C

4 I2 _4 P' V4 f7 ^这里，通解的表达式还是类似于X=X*+X0的形式，A+b相当于是特解，后面那一项就是带系数的自由解，y可以取任意数，注意维度匹配即可。

# r% m4 N  ]) r5 P3 L: b1 p! Y+ n# H; ~

所以，大家调用pinv求出M-P广义逆，然后用x=A+b + (I-A+A)y这个式子构造出通解就可以啦！8 ~. P. n0 a% i4 h1 M! p

$ @# g1 p7 V+ p6 N, M( V

* E& K' ]3 C/ c& ~! m

2 a* _' S" x, N3 l8 u4 G

! T) h5 M5 _7 h1 J; K

0 k, t) s7 O* j& G4 Y: h+ f, t