找回密码
 注册
关于网站域名变更的通知
查看: 3261|回复: 2
打印 上一主题 下一主题

Moore-Penrose广义逆:可解决MATLAB报错“矩阵接近奇异值,或者缩放错误。结果可能...

[复制链接]

该用户从未签到

跳转到指定楼层
1#
发表于 2020-6-5 15:17 | 只看该作者 |只看大图 回帖奖励 |倒序浏览 |阅读模式

EDA365欢迎您登录!

您需要 登录 才可以下载或查看,没有帐号?注册

x
上一篇讲到:方程AX=b的解的讨论(特解、通解、零空间向量等概念)及其MATLAB实现,程序中用到的是mldivide或者A\b的方法(二者相同)来解方程。/ }, Z8 a2 I! C

2 |* u% M- E3 b3 L8 P但实际上运行过程中我们会遇到:当AX=b线性方程组是一个病态方程组;或者A是奇异矩阵(即det(A)=0,不可逆),没法求逆,用不了inv(A)方法只能用A\b,此时MATLAB会报错“矩阵接近奇异值,或者缩放错误。结果可能不准确”…网络上很多人问这个问题怎么解决,其实不是MATLAB的问题,而是MATLAB内置算法的鲁棒性问题,直接用A\b方法无法处理这个棘手的问题。如果没有矩阵论或数值计算方法基础的同学可能会一头雾水。本文借用Moore-Penrose广义逆来解决这个问题,帮助大家理解带奇异矩阵的病态方程组如何解决。, c2 ~  I, d, V; P# V: k% T# n# b

) I; o- @6 D( ~: D* U% x6 d首先我们先来看下mldivide, \在MATLAB中的含义:  S2 }& L* f8 L9 F$ c
- u  v7 _  O" a4 j1 V6 j
, ?+ x, v, E: F* Q6 o- m  u% e

, A: c9 `( L7 v  Z/ X0 E5 D . z/ Z) K* z5 @  t

3 k: ?* y1 g, z* V0 S7 ]( L3 ~  U  s& T; q3 u$ O
: p0 d/ c6 ~- b
也就是说,A\b的方法是可以求包含奇异矩阵的方程组的,但是可能会出错。而且错误可能非常离谱。(这个例子告诫我们,不要以为MATLAB算出来的结果都是准确的,MATLAB也不过是调用一些算法进行运算,每个算法都可能存在一些缺陷,无法处理某些极端的情况)。* }" o' G4 c" x% b

/ v; O5 L$ I' W+ s3 O# O  B& F
& g( D0 ~7 V( k8 Z这里就涉及到数值计算方法领域矩阵的性态的问题了。我们可以直观来感受一下:
( N$ M) X6 F' V9 f0 T* N, `
9 _& P& I( J' M' y) x( m5 u假设如下方程组:
2 l* H( A; X/ `* o5 N- y/ S4 G7 ]' S

: q0 c6 w% t! r: B8 _( Y
% y" M, N7 o3 i6 _* Q9 h其精确解是(1,1)。
: r/ D* ]& ^6 p4 e6 _' q6 c* M: y
; ?8 ~+ x, q$ w9 F2 w# Z2 O9 @若对左右边都做一些非常非常微小的变化:
* x: C" _7 y& F( H7 ]7 W3 S0 ?
% V6 P+ h" d8 B8 q- M" I
- H' `- L1 P# O
1 {6 d( s/ @" w1 ^) y6 d其精确解变为:(10,-2)。
! T1 d  R% r* P5 ~7 ]" E
% Y2 v" W" i4 h一个非常非常微小的扰动就让方程的解产生巨大的变动,我们称上述方程组是病态方程组,系数矩阵A是病态矩阵。& S' z1 a/ ^. B1 M

" K9 k8 g# N+ f2 {如果我们遇到不是方阵的矩阵,或者不能求逆的方阵,要想求解AX=b,避免奇异值导致MATLAB产生错误的情况,我们可以采用“伪逆”来帮助我们解决这个问题。
5 H% J5 o, L+ F* A3 K0 f( h! Y' `6 ?6 Y9 k- k% S
广义逆矩阵:
/ a0 n2 e0 z1 H; x  u' \5 g- ^
/ [! v4 k, t) ]7 g% h对任意一个矩阵A,提出四个条件:
' ~+ o* u4 D1 z- J" [$ |7 D
! i6 \- m* k  b 7 h4 ^) q7 G' U* P0 ?7 T
  t  |' i5 z- C- h$ f7 x
如果存在矩阵G满足上述的一部分或全部条件,G就可以称为A的广义逆矩阵。最常用的四种广义逆矩阵定义如下:, `/ d0 W1 t5 A5 K1 z
5 _+ p4 x/ {3 r( {+ h& \4 W

# ~9 ?/ N3 O- l
( Z5 k0 ^' g% p3 x8 FMATLAB中自带的pinv方法,就可以求矩阵的M-P广义逆,即A+矩阵。$ S. T; f1 q/ u, [4 A; e

& R/ a4 a. {  H- m大家可以查阅官方文档看具体应用实例。- ~- l5 E' M: m9 W& K5 V
( M- E% m. p6 `5 ]+ v2 T0 _( Q
. V: [  O- u, |8 j* X

: ~% E) L" Z8 B& x- V0 z  k4 b如果我们求出A+,就可以有另一种思路来解AX=b了:
0 U! d. O$ c7 O. s0 J6 e. U3 C! C# ?" G5 {, J: c; Q

4 a+ D7 z- Z$ ], Q' w( K0 W- r& h7 f1 W* b
这里,通解的表达式还是类似于X=X*+X0的形式,A+b相当于是特解,后面那一项就是带系数的自由解,y可以取任意数,注意维度匹配即可。
* h$ j6 F6 b6 F6 `6 J* W4 D" B# a; ^1 B* v8 @
所以,大家调用pinv求出M-P广义逆,然后用x=A+b + (I-A+A)y这个式子构造出通解就可以啦!
6 p0 `* [2 o2 m; k# y) e/ X; r
; }8 H+ v, |$ H5 f) D: b, W' I) d
) k2 O# i6 Z6 p

' B5 R! n3 m: v. z6 M$ G4 L* ?  y! a% U1 n+ z: ]& d
  • TA的每日心情

    2019-11-29 15:37
  • 签到天数: 1 天

    [LV.1]初来乍到

    3#
    发表于 2020-6-8 17:03 | 只看该作者
    可解决MATLAB报错“矩阵接近奇异值,或者缩放错误
    您需要登录后才可以回帖 登录 | 注册

    本版积分规则

    关闭

    推荐内容上一条 /1 下一条

    EDA365公众号

    关于我们|手机版|EDA365电子论坛网 ( 粤ICP备18020198号-1 )

    GMT+8, 2025-7-24 05:31 , Processed in 0.125000 second(s), 26 queries , Gzip On.

    深圳市墨知创新科技有限公司

    地址:深圳市南山区科技生态园2栋A座805 电话:19926409050

    快速回复 返回顶部 返回列表