Tuesday, April 28, 2009

Moore-Penrose Pseudoinverse in JAMA

import Jama.Matrix;
import Jama.SingularValueDecomposition;

public class Matrices {
 /**
  * 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.
  */
 public static double MACHEPS = 2E-16;

 /**
  * Updates MACHEPS for the executing machine.
  */
 public static void updateMacheps() {
  MACHEPS = 1;
  do
   MACHEPS /= 2;
  while (1 + MACHEPS / 2 != 1);
 }

 /**
  * Computes the Moore–Penrose pseudoinverse using the SVD method.
  * 
  * Modified version of the original implementation by Kim van der Linde.
  */
 public static Matrix pinv(Matrix x) {
  int rows = x.getRowDimension();
  int cols = x.getColumnDimension();
  if (rows < cols) {
   Matrix result = pinv(x.transpose());
   if (result != null)
    result = result.transpose();
   return result;
  }
  SingularValueDecomposition svdX = new SingularValueDecomposition(x);
  if (svdX.rank() < 1)
   return null;
  double[] singularValues = svdX.getSingularValues();
  double tol = Math.max(rows, cols) * singularValues[0] * MACHEPS;
  double[] singularValueReciprocals = new double[singularValues.length];
  for (int i = 0; i < singularValues.length; i++)
   if (Math.abs(singularValues[i]) >= tol)
    singularValueReciprocals[i] =  1.0 / singularValues[i];
  double[][] u = svdX.getU().getArray();
  double[][] v = svdX.getV().getArray();
  int min = Math.min(cols, u[0].length);
  double[][] inverse = new double[cols][rows];
  for (int i = 0; i < cols; i++)
   for (int j = 0; j < u.length; j++)
    for (int k = 0; k < min; k++)
     inverse[i][j] += v[i][k] * singularValueReciprocals[k] * u[j][k];
  return new Matrix(inverse);
 }
}
Update 11/20/2014: As this code continues to live out there, I'm adding a basic test. Just so you know what you're getting. (As is. No warranty. Wish you the best.)
public static boolean checkEquality(Matrix A, Matrix B) {
 return A.minus(B).normInf() < 1e-9;
}
 
public static void testPinv() {
 int rows = (int) Math.floor(100 + Math.random() * 200);
 int cols = (int) Math.floor(100 + Math.random() * 200);
 double[][] mat = new double[rows][cols];
 for (int i = 0; i < rows; i++)
  for (int j = 0; j < cols; j++)
   mat[i][j] = Math.random();
 Matrix A = new Matrix(mat);
 long millis = System.currentTimeMillis();
 Matrix Aplus = pinv(A);
 millis = System.currentTimeMillis() - millis;
 if (Aplus == null) {
  System.out.println("Always check for null");
  return;
 }
 // Good to know.
 boolean c1 = checkEquality(A.times(Aplus).times(A), A);
 boolean c2 = checkEquality(Aplus.times(A).times(Aplus), Aplus);
 boolean c3 = checkEquality(A.times(Aplus), A.times(Aplus).transpose());
 boolean c4 = checkEquality(Aplus.times(A), Aplus.times(A).transpose());
 System.out.println(rows + " x " + cols + "\t" +
                    c1 + "/" + c2 + "/" + c3 + "/" + c4 + "\t" + millis);
}
 
public static void main(String[] args) {
 for (int z = 0; z < 100; z++)
  testPinv();
}

14 comments:

Anonymous said...

Thanks you saved me so much time!!

Anonymous said...

Thanks

Anonymous said...

Add this to the beginning

if (x.getColumnDimension() > x.getRowDimension())
return pinv(x.transpose()).transpose();

Ahmed Abdelkader said...

Done. Thanks!

Anonymous said...

Nice and fast. Thanks!

Anonymous said...

شكراً !!!

Kalle_Mett said...

Very nice, thank you very much!

someone u don know said...
This comment has been removed by the author.
Ismail Marmoush said...

Thanks alot, gazak allah kheer

PrasetiaUtama's said...

Thanks you make my day

Anonymous said...

Thank you so much!!

Unknown said...

Are you aware that x.rank() actually calls new SingularValueDecomposition(this).rank()?
This means that you actually calculate the SVD twice, right?

Ahmed Abdelkader said...

Thanks for pointing that out!

Oleksandr Kovalenko said...

Thank you for saving me the day! I'm grateful to you. The thing is that I'm converting some complex algorithm with pseudoinverse operation in it from Python to Java Interestingly enough, Python Numpy "pinv" uses Moore-Penrose SVD method while JAMA "pinv" employs QR decomposition. And when I load some really big matrices into JAMA "pinv" it throws OutOfMemory error. Why is that, I didn't figure it out yet... For now just using your solution! Thank you again!