org.openimaj.math.matrix.ThinSingularValueDecomposition Maven / Gradle / Ivy
/**
* Copyright (c) 2011, The University of Southampton and the individual contributors.
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without modification,
* are permitted provided that the following conditions are met:
*
* * Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* * Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* * Neither the name of the University of Southampton nor the names of its
* contributors may be used to endorse or promote products derived from this
* software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
* ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
* ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
package org.openimaj.math.matrix;
import gov.nist.math.jama.JamaMatrix;
import ch.akuhn.matrix.KuhnMatrix;
import ch.akuhn.matrix.Vector;
import ch.akuhn.matrix.eigenvalues.SingularValues;
/**
* Thin SVD based on Adrian Kuhn's wrapper around ARPACK. This can scale to
* really large matrices (bigger than RAM), given an implementation of
* {@link ch.akuhn.matrix.KuhnMatrix} that is backed by disk.
*
* Note that the current version of (Java) ARPACK is not thread-safe. Allowances
* have been made in this implementation to synchronize the call to
* {@link SingularValues#decompose()} against the
* {@link ThinSingularValueDecomposition} class. Care must be taken if you are
* using Java-ARPACK outside this class in a multi-threaded application.
*
* @author Jonathon Hare ([email protected])
*/
public class ThinSingularValueDecomposition {
/** The U matrix */
public JamaMatrix U;
/** The singular values */
public double[] S;
/** The transpose of the V matrix */
public JamaMatrix Vt;
/**
* Perform thin SVD on matrix, calculating at most ndims dimensions.
*
* @param matrix
* the matrix
* @param ndims
* the number of singular values/vectors to calculate; actual
* number may be less.
*/
public ThinSingularValueDecomposition(JamaMatrix matrix, int ndims) {
this(new JamaDenseMatrix(matrix), ndims);
}
/**
* Perform thin SVD on matrix, calculating at most ndims dimensions.
*
* @param matrix
* the matrix
* @param ndims
* the number of singular values/vectors to calculate; actual
* number may be less.
*/
public ThinSingularValueDecomposition(KuhnMatrix matrix, int ndims) {
// FIXME: I'm (Jon) not sure why this was added, but it causes problems
// with big matrices... commented it out for the time being
// if (ndims > Math.min(matrix.rowCount(), matrix.columnCount())) {
// try {
// final no.uib.cipr.matrix.DenseMatrix mjtA = new
// no.uib.cipr.matrix.DenseMatrix(matrix.asArray());
// no.uib.cipr.matrix.SVD svd;
// svd = no.uib.cipr.matrix.SVD.factorize(mjtA);
// this.S = svd.getS();
// this.U = MatrixUtils.convert(svd.getU());
// this.Vt = MatrixUtils.convert(svd.getVt());
//
// this.S = Arrays.copyOf(this.S, Math.min(ndims, this.S.length));
// this.U = U.getMatrix(0, U.getRowDimension() - 1, 0, Math.min(ndims,
// U.getColumnDimension()) - 1);
// this.Vt = Vt.getMatrix(0, Math.min(Vt.getRowDimension(), ndims) - 1,
// 0, Vt.getColumnDimension() - 1);
// } catch (NotConvergedException e) {
// throw new RuntimeException(e);
// }
// } else {
SingularValues sv = new SingularValues(matrix, ndims);
// Note: SingularValues uses JARPACK which isn't currently
// thread-safe :-(
synchronized (ThinSingularValueDecomposition.class) {
sv.decompose();
}
S = reverse(sv.value);
U = vectorArrayToMatrix(sv.vectorLeft, false);
Vt = vectorArrayToMatrix(sv.vectorRight, true);
// }
}
protected static double[] reverse(double[] vector) {
for (int i = 0; i < vector.length / 2; i++) {
double tmp = vector[i];
vector[i] = vector[vector.length - i - 1];
vector[vector.length - i - 1] = tmp;
}
return vector;
}
protected static JamaMatrix vectorArrayToMatrix(Vector[] vectors, boolean rows) {
final int m = vectors.length;
final double[][] data = new double[m][];
for (int i = 0; i < m; i++) {
data[m - i - 1] = vectors[i].unwrap();
}
JamaMatrix mat = new JamaMatrix(data);
if (!rows) {
mat = mat.transpose();
}
return mat;
}
/**
* @return The S matrix
*/
public JamaMatrix getSmatrix() {
JamaMatrix Smat = new JamaMatrix(S.length, S.length);
for (int r = 0; r < S.length; r++) {
Smat.set(r, r, S[r]);
}
return Smat;
}
/**
* @return The sqrt of the singular vals as a matrix.
*/
public JamaMatrix getSmatrixSqrt() {
JamaMatrix Smat = new JamaMatrix(S.length, S.length);
for (int r = 0; r < S.length; r++) {
Smat.set(r, r, Math.sqrt(S[r]));
}
return Smat;
}
/**
* Reduce the rank of the input matrix using the thin SVD to get a lower
* rank least-squares estimate of the input.
*
* @param m
* matrix to reduce the rank of
* @param rank
* the desired rank
* @return the rank-reduced matrix
*/
public static JamaMatrix reduceRank(JamaMatrix m, int rank) {
if (rank > Math.min(m.getColumnDimension(), m.getRowDimension())) {
return m;
}
ThinSingularValueDecomposition t = new ThinSingularValueDecomposition(m, rank);
return t.U.times(t.getSmatrix()).times(t.Vt);
}
}