手动计算伪逆矩阵

Question

我按照维基百科的公式：http://en.wikipedia.org/wiki/Pseudo_inverse，计算伪逆矩阵但是我没有得到正确的结果。例如： 我想要找出方程Y=R*theta中的theta，我在matlab上写了以下代码：

R = -[1/sqrt(2) 1 1/sqrt(2) 0;0 1/sqrt(2) 1 1/sqrt(2);-1/sqrt(2) 0 1/sqrt(2) 1];
% R is 3x4 matrix

Y = [0; -1/sqrt(2);-1]; %Y is 3x1 matrix

B1 = pinv(R);
theta1 = B1*Y;
result1 = R*theta1 - Y

B2 = R'*inv(R*R');
theta2 = B2*Y;
result2 = R*theta2 - Y

这是结果：

   result1 =
   1.0e-15 *
   -0.1110
   -0.2220
   -0.2220
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND =  1.904842e-17. 
> In pseudoinverse at 14 
result2 =
    0.1036
   -0.1768
   -0.3536

显然，theta2是错误的答案，但我不知道为什么。我读了很多书，它们给我提供了与维基百科相同的公式。 有人能帮我手动计算伪逆吗？谢谢！

- kerry_13

---

我觉得你的意思是"theta1 是错误的答案"，就像你写的那样。 - horchler

这是theta2。theta1正确，因为result1约等于0，而theta2错误，因为result2不等于0。 - kerry_13

这不是 theta2。你没有叫做 theta2 的变量。你在脚本中所称的 theta2 实际上被命名为 theta1。我只是想为其他人澄清你的问题。我无法编辑你的问题来更改一个字符 - 你需要自己去做。 - horchler

是的，我看到了，谢谢 ^^ - kerry_13

Answer 1

代数学告诉你，伪逆可以用来解决这样的方程，但是代数并没有考虑有限精度计算的影响。

事实上，使用矩阵乘法方法计算伪逆是不可行的，因为它在数值上是不稳定的。使用\运算符进行矩阵除法。

theta = R \ Y;

从代数角度看，矩阵除法和伪逆相乘是相同的。但MATLAB的实现更加稳定。

有关更多信息，包括稳定的方法，请参见

• 维基百科上的Moore-Penrose伪逆

- Ben Voigt

---

谢谢您的回答，但在我的问题中，pinv() 显示了正确的答案，而公式 B1 = R'*inv(R*R'); 是错误的。我想知道 pinv 的实现方式。 我尝试使用 theta = R \ B;，但 Matlab 给出了一个错误提示： ??? Error using ==> mldivide Matrix dimensions must agree. - kerry_13

@kerry_13：请查看我提供的链接。不清楚为什么你还没有看到它，因为它在你回答中链接了多次。 - Ben Voigt

抱歉打错了。你使用的是B伪逆，而我习惯于用Ax = b来描述方程系统，然后解决方案是x = A \ b;。 - Ben Voigt

在Matlab命令窗口中输入“edit pinv”。您会发现该函数基于“ svd”函数（奇异值分解） -正如@BenVoigt提供的维基百科链接所示。 - horchler

Answer 2

正如Ben所说，矩阵求逆是数值不稳定的。函数inv不建议使用，除非你想要实际求解矩阵的逆，例如请参考this link。误用inv是新学数值线性代数学生最常犯的错误。

在大多数线性代数计算中，可以通过使用数值稳定的算法来避免使用inv。例如，我们有LU分解用于线性求解器，以及QR或SVD方法用于普通最小二乘。

- Da Kuang

Answer 3

你没有正确计算B1。在你的情况下

B1 = inv(R'*R)*R';

因为R是领先的（传统上是相反的）。然而，这并不能解决您的奇异性问题。

pinv使用SVD来计算伪逆，您可以在此处阅读有关它的信息。

因此，基本上它以更优雅的方式执行：

[U,S,V] = svd(R);
S(abs(S)<(sum(sum(S))*1e-8)) = 0; % removes near-singular values.
Stemp = 1./S;
Stemp(isinf(Stemp)) = 0; % This take the inverse of non-zero terms... I'm sure there is faster way
B1 = V * Stemp' * U';

然后你就可以继续你的路程...

- Rasman