org.apache.commons.math.linear.SingularValueDecompositionImpl Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of commons-math Show documentation
Show all versions of commons-math Show documentation
The Math project is a library of lightweight, self-contained mathematics and statistics components addressing the most common practical problems not immediately available in the Java programming language or commons-lang.
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.linear;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.util.MathUtils;
/**
* Calculates the Singular Value Decomposition of a matrix.
* The Singular Value Decomposition of matrix A is a set of three matrices:
* U, Σ and V such that A = U × Σ × VT.
* Let A be an m × n matrix, then U is an m × m orthogonal matrix,
* Σ is a m × n diagonal matrix with positive diagonal elements,
* and V is an n × n orthogonal matrix.
*
* @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
* @since 2.0
*/
public class SingularValueDecompositionImpl implements SingularValueDecomposition {
/** Number of rows of the initial matrix. */
private int m;
/** Number of columns of the initial matrix. */
private int n;
/** Transformer to bidiagonal. */
private BiDiagonalTransformer transformer;
/** Main diagonal of the bidiagonal matrix. */
private double[] mainBidiagonal;
/** Secondary diagonal of the bidiagonal matrix. */
private double[] secondaryBidiagonal;
/** Main diagonal of the tridiagonal matrix. */
private double[] mainTridiagonal;
/** Secondary diagonal of the tridiagonal matrix. */
private double[] secondaryTridiagonal;
/** Eigen decomposition of the tridiagonal matrix. */
private EigenDecomposition eigenDecomposition;
/** Singular values. */
private double[] singularValues;
/** Cached value of U. */
private RealMatrix cachedU;
/** Cached value of UT. */
private RealMatrix cachedUt;
/** Cached value of S. */
private RealMatrix cachedS;
/** Cached value of V. */
private RealMatrix cachedV;
/** Cached value of VT. */
private RealMatrix cachedVt;
/**
* Calculates the Singular Value Decomposition of the given matrix.
* @param matrix The matrix to decompose.
* @exception InvalidMatrixException (wrapping a {@link ConvergenceException}
* if algorithm fails to converge
*/
public SingularValueDecompositionImpl(RealMatrix matrix)
throws InvalidMatrixException {
m = matrix.getRowDimension();
n = matrix.getColumnDimension();
cachedU = null;
cachedS = null;
cachedV = null;
cachedVt = null;
// transform the matrix to bidiagonal
transformer = new BiDiagonalTransformer(matrix);
mainBidiagonal = transformer.getMainDiagonalRef();
secondaryBidiagonal = transformer.getSecondaryDiagonalRef();
// compute Bt.B (if upper diagonal) or B.Bt (if lower diagonal)
mainTridiagonal = new double[mainBidiagonal.length];
secondaryTridiagonal = new double[mainBidiagonal.length - 1];
double a = mainBidiagonal[0];
mainTridiagonal[0] = a * a;
for (int i = 1; i < mainBidiagonal.length; ++i) {
final double b = secondaryBidiagonal[i - 1];
secondaryTridiagonal[i - 1] = a * b;
a = mainBidiagonal[i];
mainTridiagonal[i] = a * a + b * b;
}
// compute singular values
eigenDecomposition =
new EigenDecompositionImpl(mainTridiagonal, secondaryTridiagonal,
MathUtils.SAFE_MIN);
singularValues = eigenDecomposition.getRealEigenvalues();
for (int i = 0; i < singularValues.length; ++i) {
singularValues[i] = Math.sqrt(singularValues[i]);
}
}
/** {@inheritDoc} */
public RealMatrix getU()
throws InvalidMatrixException {
if (cachedU == null) {
if (m >= n) {
// the tridiagonal matrix is Bt.B, where B is upper bidiagonal
final double[][] eData = eigenDecomposition.getV().getData();
final double[][] iData = new double[m][];
double[] ei1 = eData[0];
iData[0] = ei1;
for (int i = 0; i < n - 1; ++i) {
// compute Bt.E.S^(-1) where E is the eigenvectors matrix
// we reuse the array from matrix E to store the result
final double[] ei0 = ei1;
ei1 = eData[i + 1];
iData[i + 1] = ei1;
for (int j = 0; j < n; ++j) {
ei0[j] = (mainBidiagonal[i] * ei0[j] +
secondaryBidiagonal[i] * ei1[j]) / singularValues[j];
}
}
// last row
final double lastMain = mainBidiagonal[n - 1];
for (int j = 0; j < n; ++j) {
ei1[j] *= lastMain / singularValues[j];
}
for (int i = n; i < m; ++i) {
iData[i] = new double[n];
}
cachedU =
transformer.getU().multiply(MatrixUtils.createRealMatrix(iData));
} else {
// the tridiagonal matrix is B.Bt, where B is lower bidiagonal
cachedU = transformer.getU().multiply(eigenDecomposition.getV());
}
}
// return the cached matrix
return cachedU;
}
/** {@inheritDoc} */
public RealMatrix getUT()
throws InvalidMatrixException {
if (cachedUt == null) {
cachedUt = getU().transpose();
}
// return the cached matrix
return cachedUt;
}
/** {@inheritDoc} */
public RealMatrix getS()
throws InvalidMatrixException {
if (cachedS == null) {
// cache the matrix for subsequent calls
cachedS = MatrixUtils.createRealDiagonalMatrix(singularValues);
}
return cachedS;
}
/** {@inheritDoc} */
public double[] getSingularValues()
throws InvalidMatrixException {
return singularValues.clone();
}
/** {@inheritDoc} */
public RealMatrix getV()
throws InvalidMatrixException {
if (cachedV == null) {
if (m >= n) {
// the tridiagonal matrix is Bt.B, where B is upper bidiagonal
cachedV = transformer.getV().multiply(eigenDecomposition.getV());
} else {
// the tridiagonal matrix is B.Bt, where B is lower bidiagonal
final double[][] eData = eigenDecomposition.getV().getData();
final double[][] iData = new double[n][];
double[] ei1 = eData[0];
iData[0] = ei1;
for (int i = 0; i < m - 1; ++i) {
// compute Bt.E.S^(-1) where E is the eigenvectors matrix
// we reuse the array from matrix E to store the result
final double[] ei0 = ei1;
ei1 = eData[i + 1];
iData[i + 1] = ei1;
for (int j = 0; j < m; ++j) {
ei0[j] = (mainBidiagonal[i] * ei0[j] +
secondaryBidiagonal[i] * ei1[j]) / singularValues[j];
}
}
// last row
final double lastMain = mainBidiagonal[m - 1];
for (int j = 0; j < m; ++j) {
ei1[j] *= lastMain / singularValues[j];
}
for (int i = m; i < n; ++i) {
iData[i] = new double[m];
}
cachedV =
transformer.getV().multiply(MatrixUtils.createRealMatrix(iData));
}
}
// return the cached matrix
return cachedV;
}
/** {@inheritDoc} */
public RealMatrix getVT()
throws InvalidMatrixException {
if (cachedVt == null) {
cachedVt = getV().transpose();
}
// return the cached matrix
return cachedVt;
}
/** {@inheritDoc} */
public RealMatrix getCovariance(final double minSingularValue) {
// get the number of singular values to consider
int dimension = 0;
while ((dimension < n) && (singularValues[dimension] >= minSingularValue)) {
++dimension;
}
if (dimension == 0) {
throw MathRuntimeException.createIllegalArgumentException(
"cutoff singular value is {0}, should be at most {1}",
minSingularValue, singularValues[0]);
}
final double[][] data = new double[dimension][n];
getVT().walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
/** {@inheritDoc} */
@Override
public void visit(final int row, final int column, final double value) {
data[row][column] = value / singularValues[row];
}
}, 0, dimension - 1, 0, n - 1);
RealMatrix jv = new Array2DRowRealMatrix(data, false);
return jv.transpose().multiply(jv);
}
/** {@inheritDoc} */
public double getNorm()
throws InvalidMatrixException {
return singularValues[0];
}
/** {@inheritDoc} */
public double getConditionNumber()
throws InvalidMatrixException {
return singularValues[0] / singularValues[singularValues.length - 1];
}
/** {@inheritDoc} */
public int getRank()
throws IllegalStateException {
final double threshold = Math.max(m, n) * Math.ulp(singularValues[0]);
for (int i = singularValues.length - 1; i >= 0; --i) {
if (singularValues[i] > threshold) {
return i + 1;
}
}
return 0;
}
/** {@inheritDoc} */
public DecompositionSolver getSolver() {
return new Solver(singularValues, getUT(), getV(),
getRank() == singularValues.length);
}
/** Specialized solver. */
private static class Solver implements DecompositionSolver {
/** Singular values. */
private final double[] singularValues;
/** UT matrix of the decomposition. */
private final RealMatrix uT;
/** V matrix of the decomposition. */
private final RealMatrix v;
/** Singularity indicator. */
private boolean nonSingular;
/**
* Build a solver from decomposed matrix.
* @param singularValues singularValues
* @param uT UT matrix of the decomposition
* @param v V matrix of the decomposition
* @param nonSingular singularity indicator
*/
private Solver(final double[] singularValues, final RealMatrix uT, final RealMatrix v,
final boolean nonSingular) {
this.singularValues = singularValues;
this.uT = uT;
this.v = v;
this.nonSingular = nonSingular;
}
/** Solve the linear equation A × X = B in least square sense.
* The m×n matrix A may not be square, the solution X is
* such that ||A × X - B|| is minimal.
* @param b right-hand side of the equation A × X = B
* @return a vector X that minimizes the two norm of A × X - B
* @exception IllegalArgumentException if matrices dimensions don't match
* @exception InvalidMatrixException if decomposed matrix is singular
*/
public double[] solve(final double[] b)
throws IllegalArgumentException, InvalidMatrixException {
if (b.length != singularValues.length) {
throw MathRuntimeException.createIllegalArgumentException(
"vector length mismatch: got {0} but expected {1}",
b.length, singularValues.length);
}
final double[] w = uT.operate(b);
for (int i = 0; i < singularValues.length; ++i) {
final double si = singularValues[i];
if (si == 0) {
throw new SingularMatrixException();
}
w[i] /= si;
}
return v.operate(w);
}
/** Solve the linear equation A × X = B in least square sense.
* The m×n matrix A may not be square, the solution X is
* such that ||A × X - B|| is minimal.
* @param b right-hand side of the equation A × X = B
* @return a vector X that minimizes the two norm of A × X - B
* @exception IllegalArgumentException if matrices dimensions don't match
* @exception InvalidMatrixException if decomposed matrix is singular
*/
public RealVector solve(final RealVector b)
throws IllegalArgumentException, InvalidMatrixException {
if (b.getDimension() != singularValues.length) {
throw MathRuntimeException.createIllegalArgumentException(
"vector length mismatch: got {0} but expected {1}",
b.getDimension(), singularValues.length);
}
final RealVector w = uT.operate(b);
for (int i = 0; i < singularValues.length; ++i) {
final double si = singularValues[i];
if (si == 0) {
throw new SingularMatrixException();
}
w.setEntry(i, w.getEntry(i) / si);
}
return v.operate(w);
}
/** Solve the linear equation A × X = B in least square sense.
* The m×n matrix A may not be square, the solution X is
* such that ||A × X - B|| is minimal.
* @param b right-hand side of the equation A × X = B
* @return a matrix X that minimizes the two norm of A × X - B
* @exception IllegalArgumentException if matrices dimensions don't match
* @exception InvalidMatrixException if decomposed matrix is singular
*/
public RealMatrix solve(final RealMatrix b)
throws IllegalArgumentException, InvalidMatrixException {
if (b.getRowDimension() != singularValues.length) {
throw MathRuntimeException.createIllegalArgumentException(
"dimensions mismatch: got {0}x{1} but expected {2}x{3}",
b.getRowDimension(), b.getColumnDimension(),
singularValues.length, "n");
}
final RealMatrix w = uT.multiply(b);
for (int i = 0; i < singularValues.length; ++i) {
final double si = singularValues[i];
if (si == 0) {
throw new SingularMatrixException();
}
final double inv = 1.0 / si;
for (int j = 0; j < b.getColumnDimension(); ++j) {
w.multiplyEntry(i, j, inv);
}
}
return v.multiply(w);
}
/**
* Check if the decomposed matrix is non-singular.
* @return true if the decomposed matrix is non-singular
*/
public boolean isNonSingular() {
return nonSingular;
}
/** Get the pseudo-inverse of the decomposed matrix.
* @return inverse matrix
* @throws InvalidMatrixException if decomposed matrix is singular
*/
public RealMatrix getInverse()
throws InvalidMatrixException {
if (!isNonSingular()) {
throw new SingularMatrixException();
}
return solve(MatrixUtils.createRealIdentityMatrix(singularValues.length));
}
}
}
© 2015 - 2024 Weber Informatics LLC | Privacy Policy