mym*_*314 14 precision matlab floating-accuracy matrix-inverse numerical-analysis
我注意到如果A是NxN矩阵并且它具有逆矩阵.但是inv()和pinv()函数输出是不同的. - 我的环境是Win7x64 SP1,Matlab R2012a,Cygwin Octave 3.6.4,FreeMat 4.2
看看Octave的例子:
A = rand(3,3)
A =
0.185987 0.192125 0.046346
0.140710 0.351007 0.236889
0.155899 0.107302 0.300623
pinv(A) == inv(A)
ans =
0 0 0
0 0 0
0 0 0
Run Code Online (Sandbox Code Playgroud)
ans
通过上面运行在Matlab同一命令的结果.inv(A)*A
或A*inv(A)
,结果是Octave和Matlab中的身份3x3矩阵.A*pinv(A)
和pinv(A)*A
在Matlab和FreeMat身份3x3矩阵.A*pinv(A)
是Octave中的身份3x3矩阵.pinv(A)*A
是不是在倍频身份3x3矩阵.我不知道原因 inv(A) != pinv(A)
,我已经考虑了矩阵中元素的细节.这似乎是导致这个问题的浮动精度问题.
点点后10位数可能不同,如下所示:
6.65858991579923298331777914427220821380615200000000
inv(A)(1,1)
反对的元素
6.65858991579923209513935944414697587490081800000000
元素 pinv(A)(1,1)
got*_*oth 15
这个问题很老了,但无论如何我都会回答它,因为它似乎在一些谷歌搜索中几乎排在最前面.
我将在我的例子中使用魔法(N)函数,该函数返回N-by-N魔方.
我将创建一个3x3魔方M3,取伪逆PI_M3并乘以它们:
prompt_$ M3 = magic(3) , PI_M3 = pinv(M3) , M3 * PI_M3
M3 = 8 1 6 3 5 7 4 9 2 PI_M3 = 0.147222 -0.144444 0.063889 -0.061111 0.022222 0.105556 -0.019444 0.188889 -0.102778 ans = 1.0000e+00 -1.2212e-14 6.3283e-15 5.5511e-17 1.0000e+00 -2.2204e-16 -5.9952e-15 1.2268e-14 1.0000e+00
正如您所看到的,答案是单位矩阵,以节省一些舍入误差.我将用4x4魔方重复操作:
prompt_$ M4 = magic(4) , PI_M4 = pinv(M4) , M4 * PI_M4
M4 = 16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1 PI_M4 = 0.1011029 -0.0738971 -0.0613971 0.0636029 -0.0363971 0.0386029 0.0261029 0.0011029 0.0136029 -0.0113971 -0.0238971 0.0511029 -0.0488971 0.0761029 0.0886029 -0.0863971 ans = 0.950000 -0.150000 0.150000 0.050000 -0.150000 0.550000 0.450000 0.150000 0.150000 0.450000 0.550000 -0.150000 0.050000 0.150000 -0.150000 0.950000
结果不是单位矩阵,这意味着4x4魔方不具有逆.我可以通过尝试Moore-Penrose伪逆的一个规则来验证这一点:
prompt_$ M4 * PI_M4 * M4
ans = 16.00000 2.00000 3.00000 13.00000 5.00000 11.00000 10.00000 8.00000 9.00000 7.00000 6.00000 12.00000 4.00000 14.00000 15.00000 1.00000
满足规则A*B*A = A. 这表明pinv在可用时返回逆矩阵,而当逆矩阵不可用时返回伪逆.这就是为什么在某些情况下你会得到一个小的差异,只是一些舍入错误,在其他情况下你会得到更大的差异.为了显示它,我将获得两个魔术象限的倒数并从伪逆中减去它们:
prompt_$ I_M3 = inv(M3) , I_M4 = inv(M4) , DIFF_M3 = PI_M3 - I_M3, DIFF_M4 = PI_M4 - I_M4
I_M3 = 0.147222 -0.144444 0.063889 -0.061111 0.022222 0.105556 -0.019444 0.188889 -0.102778 warning: inverse: matrix singular to machine precision, rcond = 1.30614e-17 I_M4 = 9.3825e+13 2.8147e+14 -2.8147e+14 -9.3825e+13 2.8147e+14 8.4442e+14 -8.4442e+14 -2.8147e+14 -2.8147e+14 -8.4442e+14 8.4442e+14 2.8147e+14 -9.3825e+13 -2.8147e+14 2.8147e+14 9.3825e+13 DIFF_M3 = 4.7184e-16 -1.0270e-15 5.5511e-16 -9.9226e-16 2.0470e-15 -1.0825e-15 5.2042e-16 -1.0270e-15 4.9960e-16 DIFF_M4 = -9.3825e+13 -2.8147e+14 2.8147e+14 9.3825e+13 -2.8147e+14 -8.4442e+14 8.4442e+14 2.8147e+14 2.8147e+14 8.4442e+14 -8.4442e+14 -2.8147e+14 9.3825e+13 2.8147e+14 -2.8147e+14 -9.3825e+13
在我看来,你在底部回答了自己的问题.原因是浮点运算.算法inv()
和pinv()
不完全相同,因为pinv()
必须能够处理非方形矩阵.因此答案将不完全相同.
如果你看一下它的值pinv(A)*A
,你会发现它非常接近单位矩阵.
我明白了:
ans =
1.0000e+00 6.1062e-16 -3.0809e-15
-5.8877e-15 1.0000e+00 6.3942e-15
2.4425e-15 -3.0184e-16 1.0000e+00
Run Code Online (Sandbox Code Playgroud)
而不是比较矩阵==
,使用< tolerance_limit
c = A*pinv(A);
d = pinv(A)*A;
(c-d) < 1e-10
Run Code Online (Sandbox Code Playgroud)
边注:
x = A^-1*b
不应该解决x = inv(A)*b;
,而是x = A \ b;
看到链接Shai发布的解释.
浮点运算有一定的精度,不能依赖等式。为了避免此类错误,您可以尝试使用 matlab 的符号工具箱。
非常简单的八度代码行来演示问题:
>>> (1/48)*48==(1/49)*49
ans = 0
>>> (1/48)*48-(1/49)*49
ans = 1.1102e-16
>>>
Run Code Online (Sandbox Code Playgroud)