|
|
EDA365欢迎您登录！
您需要 登录 才可以下载或查看，没有帐号？注册
x
上一篇讲到：方程AX=b的解的讨论（特解、通解、零空间向量等概念）及其MATLAB实现，程序中用到的是mldivide或者A\b的方法（二者相同）来解方程。# T# r8 n+ B) E% }1 \5 E" m
' t. F i. |! h# y' ?$ ^
但实际上运行过程中我们会遇到：当AX=b线性方程组是一个病态方程组；或者A是奇异矩阵（即det(A)=0，不可逆），没法求逆，用不了inv(A)方法只能用A\b，此时MATLAB会报错“矩阵接近奇异值，或者缩放错误。结果可能不准确”…网络上很多人问这个问题怎么解决，其实不是MATLAB的问题，而是MATLAB内置算法的鲁棒性问题，直接用A\b方法无法处理这个棘手的问题。如果没有矩阵论或数值计算方法基础的同学可能会一头雾水。本文借用Moore-Penrose广义逆来解决这个问题，帮助大家理解带奇异矩阵的病态方程组如何解决。
2 U5 ^$ A6 q/ \& P$ I3 c" t% X' @1 I7 s( M
首先我们先来看下mldivide, \在MATLAB中的含义：
: [4 e& q9 \& i& I2 w* C# j. j* x6 U1 ~, j2 {: `4 l$ P
% [/ _- C: W8 \" `
/ y" ~4 C2 ?1 _9 o8 S( _' s
/ d3 U- \ f- K& z4 H& T4 Y: W
1 t, n, ^2 p/ ?: j
& M7 w( ]4 H: ~2 x+ {! O
3 Y& x& y3 x2 i/ c7 ~也就是说，A\b的方法是可以求包含奇异矩阵的方程组的，但是可能会出错。而且错误可能非常离谱。（这个例子告诫我们，不要以为MATLAB算出来的结果都是准确的，MATLAB也不过是调用一些算法进行运算，每个算法都可能存在一些缺陷，无法处理某些极端的情况）。' a4 D. F/ [) j- C: @# E( R; @
- _$ F! B0 L; t( z" G
- L0 r |1 b7 Q, k" o) T. Q这里就涉及到数值计算方法领域矩阵的性态的问题了。我们可以直观来感受一下：
/ K" n. N; M6 z( G$ ~0 l- \% I4 r
% E, y% b4 W! u8 M假设如下方程组：
7 f5 t( y& k0 D- ]7 r$ I* C7 P/ D0 W
% }* h, _* q4 \- s3 e
; P6 A* H" [2 V; S t/ G( j其精确解是(1,1)。/ w) L, X1 j+ P
6 i5 z; e+ t" U+ T- }- N
若对左右边都做一些非常非常微小的变化：" A" y" P% t& Y. M
* k( u1 W! z3 X. @ C @
5 F6 ^+ ~5 G' F
0 ~& p5 }' C3 }! _3 w: T. E其精确解变为：(10,-2)。
6 Z' m( G% _. R( G3 Z) r- [3 a9 ~7 Q8 a& X; p( |
一个非常非常微小的扰动就让方程的解产生巨大的变动，我们称上述方程组是病态方程组，系数矩阵A是病态矩阵。' o5 u: \. ^1 B) f
) U0 H; L1 ]. J0 n# M4 D
如果我们遇到不是方阵的矩阵，或者不能求逆的方阵，要想求解AX=b，避免奇异值导致MATLAB产生错误的情况，我们可以采用“伪逆”来帮助我们解决这个问题。
( c' C; T; a" m$ f$ C" A
, Q/ q7 V7 D9 Q H7 Z& B广义逆矩阵：6 T" c% }, n+ A1 V
# X1 x& F. P) @$ {对任意一个矩阵A，提出四个条件：
a3 O0 B! i( w
9 H" J& p% H. S' F' I3 p
4 R3 v& ]: d# x4 ]3 w( L
0 K1 b3 a2 R9 n6 a1 m5 d6 N* N
如果存在矩阵G满足上述的一部分或全部条件，G就可以称为A的广义逆矩阵。最常用的四种广义逆矩阵定义如下：
1 r7 l) `4 ~- x7 K) |& X/ n2 F; |: |" h8 ~' l. o1 W
1 H* ^6 v8 p) u0 ]/ k' b% t3 f* g- d. q3 U& j m9 g* v) e
MATLAB中自带的pinv方法，就可以求矩阵的M-P广义逆，即A+矩阵。. v0 R4 l- \3 T8 b. M' ~. q. k$ ~
; z, f) A) q" [& M9 O大家可以查阅官方文档看具体应用实例。* x! u( ~8 ^9 i! a$ g: J) h& C, n7 l. ]
' [9 |/ }2 H# W5 k' |
0 L& M: j) G9 p9 ?. E8 A
) d" @ Z, |0 a
如果我们求出A+，就可以有另一种思路来解AX=b了：8 T; a* b1 B' @$ o/ S
+ v# r4 R. V0 b2 G$ N' K: ^$ P
6 Y: R5 v& Y/ Y7 E& s0 f6 ~$ ]7 P! g! F1 o
这里，通解的表达式还是类似于X=X*+X0的形式，A+b相当于是特解，后面那一项就是带系数的自由解，y可以取任意数，注意维度匹配即可。( y2 v% {1 r: v/ B/ V2 ]
u) n2 L" R- ]7 } i& L% y6 Y所以，大家调用pinv求出M-P广义逆，然后用x=A+b + (I-A+A)y这个式子构造出通解就可以啦！+ ~4 I" J/ {7 q1 q- r
% z5 x/ H, `0 w2 s0 G6 l
' h9 r* R; s& o8 ^ @- s8 K6 i2 H5 ]7 {
* n! b6 J8 R4 {- b8 h. V2 W
# O1 M) {. \) G- m9 P4 v |
|
/1 
发表于 2020-6-5 15:17
收藏
淘帖
支持！
反对！