Matlab: Moore-Penrose pseudo inverse algorithm implementation

Question

I am searching for a Matlab implementation of the Moore-Penrose algorithm computing pseudo-inverse matrix.

I tried several algoithm, this one

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

appeared good at the first look.

However, the problem it, that for large elements it produces badly scaled matrices and some internal operations fail. It concerns the following steps:

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

I am trying to find a more robust solution which is easily implementable in my other SW. Thanks for your help.

You may want to take this question to either http://math.stackexchange.com/ or http://dsp.stackexchange.com/ — Ali, Nov 11 '12 at 10:54

score 4 · Answer 1 · answered Nov 11 '12 at 11:19

4

What is wrong with using the built-in pinv?

Otherwise, you could take a look at the implementation used in Octave. It is not in Octave/MATLAB syntax, but I guess you should be able to port it without much problems.

answered Nov 11 '12 at 11:19

Egon

4,757
1
23
38

Yes, pinv is OK. But I would like to use the code in another SW written in different language. – justik Nov 11 '12 at 11:28
I think the pseudo-inverse should be available for almost any decent programming language (e.g. using the LAPACK library). In general I wouldn't recommend implementing numerical algorithms yourself for anything that should be reliable (unless of course you know what you are doing). – Egon Nov 11 '12 at 22:41

score 4 · Accepted Answer · answered Dec 08 '12 at 15:22

I re-implemented one in C# using the Mapack matrix library by Lutz Roeder. Perhaps this, or the Java version, will be useful to you.

/// <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;
}

score 1 · Answer 3 · answered Sep 26 '14 at 18:11

Here is the R code [I][1] have written to compute M-P pseudoinverse. I think that is simple enough to be translated into matlab code.

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

Matlab: Moore-Penrose pseudo inverse algorithm implementation

3 Answers3