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

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

[复制链接]

该用户从未签到

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

EDA365欢迎您登录!

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

x
上一篇讲到:方程AX=b的解的讨论(特解、通解、零空间向量等概念)及其MATLAB实现,程序中用到的是mldivide或者A\b的方法(二者相同)来解方程。& t, E1 A9 H: R3 X  A9 ~
" }5 I; s: K$ W& V7 i7 U
但实际上运行过程中我们会遇到:当AX=b线性方程组是一个病态方程组;或者A是奇异矩阵(即det(A)=0,不可逆),没法求逆,用不了inv(A)方法只能用A\b,此时MATLAB会报错“矩阵接近奇异值,或者缩放错误。结果可能不准确”…网络上很多人问这个问题怎么解决,其实不是MATLAB的问题,而是MATLAB内置算法的鲁棒性问题,直接用A\b方法无法处理这个棘手的问题。如果没有矩阵论或数值计算方法基础的同学可能会一头雾水。本文借用Moore-Penrose广义逆来解决这个问题,帮助大家理解带奇异矩阵的病态方程组如何解决。
% {& I" i/ F/ `% w
8 U$ y, R0 V, Z, {% x首先我们先来看下mldivide, \在MATLAB中的含义:) y' D& k) m( S
; q' Z1 _1 U: s! x. n. D7 T

5 {0 G4 y# |0 e: \2 n
% S* o6 B6 m# Q; W" @
! K1 p8 V; H" T% U
0 i4 I4 \5 ?% @+ u5 u2 X+ L$ j+ i1 t' {

7 R& x0 F: O$ V1 T也就是说,A\b的方法是可以求包含奇异矩阵的方程组的,但是可能会出错。而且错误可能非常离谱。(这个例子告诫我们,不要以为MATLAB算出来的结果都是准确的,MATLAB也不过是调用一些算法进行运算,每个算法都可能存在一些缺陷,无法处理某些极端的情况)。, e* A- ~% m* J8 |' v* z' G

! u% s* D. q& X8 ]" e  Y5 D) ?. {( q1 o
这里就涉及到数值计算方法领域矩阵的性态的问题了。我们可以直观来感受一下:
- A4 D" k- N% m& @( k7 {) W
$ h& p# R) ?3 p2 u! x6 _6 c假设如下方程组:. i( q' [& `) ?6 `, M: o; d
( i% S% o2 z9 B4 o) X* R
( O9 C  Y! q% J& V6 ^7 X: \
0 {% d; y. g+ {' P
其精确解是(1,1)。" Y% K  L, e8 q. j

0 w, ~% P) s7 D/ R- t, i若对左右边都做一些非常非常微小的变化:
2 R$ G* y! L$ `# @/ F  z3 ^% C. e
: ~. r3 C( Q- `) {4 P

1 Y, g  y2 k7 V; N( Y其精确解变为:(10,-2)。
4 t+ g3 v( Z% Q" C6 O$ }& c- H) J6 y0 O
一个非常非常微小的扰动就让方程的解产生巨大的变动,我们称上述方程组是病态方程组,系数矩阵A是病态矩阵。
: ?* m# ~5 D% r8 A% Z
* U( Y7 l5 E$ r* r4 N9 A# `* d. B  q如果我们遇到不是方阵的矩阵,或者不能求逆的方阵,要想求解AX=b,避免奇异值导致MATLAB产生错误的情况,我们可以采用“伪逆”来帮助我们解决这个问题。& s7 e% @3 P6 L: w3 b& e0 N
' V4 a, N4 l0 [7 K4 b
广义逆矩阵:) ?/ k6 P3 G3 m& W

3 P- M1 \7 a; ~1 D对任意一个矩阵A,提出四个条件:( E9 [2 L$ A1 G4 W
5 q9 O- u2 F) H8 e$ A* Y$ s! P( q

) M- z' d7 I) r: q3 C  X' q" T2 ~, e0 q* T
如果存在矩阵G满足上述的一部分或全部条件,G就可以称为A的广义逆矩阵。最常用的四种广义逆矩阵定义如下:( j( V- G7 U  d- s& }, o

$ J+ q1 v" F/ h: A3 R3 { 0 a3 t: I% }6 }- \/ t5 K
6 T8 ]6 g/ H& M% f" ~. U3 x' _
MATLAB中自带的pinv方法,就可以求矩阵的M-P广义逆,即A+矩阵。0 Y8 D, g3 v/ }. {. ~
1 `  H7 I# ?8 ~2 z$ D. U9 P  R
大家可以查阅官方文档看具体应用实例。
  X2 B) ]+ z$ E/ K: r! {5 L0 Z; U9 y6 s$ G: z8 f! M& J& P

  e; h3 k6 s; U' m7 j( q- R1 E/ {! K6 o
$ N7 h$ ^8 f0 E3 ?# [: n; @如果我们求出A+,就可以有另一种思路来解AX=b了:
% V- {* E$ L6 [# D% O# c! D' N8 ~1 f. K  v) O, Q
0 r+ H+ u2 h" H+ {3 B2 d

2 S' v& }& o  ?5 k' E, g5 G这里,通解的表达式还是类似于X=X*+X0的形式,A+b相当于是特解,后面那一项就是带系数的自由解,y可以取任意数,注意维度匹配即可。
* c' A9 l. B# G5 h: l! c8 x+ p9 M
所以,大家调用pinv求出M-P广义逆,然后用x=A+b + (I-A+A)y这个式子构造出通解就可以啦!# ^3 f* S- j! e. H# w0 ^' k# X2 w

2 f! S$ Z- |# z# y3 W# z- z
2 v$ U5 v) h9 {+ X3 z$ r! d' R! e$ t9 F
! B# Z0 X4 b: }0 c

  J! ]  w( a4 L: [
  • 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-11-24 15:39 , Processed in 0.203125 second(s), 26 queries , Gzip On.

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

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

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