|
|
EDA365欢迎您登录！
您需要 登录 才可以下载或查看，没有帐号？注册
x
上一篇讲到：方程AX=b的解的讨论（特解、通解、零空间向量等概念）及其MATLAB实现，程序中用到的是mldivide或者A\b的方法（二者相同）来解方程。( T9 b5 x' _& y% m
# K# h' I( ]* I4 R1 c
但实际上运行过程中我们会遇到：当AX=b线性方程组是一个病态方程组；或者A是奇异矩阵（即det(A)=0，不可逆），没法求逆，用不了inv(A)方法只能用A\b，此时MATLAB会报错“矩阵接近奇异值，或者缩放错误。结果可能不准确”…网络上很多人问这个问题怎么解决，其实不是MATLAB的问题，而是MATLAB内置算法的鲁棒性问题，直接用A\b方法无法处理这个棘手的问题。如果没有矩阵论或数值计算方法基础的同学可能会一头雾水。本文借用Moore-Penrose广义逆来解决这个问题，帮助大家理解带奇异矩阵的病态方程组如何解决。, w& D0 a! G& H. S9 S+ z6 n
+ B4 N( O# @0 [; w3 g( Q
首先我们先来看下mldivide, \在MATLAB中的含义：* }$ c) A0 z" T! l ^+ ~
$ J/ R( {6 V9 v8 L9 E
" i, S2 Z8 Z8 c1 p0 `
( Z* h1 T! t/ }" I4 B6 E
6 i1 Y- O1 e( M4 T, L4 I# G) ?- l8 a; j9 G+ I, O j
7 r p& j! E- g& |( Z( U3 i# n5 E. E
也就是说，A\b的方法是可以求包含奇异矩阵的方程组的，但是可能会出错。而且错误可能非常离谱。（这个例子告诫我们，不要以为MATLAB算出来的结果都是准确的，MATLAB也不过是调用一些算法进行运算，每个算法都可能存在一些缺陷，无法处理某些极端的情况）。
3 T' n! S3 I9 M( v8 Z4 L
* ^$ e2 O' g- h# v3 h
! I4 E7 z* d: R _1 l( Q这里就涉及到数值计算方法领域矩阵的性态的问题了。我们可以直观来感受一下：8 ]: j7 I" h9 S8 e0 F7 B
5 k" n1 M* a8 {( O/ n# p0 b
假设如下方程组： f; d; ?4 M @# M% C F
+ y$ j/ Q/ }8 I, N% p/ o! ~% F2 I
- `" }0 N8 y' _0 A1 ?- \
# k4 j7 n0 ?) Q8 D8 I; ^其精确解是(1,1)。! C8 P% C# H" r0 D
; D/ F" z& N& `9 O, |1 U( U
若对左右边都做一些非常非常微小的变化：
& |. }, h& T' z/ I, j( ]
3 k: V' y6 E& v; l# g7 ?
- C/ g3 h* v$ l! A5 e- |
! v7 o) R3 a# [1 r, A: a
其精确解变为：(10,-2)。
6 k1 g! L' i; w9 ^6 c0 i i. R0 v$ T4 g9 `' ^/ X' _; M) s
一个非常非常微小的扰动就让方程的解产生巨大的变动，我们称上述方程组是病态方程组，系数矩阵A是病态矩阵。
Y1 Z2 z! i @. F) _1 y) M
4 L, g' D8 B' \. N/ j# e/ o如果我们遇到不是方阵的矩阵，或者不能求逆的方阵，要想求解AX=b，避免奇异值导致MATLAB产生错误的情况，我们可以采用“伪逆”来帮助我们解决这个问题。
: a/ t0 N' D" k f( \6 H* b. N. D) I+ B
广义逆矩阵：
* N/ |- v4 @" k# ~ t4 u
- z, \6 Q; J3 [- h对任意一个矩阵A，提出四个条件：5 L0 }' u4 e' F
2 ?9 P$ _3 X) Q+ q7 W) j+ i
; C, k4 f% M2 \ e& Z" H# }0 Z2 F1 L2 @
如果存在矩阵G满足上述的一部分或全部条件，G就可以称为A的广义逆矩阵。最常用的四种广义逆矩阵定义如下：
: D o7 }2 o1 e+ I# P: w7 w: `$ y' }" K
/ {2 U* T9 g# ~" p: H
0 g/ k, {* W8 [- @$ { `MATLAB中自带的pinv方法，就可以求矩阵的M-P广义逆，即A+矩阵。0 s5 O0 v6 @3 O) z8 q2 w
7 u) u$ ?) `2 q1 k+ ]$ E k
大家可以查阅官方文档看具体应用实例。' T2 p% \' n6 C0 B' `
1 l O) C. c$ z
& ~0 w" P/ K$ _: Q6 [& g- T
V+ r1 h2 m0 z- B6 }0 `
如果我们求出A+，就可以有另一种思路来解AX=b了：
. U7 N' K7 Z5 j1 L# `
! Y H3 I/ X6 v8 M0 N2 B
: [. i5 G3 [ Z+ M- _* L( S
4 X# D8 i- Y; X+ y) H: ^& A这里，通解的表达式还是类似于X=X*+X0的形式，A+b相当于是特解，后面那一项就是带系数的自由解，y可以取任意数，注意维度匹配即可。. H4 Z/ [9 I4 e6 P2 y
/ a$ l+ Q' {( K( q* H4 z6 o- f; y所以，大家调用pinv求出M-P广义逆，然后用x=A+b + (I-A+A)y这个式子构造出通解就可以啦！9 Z& v* s% _2 v& [2 N& [
, j5 H. U! M' {
# A1 x) q: {7 q! E9 T
5 w- b2 ~9 t' s6 e
4 a9 _7 I1 s; G! @0 w4 Z/ @
" y1 O% T3 k, c/ N |
|
/1 
发表于 2020-6-5 15:17
收藏
淘帖
支持！
反对！