|
|
EDA365欢迎您登录！
您需要 登录 才可以下载或查看，没有帐号？注册
x
上一篇讲到：方程AX=b的解的讨论（特解、通解、零空间向量等概念）及其MATLAB实现，程序中用到的是mldivide或者A\b的方法（二者相同）来解方程。
. `: x& J/ q0 } a+ M" Q
8 V; F' M, o/ ?' G' Y1 o6 @但实际上运行过程中我们会遇到：当AX=b线性方程组是一个病态方程组；或者A是奇异矩阵（即det(A)=0，不可逆），没法求逆，用不了inv(A)方法只能用A\b，此时MATLAB会报错“矩阵接近奇异值，或者缩放错误。结果可能不准确”…网络上很多人问这个问题怎么解决，其实不是MATLAB的问题，而是MATLAB内置算法的鲁棒性问题，直接用A\b方法无法处理这个棘手的问题。如果没有矩阵论或数值计算方法基础的同学可能会一头雾水。本文借用Moore-Penrose广义逆来解决这个问题，帮助大家理解带奇异矩阵的病态方程组如何解决。
& y: `; F3 x2 D/ K0 u# d
: y1 o! |6 W3 R' v* {5 A7 ?首先我们先来看下mldivide, \在MATLAB中的含义： v# N3 Z! {% \
, m$ I6 _9 i% U& \
" r1 i. X- ~; o; u8 x2 F
0 m u1 {' l; G( j5 ~0 y
( O& Z9 K+ p6 B
. d8 h4 ^ v# L7 d
/ z0 A9 H4 @" v+ w
0 d7 M8 g- N3 L: w* d( _4 \也就是说，A\b的方法是可以求包含奇异矩阵的方程组的，但是可能会出错。而且错误可能非常离谱。（这个例子告诫我们，不要以为MATLAB算出来的结果都是准确的，MATLAB也不过是调用一些算法进行运算，每个算法都可能存在一些缺陷，无法处理某些极端的情况）。& q- X+ [% w: }# T' C ^
& g& p3 ?6 i& M2 p
3 A! m/ ~" m) O/ W, W
这里就涉及到数值计算方法领域矩阵的性态的问题了。我们可以直观来感受一下：
0 U, p# f2 f4 T; L) M: O7 l7 B0 `4 N& z
假设如下方程组：
& x. w+ P1 l6 {. k. |1 O Q* j5 _5 e
7 u5 t6 O2 b; b7 r6 `
7 ]* V# M1 ~! [5 U1 ?9 c其精确解是(1,1)。
6 f' K6 ~5 F2 W% |7 G* ~% K
1 z$ j- T+ {/ E: o: K4 i若对左右边都做一些非常非常微小的变化：
- i ]0 K) z/ d5 j4 {
" k: G! J, B5 y: S& `0 z% H
* z( ]3 P5 X$ `, a' B
: c& i" x6 R. W1 q }; Q, o其精确解变为：(10,-2)。# |+ @8 }0 J8 L, @* V
/ J h1 U6 K, o: s' o9 L: g6 V: D
一个非常非常微小的扰动就让方程的解产生巨大的变动，我们称上述方程组是病态方程组，系数矩阵A是病态矩阵。
. f" v6 w9 G: r8 w; D8 P8 N; [5 r. A+ O4 k& r$ }& X* O
如果我们遇到不是方阵的矩阵，或者不能求逆的方阵，要想求解AX=b，避免奇异值导致MATLAB产生错误的情况，我们可以采用“伪逆”来帮助我们解决这个问题。1 l8 ?% N& z1 I0 w' G: v+ s
: W# q& ^% P! f [) Y" k* l
广义逆矩阵：, O( r" e$ y7 K/ I5 U
# L6 I& K, m( `. r
对任意一个矩阵A，提出四个条件：
* Y, d' h4 \& h" R* N. E2 W- P
}/ B4 y# a; m8 H& E
) b% k8 {2 H( C- q1 ^" j$ s
6 A: j7 V$ l' U7 x9 Y
如果存在矩阵G满足上述的一部分或全部条件，G就可以称为A的广义逆矩阵。最常用的四种广义逆矩阵定义如下：% v$ r- J' c7 E9 R9 x! `4 h' x/ i
! }5 M4 z$ O* ?. v) K% u
( F$ m/ A, _$ q: m
* X/ d) Q9 a' e3 N
MATLAB中自带的pinv方法，就可以求矩阵的M-P广义逆，即A+矩阵。& E+ f3 J. I8 o3 i
; Z# q. s# `3 A/ A+ Y
大家可以查阅官方文档看具体应用实例。2 d: Q2 K9 j( Q; H v* M
- f8 V# ^0 r1 a5 n3 Z# G0 U5 z
2 F! i" C; i0 v) _% l
2 N( H* G, C( `* c4 d5 M如果我们求出A+，就可以有另一种思路来解AX=b了：( w# ]+ D" w1 }# h5 E. n2 s
8 X; y- D3 p- ^* l8 Y. ]9 c: R
/ L% _ A0 e3 B7 ^" h9 M& w$ k/ P
6 \1 U% U5 u; t" @8 f; v- n
这里，通解的表达式还是类似于X=X*+X0的形式，A+b相当于是特解，后面那一项就是带系数的自由解，y可以取任意数，注意维度匹配即可。2 x: m9 I4 b6 H+ E; b: I( ]" k, c
9 R" j1 }& G h) Z; v0 K+ s
所以，大家调用pinv求出M-P广义逆，然后用x=A+b + (I-A+A)y这个式子构造出通解就可以啦！+ i+ i5 l, q* S$ U6 n9 f O& f
, w( Y$ f6 Z8 S" G/ F5 q9 E! T4 H9 p$ U" H
) N) D; `0 S: B& x" f4 p. @3 L) E; D# p
: L( M8 `) B. N/ Q U
|
|
/1 
发表于 2020-6-5 15:17
收藏
淘帖
支持！
反对！