为什么 inv() 和 pinv() 在 Matlab 和 Octave 中的输出不相等？答案

【问题标题】：Why is the output of inv() and pinv() not equal in Matlab and Octave?为什么 inv() 和 pinv() 在 Matlab 和 Octave 中的输出不相等？
【发布时间】：2013-10-25 17:38:38
【问题描述】：

我注意到如果 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

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 的结果不是 Octave 中的身份 3x3 矩阵。

不知道inv(A) != pinv(A)的原因，我考虑过矩阵中元素的细节。似乎是导致此问题的浮动精度问题。

点后的10+位可能不同如下：

inv(A)(1,1) 中的6.65858991579923298331777914427220821380615200000000 元素反对
6.65858991579923209513935944414697587490081800000000 元素在pinv(A)(1,1)

【问题讨论】：

Why is Matlab's inv slow and inaccurate?的可能重复
@Shai，我相信 OP 可能会从阅读您链接到的问题的答案中受益（至少如果 OP 使用inv 解决x = A^-1*b），但IMO 这不是重复的.

标签： matlab precision floating-accuracy matrix-inverse numerical-analysis

【解决方案1】：

这个问题已经很老了，但我还是会回答它，因为它在一些谷歌搜索中几乎出现在首位。

我将在我的示例中使用返回 N×N 魔方的 magic(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 答案= 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 答案= 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 伪逆规则之一来验证这一点：

 提示_$ M4 * PI_M4 * M4

答案= 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 警告：逆：矩阵奇异到机器精度，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

【讨论】：

【解决方案2】：

在我看来，您似乎在底部回答了您自己的问题。原因是浮点运算。 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; 参见the link Shai posted 的解释。

【讨论】：

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

【解决方案3】：

浮点运算有一定的精度，不能依赖相等。为避免此类错误，您可以尝试使用 matlab 的符号工具箱。

非常简单的八度代码行来演示问题：

>>> (1/48)*48==(1/49)*49
ans = 0
>>> (1/48)*48-(1/49)*49
ans =  1.1102e-16
>>>

【讨论】：

【解决方案4】：

我使用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) 本身可能与单位矩阵不同，因为它是非满秩的。可能不依赖于pinv(A)的Matlab算法。

【讨论】：

A 没有逆矩阵，这意味着您找不到矩阵B 使得A*B 是单位矩阵。所以难怪A*pinv(A)不是单位矩阵。