EDA365电子论坛网

标题: Moore-Penrose广义逆：可解决MATLAB报错“矩阵接近奇异值，或者缩放错误。结果可能... [打印本页]
作者: haidaowang    时间: 2020-6-5 15:17
标题: Moore-Penrose广义逆：可解决MATLAB报错“矩阵接近奇异值，或者缩放错误。结果可能...
上一篇讲到：方程AX=b的解的讨论（特解、通解、零空间向量等概念）及其MATLAB实现，程序中用到的是mldivide或者A\b的方法（二者相同）来解方程。
5 J) I, T3 j  Q" ]6 @
# s9 y+ R* Y% t% J3 o- S0 b: ?6 ?但实际上运行过程中我们会遇到：当AX=b线性方程组是一个病态方程组；或者A是奇异矩阵（即det(A)=0，不可逆），没法求逆，用不了inv(A)方法只能用A\b，此时MATLAB会报错“矩阵接近奇异值，或者缩放错误。结果可能不准确”…网络上很多人问这个问题怎么解决，其实不是MATLAB的问题，而是MATLAB内置算法的鲁棒性问题，直接用A\b方法无法处理这个棘手的问题。如果没有矩阵论或数值计算方法基础的同学可能会一头雾水。本文借用Moore-Penrose广义逆来解决这个问题，帮助大家理解带奇异矩阵的病态方程组如何解决。* R# g2 {+ S9 L) B* x, w( _* v
: T9 `- T' _; U; N1 x3 B
首先我们先来看下mldivide, \在MATLAB中的含义：/ p1 G% D2 h9 z  X0 {/ M# P" Z) S
# M' o& `0 s6 R. d2 E
# N# i; e7 r0 F- v4 K9 p3 g  k

. L! E8 R7 `9 j$ R! Q4 ?1 ] 0 d. V% `! X( |5 h0 M

+ g: [) X1 Z1 d) l! h
1 R- s7 n2 x! s# X) a0 d1 Q" b/ ~" k% ^. M5 ]
也就是说，A\b的方法是可以求包含奇异矩阵的方程组的，但是可能会出错。而且错误可能非常离谱。（这个例子告诫我们，不要以为MATLAB算出来的结果都是准确的，MATLAB也不过是调用一些算法进行运算，每个算法都可能存在一些缺陷，无法处理某些极端的情况）。0 y1 `1 @, ~8 }/ x
, f( V! w& U% j( c5 y: H! [
7 C, K6 P6 i9 e: V# R
这里就涉及到数值计算方法领域矩阵的性态的问题了。我们可以直观来感受一下：
( L( V, C0 Q9 r/ U  Y. }7 h, E' ?5 R/ m0 g) K) v  C4 d5 v
假设如下方程组：' j. \: m4 f5 E" G( j0 g

# J  n# l9 K; R
# a6 S7 r* B5 A( p
: m. n# T/ Q& i# @) y! `8 x其精确解是(1,1)。
! n) t* I  W- I9 K4 r
6 G3 h9 A! T  `+ t6 z若对左右边都做一些非常非常微小的变化：
1 D5 {/ ^3 }- T& c/ G
2 _( l5 H& r$ O+ G
' M& K# w, A4 c( ~7 C
; u' R% b: l& G' c4 _; _4 u5 c( T, m其精确解变为：(10,-2)。" `5 e" h8 ]% Q
" f' K5 W' H4 q
一个非常非常微小的扰动就让方程的解产生巨大的变动，我们称上述方程组是病态方程组，系数矩阵A是病态矩阵。- |+ S2 P+ p. c* H* |# C3 z

, Y! X* t+ G8 ^如果我们遇到不是方阵的矩阵，或者不能求逆的方阵，要想求解AX=b，避免奇异值导致MATLAB产生错误的情况，我们可以采用“伪逆”来帮助我们解决这个问题。
6 r2 P% |& _+ g3 Z* K! h  u! _4 V/ E1 x7 ^' F+ c4 Q
广义逆矩阵：
2 A' v4 x7 @! k: A( i, F" o+ i) \7 w& N! w. t) e2 e1 Y- `
对任意一个矩阵A，提出四个条件：
4 R( I  j1 W6 h, `( ~! e* h& F4 S; s9 v. Q

( v9 }" t2 ^8 Q7 C" y6 r& c" e
( l7 p/ g/ C. J( @; x如果存在矩阵G满足上述的一部分或全部条件，G就可以称为A的广义逆矩阵。最常用的四种广义逆矩阵定义如下：
( z/ R/ d/ x8 b( }. |" |% v4 H$ @  {; O
) C: o& j) }- n2 L- K: u3 W: R
* W, Z9 ^6 C5 C& q3 M0 r6 V
MATLAB中自带的pinv方法，就可以求矩阵的M-P广义逆，即A+矩阵。
3 }3 }" |8 u+ O5 f' n: ]- n% T1 T
8 e1 g) j! |7 V" [" Y: y  A% Z大家可以查阅官方文档看具体应用实例。0 H) I& }/ V0 x) A
  b5 ], U* g3 p  g- ]8 R

) b9 H; M, q. Q" L; ~7 g7 z
" n! U! I( Z, O% Y& V如果我们求出A+，就可以有另一种思路来解AX=b了：% |: g* O$ V8 d$ u. \
  h& Y* u; Q- I

# @* S( S% X6 r5 g6 L! N% D/ m0 G
这里，通解的表达式还是类似于X=X*+X0的形式，A+b相当于是特解，后面那一项就是带系数的自由解，y可以取任意数，注意维度匹配即可。" |0 A$ ?7 V2 _; A2 @- s/ T

+ q4 i! A  t4 [所以，大家调用pinv求出M-P广义逆，然后用x=A+b + (I-A+A)y这个式子构造出通解就可以啦！# b& w7 M9 T1 D( H

4 k; d" {6 d4 V$ n8 }  Q
' _- Y# U- L& v8 {5 C: O+ k; A7 f$ b6 {4 ]7 t" _

' ~4 `% }5 |$ {! ?  ]: J
! E8 d, o! s" R) W
作者: youOK    时间: 2020-6-5 16:22
太实用了
作者: yin123    时间: 2020-6-8 17:03
可解决MATLAB报错“矩阵接近奇异值，或者缩放错误



欢迎光临 EDA365电子论坛网 (https://bbs.eda365.com/) Powered by Discuz! X3.2