logo标签列表

为什么在Matlab和Octave中，inv()和pinv()的输出结果不相等？

matlabprecisionfloating-accuracymatrix-inversenumerical-analysis
15

我注意到如果A是一个NxN矩阵，且它有逆矩阵。但是inv()和pinv()函数的输出结果不同。

以下是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
  • 在Matlab中运行相同的命令，ans结果是相同的。
  • 当我计算inv(A)*A或A*inv(A)时，在Octave和Matlab中都得到了3x3单位矩阵的结果。
  • 在Matlab和FreeMat中，A*pinv(A)和pinv(A)*A的结果都是3x3单位矩阵。
  • 在Octave中，A*pinv(A)的结果是3x3单位矩阵。
  • 在Octave中，pinv(A)*A的结果不是3x3单位矩阵。

我不知道为什么inv(A) != pinv(A)，我考虑了矩阵元素的细节。这似乎是浮点精度问题造成的问题。

小数点后的10+位数字可能会有所不同，如下所示：

  • inv(A)(1,1)中的元素为6.65858991579923298331777914427220821380615200000000
  • pinv(A)(1,1)中的元素为6.65858991579923209513935944414697587490081800000000
可能是为什么Matlab的inv函数又慢又不准确？的重复问题。 - Shai
2
@Shai，我认为OP可能会从你链接的问题的答案中受益（至少如果OP正在使用inv来解决x = A^-1*b），但在我看来这不是重复的。 - Stewie Griffin
4个回答
16
这个问题很旧，但我还是会回答它，因为它在某些谷歌搜索中排名很高。
我将使用magic(N)函数作为我的示例，该函数返回一个N * 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

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
warning: inverse: matrix singular to machine precision, rcond = 1.30614e-17
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

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

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
9

在底部，你似乎已经回答了自己的问题。原因是浮点运算。 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

不要使用 == 来比较矩阵，而应该使用 < tolerance_limit

c = A*pinv(A);
d = pinv(A)*A;

(c-d) < 1e-10

附注:

x = A^-1*b 不应该用 x = inv(A)*b; 来解决，而是要用 x = A \ b;。请参考Shai发布的链接进行解释。

1
我的意思是pinv(A)*A应该是单位矩阵，但现在Octave没有给我这个答案。我知道为什么会出现这个微小的结果，但我们不应该总是得到矩阵乘以它的逆矩阵等于单位矩阵，就像Matlab和Freemat现在正在做的那样吗？ - myme5261314
Octave 给你什么？我猜是一个非常接近单位矩阵的结果（就像我的例子一样）。要理解这是为什么，Daniel 的回答非常详细。Octave 首先计算逆矩阵，然后再与原始矩阵相乘。如果你得到完全不同的结果，请举个例子。 - Stewie Griffin
我得到的ans矩阵与你在答案中的ans矩阵非常相似。我的意思是，作为一个产品或者说作为一个库，Octave应该让pinv(A)*A=单位矩阵，这是有意义的，但现在它返回的ans矩阵没有意义，在Matlab和FreeMat中，pinv(A)*A等于单位矩阵。 - myme5261314
1
@myme5261314 在Octave、Matlab和FreeMat中，逆矩阵的计算方法可能不同，因此答案也会有所不同。如果你将其视为两个单独的计算，先计算A的逆矩阵，然后再乘以A，这就是有意义的。由于第一个计算结果不准确，所以答案也不准确。Octave并不像我们在数学上那样将pinv(A)*A识别为单位矩阵，它实际上进行了计算。 - Stewie Griffin
2

浮点算术有一定的精度，不能依赖于等号。为避免这种错误，您可以尝试使用matlab的符号工具箱来处理。

以下是在Octave中演示此问题的非常简单的代码行：

>>> (1/48)*48==(1/49)*49
ans = 0
>>> (1/48)*48-(1/49)*49
ans =  1.1102e-16
>>>
0
我使用pinv(A)计算了伪逆矩阵A=[1,1,0;1,0,1;2,1,1]（其秩为2）。A*pinv(A)得到一个非单位矩阵，即A*pinv(A)=[0.667, -0.333, 0.333; -0.333, 0.667,0.333; 0.333, 0.333, 0.667]。我认为对于奇异矩阵，最好手动使用svd()计算伪逆矩阵。
这里有一些更新：由于它是非满秩的，所以A*pinv(A)本身可能与单位矩阵不同。这可能与Matlab的pinv(A)算法无关。
矩阵A没有逆矩阵，这意味着你无法找到一个矩阵B，使得AB等于单位矩阵。因此，Apinv(A)不是单位矩阵也就不足为奇了。 - Cris Luengo
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接