|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
上一篇讲到:方程AX=b的解的讨论(特解、通解、零空间向量等概念)及其MATLAB实现,程序中用到的是mldivide或者A\b的方法(二者相同)来解方程。4 k- D( S! I9 }0 P
1 \, O, u! E$ F5 s1 |& H8 o0 Q
但实际上运行过程中我们会遇到:当AX=b线性方程组是一个病态方程组;或者A是奇异矩阵(即det(A)=0,不可逆),没法求逆,用不了inv(A)方法只能用A\b,此时MATLAB会报错“矩阵接近奇异值,或者缩放错误。结果可能不准确”…网络上很多人问这个问题怎么解决,其实不是MATLAB的问题,而是MATLAB内置算法的鲁棒性问题,直接用A\b方法无法处理这个棘手的问题。如果没有矩阵论或数值计算方法基础的同学可能会一头雾水。本文借用Moore-Penrose广义逆来解决这个问题,帮助大家理解带奇异矩阵的病态方程组如何解决。( s" `2 @$ ^- u
: ~, J" C# W% j% y1 d/ S* s首先我们先来看下mldivide, \在MATLAB中的含义:
+ S# `: N- e* F2 g. v$ V
N3 l+ F* M; r, Q* I
0 E! {# p$ E6 P
- w3 | m# C4 v* f. V
) @# v" s* h! H! Y
- z0 c% n( c- N$ m& K% }/ p7 O
* c' ~/ w6 Q! ?! q- j/ Y W
^6 q; j4 u# W) o8 G$ x也就是说,A\b的方法是可以求包含奇异矩阵的方程组的,但是可能会出错。而且错误可能非常离谱。(这个例子告诫我们,不要以为MATLAB算出来的结果都是准确的,MATLAB也不过是调用一些算法进行运算,每个算法都可能存在一些缺陷,无法处理某些极端的情况)。
- E1 b' \& s) w3 d& O, t% O2 D0 {* _( l9 d" x% W
2 q7 P2 m, H O0 S& U, d
这里就涉及到数值计算方法领域矩阵的性态的问题了。我们可以直观来感受一下:, O2 o2 j5 _0 m" y3 \
: x5 N1 I2 s8 [/ j1 M+ V' j) S& h; U5 G假设如下方程组:) C4 z g9 k# O) U) t; @
. D) Q5 B; E. a5 C0 @" @
$ Y' j4 Y7 f' \% Q" D9 q
' Q" F9 Z) c$ l3 r% }其精确解是(1,1)。
! T& J. A) F" C/ e5 x9 H k" `5 g/ D# Q
若对左右边都做一些非常非常微小的变化:
% i( S$ k* [; J) r8 }' y [& P" Q/ @" Y
; K) _( y" N r0 O b
; Z( e* ]/ u6 x; X+ U/ H% s其精确解变为:(10,-2)。0 m( m; l7 D0 g) v+ @. H
& A' E% C/ ~/ O5 T- M一个非常非常微小的扰动就让方程的解产生巨大的变动,我们称上述方程组是病态方程组,系数矩阵A是病态矩阵。
# ~* [4 ]' A/ b: F) Z' o
7 ~0 B4 f- [1 P0 V如果我们遇到不是方阵的矩阵,或者不能求逆的方阵,要想求解AX=b,避免奇异值导致MATLAB产生错误的情况,我们可以采用“伪逆”来帮助我们解决这个问题。
8 Z# p# r/ E: s- [( h
8 |6 `% w" f. x: g广义逆矩阵:) x3 N/ V1 u v- @1 c3 \; f
: [* u# P, q/ q6 j+ ~3 ^
对任意一个矩阵A,提出四个条件:
, q0 c$ c2 l- a4 B
9 ^4 s/ Q' W- O0 \" n- y$ x
+ K3 N& o$ u2 ]$ b' u( V
. G# d. f$ \+ r# ~2 L8 x* U如果存在矩阵G满足上述的一部分或全部条件,G就可以称为A的广义逆矩阵。最常用的四种广义逆矩阵定义如下:
4 \( i, \/ b' t! V8 o! G/ X6 p, S
* T r6 {8 }1 G0 s" ^, Z" C0 F5 W
9 X9 m% \! O Q+ p+ G8 e$ N6 d1 B3 p* [ T, p7 O
MATLAB中自带的pinv方法,就可以求矩阵的M-P广义逆,即A+矩阵。4 q) r T; {* q
& P$ v8 G& V+ \- e4 ~$ m$ B
大家可以查阅官方文档看具体应用实例。
' U2 Q) ]" h: E/ @. b% w% ?: b3 d! x# g; W
) N; r% C6 }' B, r
9 n, i) o# n3 t( H+ I如果我们求出A+,就可以有另一种思路来解AX=b了:4 Q. w/ h# |6 \. M& t7 [
: {# ~ n% l/ ]
4 b- I" z e3 P! x O u3 a
/ O% L' K! Z4 K" [
这里,通解的表达式还是类似于X=X*+X0的形式,A+b相当于是特解,后面那一项就是带系数的自由解,y可以取任意数,注意维度匹配即可。
' H9 V3 _: q% {
2 _ f$ z. R$ C- _0 i$ Y所以,大家调用pinv求出M-P广义逆,然后用x=A+b + (I-A+A)y这个式子构造出通解就可以啦!
5 `/ [! ?; ]; q" z
2 G7 e! U# I6 D& |+ G0 v4 m' ?; w# J/ D* j/ h0 b1 `
1 c. m) ^& w" n# o3 d- H8 M- A$ q4 g
1 ^$ I- r8 x# p; j) m+ I# b; V# I' m" W! q7 X1 X) d+ c; [
|
|