|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
上一篇讲到:方程AX=b的解的讨论(特解、通解、零空间向量等概念)及其MATLAB实现,程序中用到的是mldivide或者A\b的方法(二者相同)来解方程。
! P3 i" U6 X2 }4 t' q
8 V! i2 y# r$ }* u6 d但实际上运行过程中我们会遇到:当AX=b线性方程组是一个病态方程组;或者A是奇异矩阵(即det(A)=0,不可逆),没法求逆,用不了inv(A)方法只能用A\b,此时MATLAB会报错“矩阵接近奇异值,或者缩放错误。结果可能不准确”…网络上很多人问这个问题怎么解决,其实不是MATLAB的问题,而是MATLAB内置算法的鲁棒性问题,直接用A\b方法无法处理这个棘手的问题。如果没有矩阵论或数值计算方法基础的同学可能会一头雾水。本文借用Moore-Penrose广义逆来解决这个问题,帮助大家理解带奇异矩阵的病态方程组如何解决。
% M5 {' A; p* t. F$ ^
/ ~- G4 b& A. F首先我们先来看下mldivide, \在MATLAB中的含义:
8 D! c9 H" Z, {* R, y" w! V& E6 b
7 P, \8 l9 N. a) Z! J
2 L9 x6 w2 P, s- V5 T0 h% j1 a5 z7 ~0 M. P5 n5 `$ L( E" U* y
+ {2 U$ k2 ]) }. _0 E+ w3 ]! X8 h( g! V) h, [! r+ K \, C1 }3 C% |5 `$ E
4 @- q, N k0 u! O9 l' j8 H
- k* K& y% @% \! ]# Q. |/ Y Z也就是说,A\b的方法是可以求包含奇异矩阵的方程组的,但是可能会出错。而且错误可能非常离谱。(这个例子告诫我们,不要以为MATLAB算出来的结果都是准确的,MATLAB也不过是调用一些算法进行运算,每个算法都可能存在一些缺陷,无法处理某些极端的情况)。
/ d& {0 s5 s9 O) k
^2 b! s) N9 |& o- M4 v7 V: q, z* b$ r7 s2 M; m
这里就涉及到数值计算方法领域矩阵的性态的问题了。我们可以直观来感受一下:
8 M9 h/ _. q7 k7 H6 m- j& \2 ?; a; p
假设如下方程组: z5 m# R! k. |# C( ]8 s+ p! ]
. Q4 \* p; D% ?9 H
: _1 S8 @3 F6 H, ^
+ k* Q' x B( }其精确解是(1,1)。
* A N& z7 l% m, U- A" B& u* |/ V! K' P. |" ?3 C
若对左右边都做一些非常非常微小的变化:
3 g& k$ F1 B, j& W6 p5 z6 o8 t6 v+ y$ d/ A
7 q) g& s4 k$ ?2 u& b5 E0 C1 q' e% j) {3 }7 g; m
其精确解变为:(10,-2)。
' W6 Y2 f4 j( E7 C
& P) D+ a5 `) a& ?6 z一个非常非常微小的扰动就让方程的解产生巨大的变动,我们称上述方程组是病态方程组,系数矩阵A是病态矩阵。
& D t V) |7 q/ @& A- Y0 i2 [# Y+ |$ N/ I
如果我们遇到不是方阵的矩阵,或者不能求逆的方阵,要想求解AX=b,避免奇异值导致MATLAB产生错误的情况,我们可以采用“伪逆”来帮助我们解决这个问题。. w: |9 ~$ k5 [+ S( R, |: W
9 w3 N5 S" M2 V' {4 J: O
广义逆矩阵:
" M: B' T' U( ]; R# V* G$ ~) g, D ?1 l' E: S- G `. w
对任意一个矩阵A,提出四个条件:9 o/ y1 Y8 J( j/ O* T& c8 n+ c
' U; D# l4 n- q* N
- p9 M9 G1 |" P7 i. S1 Y. Q
8 s* q4 P* u* W- O- k% A如果存在矩阵G满足上述的一部分或全部条件,G就可以称为A的广义逆矩阵。最常用的四种广义逆矩阵定义如下:
2 ^8 B) y8 _4 g! A) r
1 W' J U% D x. Y( B, c8 @' t
" {* P6 W' S$ n2 _$ \6 H9 b* _7 L7 d% P, ^0 v3 i
MATLAB中自带的pinv方法,就可以求矩阵的M-P广义逆,即A+矩阵。' k3 `* t& \; i( B# b
9 R8 g2 ^3 Z S& D大家可以查阅官方文档看具体应用实例。8 M# n: i) r1 a" s4 ^
6 V, P: Z, `; b' n' f+ Z! y
* N6 x9 }! G3 l% Z8 d, F: M
& [: K" H5 ?& R4 q' u
如果我们求出A+,就可以有另一种思路来解AX=b了:( Q9 X. w1 z0 a7 v# F/ @
: S% |0 k. M7 }) v
( @5 x; R7 v8 v3 H& f
. I8 L% Q* b' D+ W% {) }; G
这里,通解的表达式还是类似于X=X*+X0的形式,A+b相当于是特解,后面那一项就是带系数的自由解,y可以取任意数,注意维度匹配即可。
- r" ~; z2 V% H0 _; |! b4 _
f$ L2 G3 M: \5 q" K+ \所以,大家调用pinv求出M-P广义逆,然后用x=A+b + (I-A+A)y这个式子构造出通解就可以啦!/ }9 e. F0 @/ Y, ~: ?8 h
8 T9 X8 W3 u5 k0 f7 w( a% B& s# ]4 p3 F
+ M4 ?4 b7 P7 T! |( V9 @& L' k( e& B* a" n# \
% M! a! S! k4 }/ V" S( u |
|