|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
上一篇讲到:方程AX=b的解的讨论(特解、通解、零空间向量等概念)及其MATLAB实现,程序中用到的是mldivide或者A\b的方法(二者相同)来解方程。
8 J# O6 G" n& O1 H
_( D( z' m) {$ b' l6 `) b" b但实际上运行过程中我们会遇到:当AX=b线性方程组是一个病态方程组;或者A是奇异矩阵(即det(A)=0,不可逆),没法求逆,用不了inv(A)方法只能用A\b,此时MATLAB会报错“矩阵接近奇异值,或者缩放错误。结果可能不准确”…网络上很多人问这个问题怎么解决,其实不是MATLAB的问题,而是MATLAB内置算法的鲁棒性问题,直接用A\b方法无法处理这个棘手的问题。如果没有矩阵论或数值计算方法基础的同学可能会一头雾水。本文借用Moore-Penrose广义逆来解决这个问题,帮助大家理解带奇异矩阵的病态方程组如何解决。
" Q% ]- M1 E+ A, ?# x1 V) @ @" v* F
首先我们先来看下mldivide, \在MATLAB中的含义:5 w9 r4 N6 q: c
- w' t) ~ o# D5 ~1 B9 ^; l8 S
& P& U8 u* u1 g/ a
) m- p) r) y; `$ s. r* h
4 n' C# c& B' Y6 E7 Y" ?
) L5 @6 n: H/ q; x, b: K; _+ N, b, I9 ~ R; _* o
# L) O+ d8 x( i+ @5 T1 E4 j
也就是说,A\b的方法是可以求包含奇异矩阵的方程组的,但是可能会出错。而且错误可能非常离谱。(这个例子告诫我们,不要以为MATLAB算出来的结果都是准确的,MATLAB也不过是调用一些算法进行运算,每个算法都可能存在一些缺陷,无法处理某些极端的情况)。
" L9 S% G- p! ^5 S5 F& w8 V) l
" Y) E# Q3 c1 w$ f! G2 i" @
# w6 y6 B; `6 X这里就涉及到数值计算方法领域矩阵的性态的问题了。我们可以直观来感受一下:' U. ]$ l7 E7 g3 r1 f j% H1 K
; e" T4 m, j1 A1 F% N假设如下方程组:
8 G2 S7 N9 S( }# i* [5 V% k! G) q7 U" S" j c- G7 J6 U4 C
0 y/ s) h% n' h9 {. k! z/ S
+ p+ w& s( @9 M+ A4 ^, _" R其精确解是(1,1)。
! ]' x8 \7 Z+ D2 t8 s7 f: G; X" N/ B2 H6 D
若对左右边都做一些非常非常微小的变化:
! H t/ N& l$ X5 l
' H, R c; }$ y0 y! {
+ z4 ]# o7 ^, ?$ L0 |. I1 R% ]
6 p; e8 L. q9 ]6 U: q T# p其精确解变为:(10,-2)。
- a5 f- W' j5 U& M4 d' @. w- w
4 l; ~2 z9 G: Y0 w& {$ _: o一个非常非常微小的扰动就让方程的解产生巨大的变动,我们称上述方程组是病态方程组,系数矩阵A是病态矩阵。3 R! R3 @+ ?$ B1 Z) G6 I
3 d( b# B4 z; O/ i: R
如果我们遇到不是方阵的矩阵,或者不能求逆的方阵,要想求解AX=b,避免奇异值导致MATLAB产生错误的情况,我们可以采用“伪逆”来帮助我们解决这个问题。
- j) r1 D& c! u) D; Q) G; M2 d" j6 r% p/ H: ]+ b2 Y- K9 N8 F) `
广义逆矩阵:2 Q) S1 G8 k# B& l" ^0 U* ~
9 l; c; @7 X. |3 C
对任意一个矩阵A,提出四个条件:
, @- G7 g4 v1 l1 h1 R/ r6 W! h; X! P# K( ^5 f
, W) E7 r2 b( w0 d2 I6 }" I: g ?( _) @+ a: _ V# K% j" Z {
如果存在矩阵G满足上述的一部分或全部条件,G就可以称为A的广义逆矩阵。最常用的四种广义逆矩阵定义如下:" C. H/ X( A4 h' F" f
* f9 ~* R( b. H) W$ z! ^ _( ~
, _9 P4 u' B3 b- d- E- i" n
) B$ X2 C9 g. v# A2 K
MATLAB中自带的pinv方法,就可以求矩阵的M-P广义逆,即A+矩阵。; Y, I8 s3 \1 ]2 o. F2 C/ s# H
/ M* w! c9 ^) Y" Z% |
大家可以查阅官方文档看具体应用实例。
9 G# a, \: \- ?- v* x t8 X# j( b6 D1 o" G2 {( v
9 h9 @) C/ h N) M
7 L [7 i* b" B' P
如果我们求出A+,就可以有另一种思路来解AX=b了:
. n: \/ w0 {( n3 L/ l+ p- l
' e1 K. ]+ v2 y: p& N$ L7 s
7 c+ U& |2 M- P* _; |* H+ i. A7 _* k
这里,通解的表达式还是类似于X=X*+X0的形式,A+b相当于是特解,后面那一项就是带系数的自由解,y可以取任意数,注意维度匹配即可。8 K: k. f- Q' B8 x7 p
& u+ q' p) I7 o& O) D所以,大家调用pinv求出M-P广义逆,然后用x=A+b + (I-A+A)y这个式子构造出通解就可以啦!+ J$ P% |( F" u% n- ?1 S
+ d3 l) d" I+ n& V
* o4 K7 X4 z8 _$ i" s
; k' i8 g* c2 E% l; @. G7 \! c3 Z; D& p% M
1 Y$ P4 I6 D U: y
|
|