|
|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
上一篇讲到:方程AX=b的解的讨论(特解、通解、零空间向量等概念)及其MATLAB实现,程序中用到的是mldivide或者A\b的方法(二者相同)来解方程。
) @4 f$ {, u% t# Z6 h, ~3 a& ]( n! T/ Y& w
但实际上运行过程中我们会遇到:当AX=b线性方程组是一个病态方程组;或者A是奇异矩阵(即det(A)=0,不可逆),没法求逆,用不了inv(A)方法只能用A\b,此时MATLAB会报错“矩阵接近奇异值,或者缩放错误。结果可能不准确”…网络上很多人问这个问题怎么解决,其实不是MATLAB的问题,而是MATLAB内置算法的鲁棒性问题,直接用A\b方法无法处理这个棘手的问题。如果没有矩阵论或数值计算方法基础的同学可能会一头雾水。本文借用Moore-Penrose广义逆来解决这个问题,帮助大家理解带奇异矩阵的病态方程组如何解决。, G6 e5 U! B: |/ u3 S4 u% w
. l5 Y2 p* u: m首先我们先来看下mldivide, \在MATLAB中的含义:
3 b3 r5 ~( |" t- p$ ?1 v. l# I; ~# k4 W/ N: O! \; m1 m
b( h; _7 I* o4 }+ p
/ d0 Z8 t8 {4 f$ V- y, ~
0 _( v+ @- [( I7 w. E) w% n9 Q! V4 R! ]) j3 k
* d8 J2 s7 M; D9 h
2 W- H9 j" [1 V% M
也就是说,A\b的方法是可以求包含奇异矩阵的方程组的,但是可能会出错。而且错误可能非常离谱。(这个例子告诫我们,不要以为MATLAB算出来的结果都是准确的,MATLAB也不过是调用一些算法进行运算,每个算法都可能存在一些缺陷,无法处理某些极端的情况)。$ e8 _$ ? g) I4 p
+ J- u& {) l0 B
( \0 V- L( n6 `6 D; N0 @3 O8 h这里就涉及到数值计算方法领域矩阵的性态的问题了。我们可以直观来感受一下:
4 S( A+ v$ U9 ^8 ?
8 l% w# b( h- o/ R2 z Z+ B2 @5 v' X假设如下方程组:8 `' d. R# p3 O1 P
$ Y& G y( [1 |5 s, @$ D! Y: R
8 }- s5 e. \% ~3 l$ I+ k; ?5 |0 \8 X: k: n% M1 e; y# ?; t) s* k! T3 @4 q
其精确解是(1,1)。
. D( C' y2 j5 g" ]3 W2 f- o; i
7 X2 G: X" {/ u" m, b) t% T2 A1 L若对左右边都做一些非常非常微小的变化:
4 B# h( I# {$ E7 o* B, d- f6 h" S7 P7 P- i+ Y6 }
0 B8 [# u# S: E. [$ u( j
( y5 r( \8 A- n8 o, A( x其精确解变为:(10,-2)。
, Y& o8 T! k* t; u7 f8 ?: i" j' G/ N" ~: m9 S4 r& k# {3 y
一个非常非常微小的扰动就让方程的解产生巨大的变动,我们称上述方程组是病态方程组,系数矩阵A是病态矩阵。
A3 y4 P- l5 R: o& M: i! e9 S+ m& x# K: @- C& H- Z
如果我们遇到不是方阵的矩阵,或者不能求逆的方阵,要想求解AX=b,避免奇异值导致MATLAB产生错误的情况,我们可以采用“伪逆”来帮助我们解决这个问题。) ^& P. J2 e; P9 G; A2 j! Y
4 H3 ^, Y! b5 {4 c! [
广义逆矩阵:5 O% a- z4 y+ G; A3 {
: W u3 o) h8 N8 V& m
对任意一个矩阵A,提出四个条件:
6 u# h8 R, Q: i
$ e: L9 o t8 D: E
& C+ Z# e( j2 {/ T
! d" I' f. u" H' i
如果存在矩阵G满足上述的一部分或全部条件,G就可以称为A的广义逆矩阵。最常用的四种广义逆矩阵定义如下:9 N# ~& a- N7 f. d8 `( C: M( U
0 P$ t1 U5 Q7 w5 T
* a- X- j0 s$ w1 e1 E1 v1 \: J: ]; J8 n( p |
MATLAB中自带的pinv方法,就可以求矩阵的M-P广义逆,即A+矩阵。6 W* L) b0 v: Q w- J5 b. o
( d1 h9 w* Q7 V: f大家可以查阅官方文档看具体应用实例。
1 T4 a- B+ ?5 a9 N
! ~7 q" Z6 k) d+ U2 B3 x
9 y2 P! l$ }9 n% {9 r7 P" R2 |9 Z+ S; G' Q! X& P) C1 V
如果我们求出A+,就可以有另一种思路来解AX=b了:4 b5 E Z6 Q2 N- v5 o2 j3 \
. F5 O8 [+ W/ G* {* Z
, ~# f4 l8 j2 V5 k$ v/ F& Q! ~
这里,通解的表达式还是类似于X=X*+X0的形式,A+b相当于是特解,后面那一项就是带系数的自由解,y可以取任意数,注意维度匹配即可。
+ d8 A9 O6 r; c+ A7 o* \" z' [7 {! V
所以,大家调用pinv求出M-P广义逆,然后用x=A+b + (I-A+A)y这个式子构造出通解就可以啦!: B3 W* u$ U. {. m" {' i
4 ~9 [8 G! |- u+ ?% m; M
9 M/ x+ s) j8 y! g m5 R9 Y, p
. S( ^" `# v1 h5 G# I D! O6 K& T- d& t6 ?- ?; k/ b
8 j+ q: ^/ B9 R
|
|