http://openimaj.org/apidocs/src-html/org/openimaj/math/matrix/PseudoInverse.html

//Copyright (c) 2011, The University of Southampton and the individual contributors.

package org.openimaj.math.matrix;
import no.uib.cipr.matrix.NotConvergedException;
import Jama.Matrix;
/**
* Methods for calculating the Moore-Penrose Pseudo-Inverse
* @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
*/
public class PseudoInverse {
/**
* Compute the Moore-Penrose Pseudo-Inverse.
*@param matrix: The matrix to invert.
* @return the pseudo-inverse.
*/
public static Matrix pseudoInverse(Matrix matrix) {
final no.uib.cipr.matrix.DenseMatrix mjtA =
new no.uib.cipr.matrix.DenseMatrix(matrix.getArray());
no.uib.cipr.matrix.SVD svd;

try {
svd = no.uib.cipr.matrix.SVD.factorize(mjtA);
} catch (final NotConvergedException e) {
throw new RuntimeException(e);
}

final Matrix Sinv = new Matrix(matrix.getColumnDimension(),
matrix.getRowDimension());
final double[] Sarr = svd.getS();
for (int i = 0; i < svd.getS().length; i++) {
if (Sarr[i] != 0)
Sinv.set(i, i, 1.0 / Sarr[i]);
}

final Matrix Vt = new Matrix(svd.getVt().numRows(),
svd.getVt().numColumns());
for (int r = 0; r < svd.getVt().numRows(); r++) {
for (int c = 0; c < svd.getVt().numColumns(); c++) {
Vt.set(r, c, svd.getVt().get(r, c));
}
}

final Matrix U = new Matrix(svd.getU().numRows(),
svd.getU().numColumns());
for (int r = 0; r < svd.getU().numRows(); r++) {
for (int c = 0; c < svd.getU().numColumns(); c++) {
U.set(r, c, svd.getU().get(r, c));
}
}

final Matrix pinv = Vt.transpose().times(Sinv).times(U.transpose());
return pinv;
}

/**
* Compute the lower-rank approximation of the Moore-Penrose Pseudo-Inverse.
*  @param matrix: The matrix to invert.
* @param rank:the desired rank.
* @return the pseudo-inverse.
*/

public static Matrix pseudoInverse(Matrix matrix, int rank) {
return pseudoInverse(new JamaDenseMatrix(matrix), rank);
}

/**
* Compute the lower-rank approximation of the Moore-Penrose Pseudo-Inverse.
* @param matrix: The matrix to invert.
* @param rank: the desired rank.
* @return the pseudo-inverse.
*/
public static Matrix pseudoInverse(ch.akuhn.matrix.Matrix matrix,int rank) {
final ThinSingularValueDecomposition tsvd =
new ThinSingularValueDecomposition(matrix, rank);

final Matrix Sinv = new Matrix(tsvd.S.length, tsvd.S.length);
for (int i = 0; i < tsvd.S.length; i++) {
if (tsvd.S[i] != 0)
Sinv.set(i, i, 1.0 / tsvd.S[i]);
}

final Matrix pinv =                                                                                                                                                                           tsvd.Vt.transpose().times(Sinv).times(tsvd.U.transpose());                                return pinv;
}
}

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);
 }
}
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. wiki
 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();
}
Advertisements