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