Matlab：Moore-Penrose伪逆算法实现答案

【问题标题】：Matlab: Moore-Penrose pseudo inverse algorithm implementationMatlab：Moore-Penrose伪逆算法实现
【发布时间】：2012-11-11 10:52:08
【问题描述】：

我正在寻找一个用于计算伪逆矩阵的 Moore-Penrose 算法的 Matlab 实现。

我尝试了几种算法，这个

http://arxiv.org/ftp/arxiv/papers/0804/0804.4809.pdf

乍一看还不错。

但是，问题在于，对于大型元素，它会产生缩放不良的矩阵，并且某些内部操作会失败。它涉及以下步骤：

L=L(:,1:r);
M=inv(L'*L);

我正在尝试找到一个更强大的解决方案，该解决方案可以在我的其他软件中轻松实施。感谢您的帮助。

【问题讨论】：

您可以将此问题发送给math.stackexchange.com 或dsp.stackexchange.com

标签： matlab matrix-inverse

【解决方案1】：

我使用 Lutz Roeder 的 Mapack 矩阵库在 C# 中重新实现了一个。也许这个或 Java 版本对您有用。

/// <summary>
/// The difference between 1 and the smallest exactly representable number
/// greater than one. Gives an upper bound on the relative error due to
/// rounding of floating point numbers.
/// </summary>
const double MACHEPS = 2E-16;

// NOTE: Code for pseudoinverse is from:
// http://the-lost-beauty.blogspot.com/2009/04/moore-penrose-pseudoinverse-in-jama.html

/// <summary>
/// Computes the Moore–Penrose pseudoinverse using the SVD method.
/// Modified version of the original implementation by Kim van der Linde.
/// </summary>
/// <param name="x"></param>
/// <returns>The pseudoinverse.</returns>
public static Matrix MoorePenrosePsuedoinverse(Matrix x)
{
    if (x.Columns > x.Rows)
        return MoorePenrosePsuedoinverse(x.Transpose()).Transpose();
    SingularValueDecomposition svdX = new SingularValueDecomposition(x);
    if (svdX.Rank < 1)
        return null;
    double[] singularValues = svdX.Diagonal;
    double tol = Math.Max(x.Columns, x.Rows) * singularValues[0] * MACHEPS;
    double[] singularValueReciprocals = new double[singularValues.Length];
    for (int i = 0; i < singularValues.Length; ++i)
        singularValueReciprocals[i] = Math.Abs(singularValues[i]) < tol ? 0 : (1.0 / singularValues[i]);
    Matrix u = svdX.GetU();
    Matrix v = svdX.GetV();
    int min = Math.Min(x.Columns, u.Columns);
    Matrix inverse = new Matrix(x.Columns, x.Rows);
    for (int i = 0; i < x.Columns; i++)
        for (int j = 0; j < u.Rows; j++)
            for (int k = 0; k < min; k++)
                inverse[i, j] += v[i, k] * singularValueReciprocals[k] * u[j, k];
    return inverse;
}

【讨论】：

【解决方案2】：

使用内置pinv有什么问题？

否则，您可以查看implementation used in Octave。它不是 Octave/MATLAB 语法，但我想你应该能够毫无问题地移植它。

【讨论】：

是的，pinv 没问题。但我想在另一个用不同语言编写的软件中使用代码。
我认为伪逆应该适用于几乎任何体面的编程语言（例如使用 LAPACK 库）。一般来说，我不建议自己为任何应该可靠的东西实现数值算法（当然，除非你知道自己在做什么）。

【解决方案3】：

这是 [I][1] 编写的用于计算 M-P 伪逆的 R 代码。我认为这很简单，可以翻译成matlab代码。

pinv<-function(H){
  x=t(H) %*% H
  s=svd(x)
  xp=s$d
  for (i in 1:length(xp)){
    if (xp[i] != 0){
      xp[i]=1/xp[i]
    }
    else{
      xp[i]=0
    }
  }
  return(s$u %*% diag(xp) %*% t(s$v) %*% t(H))
}

[1]：http://hamedhaseli.webs.com/downloads

【讨论】：