|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
上一篇讲到:方程AX=b的解的讨论(特解、通解、零空间向量等概念)及其MATLAB实现,程序中用到的是mldivide或者A\b的方法(二者相同)来解方程。
: X% b5 P* u+ a( S1 }& f: w& s: ^+ @6 F. V' e# z; W6 ?& z
但实际上运行过程中我们会遇到:当AX=b线性方程组是一个病态方程组;或者A是奇异矩阵(即det(A)=0,不可逆),没法求逆,用不了inv(A)方法只能用A\b,此时MATLAB会报错“矩阵接近奇异值,或者缩放错误。结果可能不准确”…网络上很多人问这个问题怎么解决,其实不是MATLAB的问题,而是MATLAB内置算法的鲁棒性问题,直接用A\b方法无法处理这个棘手的问题。如果没有矩阵论或数值计算方法基础的同学可能会一头雾水。本文借用Moore-Penrose广义逆来解决这个问题,帮助大家理解带奇异矩阵的病态方程组如何解决。
& o0 b& x) R8 `) s
: j9 h7 _+ e p首先我们先来看下mldivide, \在MATLAB中的含义:
/ {5 |& @5 b/ t1 U1 S1 P: k2 A/ P5 l0 l5 ?4 r H! n
4 U, Q1 b8 T, Q; s) {( ^6 N" {
4 ^: x J$ g8 v
) t# e6 H9 \ j0 \" i) }
( T+ o9 e; T D6 _. G6 w0 [
5 a* A% k! a6 n; T0 E' {! ^
$ T' o% O7 B" ], ^: m也就是说,A\b的方法是可以求包含奇异矩阵的方程组的,但是可能会出错。而且错误可能非常离谱。(这个例子告诫我们,不要以为MATLAB算出来的结果都是准确的,MATLAB也不过是调用一些算法进行运算,每个算法都可能存在一些缺陷,无法处理某些极端的情况)。" W+ U4 C3 z6 @0 V: R7 J8 d
2 @0 o* O+ p. y( L- a
$ b! Q# B) S: l4 ~9 Y0 q- @这里就涉及到数值计算方法领域矩阵的性态的问题了。我们可以直观来感受一下:: u! r( O( X( f% p
, n0 B& x0 N; \. t0 ]1 D
假设如下方程组:$ t( M! \+ r7 p
* ^. J0 N3 j7 s- e+ o' _& e1 E. J# k. }
5 l/ W, B2 h7 z0 C+ c4 ~) a% i3 R! z
! B, B* K: m; s: U' U其精确解是(1,1)。7 |, p/ U8 b# K, Y- g) s# y+ r
) |7 g& ~2 W d
若对左右边都做一些非常非常微小的变化:% C4 q8 h) Z0 Y+ h' P: b$ {
' ~# x+ P- d2 y3 |- L. E1 V
' {# G! m+ y9 {+ L% }; l _* Y l5 G- ]4 a2 W5 ~# I
其精确解变为:(10,-2)。
( ?6 d o2 S7 X1 X# l- P3 _* a3 b d6 K! Z# j2 U. ^
一个非常非常微小的扰动就让方程的解产生巨大的变动,我们称上述方程组是病态方程组,系数矩阵A是病态矩阵。7 w# _7 k: o: k, I) X) r
, \$ S; J( N b* N% ^( }如果我们遇到不是方阵的矩阵,或者不能求逆的方阵,要想求解AX=b,避免奇异值导致MATLAB产生错误的情况,我们可以采用“伪逆”来帮助我们解决这个问题。, o: r i* l$ b' l P: j; o
9 n) n r5 q& a" o+ S) ]+ f2 e
广义逆矩阵:/ Z5 U; F4 d: s6 ~( s. X0 W1 b& X/ i
$ `0 r9 @/ V! [
对任意一个矩阵A,提出四个条件:
?' ]" a1 I$ t$ s3 m4 f4 L- q2 L1 C4 E) C3 B* c( \# I3 d
$ w2 W5 L+ B; u# q$ y0 a; [# U
3 x7 Y* B/ }2 ]# C如果存在矩阵G满足上述的一部分或全部条件,G就可以称为A的广义逆矩阵。最常用的四种广义逆矩阵定义如下:
2 g' s9 ^- k% i h" ~; C% T' p# U# R0 c1 \9 v3 U7 n
$ I% h3 P. }% I; G9 d# ~7 o' a: s& c
MATLAB中自带的pinv方法,就可以求矩阵的M-P广义逆,即A+矩阵。3 M/ U( p1 E4 N# k. F
( D# D2 Z0 m8 O' U大家可以查阅官方文档看具体应用实例。
5 v# w4 C- P1 p1 B% J0 n( p2 a* q) o; f# V. U( w
" |, w$ ~2 [/ p8 _) q+ |% f6 b! }( e& @8 p' _0 J; ?9 b
如果我们求出A+,就可以有另一种思路来解AX=b了:
& Z: j4 g1 F9 T* {& X1 @
4 b, k" Y, H$ p# B3 Z6 V
8 D& d y& v" F' O; {% t8 Z" |
8 I9 l! f3 q2 j' T3 |# f
这里,通解的表达式还是类似于X=X*+X0的形式,A+b相当于是特解,后面那一项就是带系数的自由解,y可以取任意数,注意维度匹配即可。# `# A7 E$ | W9 ~- d" j% u
3 E/ {/ Y9 t& s, j
所以,大家调用pinv求出M-P广义逆,然后用x=A+b + (I-A+A)y这个式子构造出通解就可以啦!4 C# y+ d& H f7 K5 ]& {
4 k1 _1 Z X$ ]
. j- x$ r& d; t5 p
2 {+ Z( q ~6 b) x
" x3 `9 _. H3 g7 k* {1 b' s
! w d: r$ d3 h& P7 i |
|