/*
* 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.exception.NumberIsTooLargeException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.util.FastMath;
/**
* Calculates the compact Singular Value Decomposition of a matrix.
* <p>
* The Singular Value Decomposition of matrix A is a set of three matrices: U,
* Σ and V such that A = U × Σ × V<sup>T</sup>. Let A be
* a m × n matrix, then U is a m × p orthogonal matrix, Σ is a
* p × p diagonal matrix with positive or null elements, V is a p ×
* n orthogonal matrix (hence V<sup>T</sup> is also orthogonal) where
* p=min(m,n).
* </p>
* @version $Id: SingularValueDecompositionImpl.java 1131229 2011-06-03 20:49:25Z luc $
* @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;
/** Eigen decomposition of the tridiagonal matrix. */
private EigenDecomposition eigenDecomposition;
/** Singular values. */
private double[] singularValues;
/** Cached value of U. */
private RealMatrix cachedU;
/** Cached value of U<sup>T</sup>. */
private RealMatrix cachedUt;
/** Cached value of S. */
private RealMatrix cachedS;
/** Cached value of V. */
private RealMatrix cachedV;
/** Cached value of V<sup>T</sup>. */
private RealMatrix cachedVt;
/**
* Calculates the compact Singular Value Decomposition of the given matrix.
*
* @param matrix Matrix to decompose.
*/
public SingularValueDecompositionImpl(final RealMatrix matrix) {
m = matrix.getRowDimension();
n = matrix.getColumnDimension();
cachedU = null;
cachedS = null;
cachedV = null;
cachedVt = null;
double[][] localcopy = matrix.getData();
double[][] matATA = new double[n][n];
//
// create A^T*A
//
for (int i = 0; i < n; i++) {
for (int j = i; j < n; j++) {
matATA[i][j] = 0.0;
for (int k = 0; k < m; k++) {
matATA[i][j] += localcopy[k][i] * localcopy[k][j];
}
matATA[j][i] = matATA[i][j];
}
}
double[][] matAAT = new double[m][m];
//
// create A*A^T
//
for (int i = 0; i < m; i++) {
for (int j = i; j < m; j++) {
matAAT[i][j] = 0.0;
for (int k = 0; k < n; k++) {
matAAT[i][j] += localcopy[i][k] * localcopy[j][k];
}
matAAT[j][i] = matAAT[i][j];
}
}
int p;
if (m >= n) {
p = n;
// compute eigen decomposition of A^T*A
eigenDecomposition
= new EigenDecompositionImpl(new Array2DRowRealMatrix(matATA), 1);
singularValues = eigenDecomposition.getRealEigenvalues();
cachedV = eigenDecomposition.getV();
// compute eigen decomposition of A*A^T
eigenDecomposition
= new EigenDecompositionImpl(new Array2DRowRealMatrix(matAAT), 1);
cachedU = eigenDecomposition.getV().getSubMatrix(0, m - 1, 0, p - 1);
} else {
p = m;
// compute eigen decomposition of A*A^T
eigenDecomposition
= new EigenDecompositionImpl(new Array2DRowRealMatrix(matAAT), 1);
singularValues = eigenDecomposition.getRealEigenvalues();
cachedU = eigenDecomposition.getV();
// compute eigen decomposition of A^T*A
eigenDecomposition
= new EigenDecompositionImpl(new Array2DRowRealMatrix(matATA), 1);
cachedV = eigenDecomposition.getV().getSubMatrix(0, n - 1 , 0, p - 1);
}
for (int i = 0; i < p; i++) {
singularValues[i] = FastMath.sqrt(FastMath.abs(singularValues[i]));
}
// Up to this point, U and V are computed independently of each other.
// There still a sign indetermination of each column of, say, U.
// The sign is set such that A.V_i=sigma_i.U_i (i<=p)
// The right sign corresponds to a positive dot product of A.V_i and U_i
for (int i = 0; i < p; i++) {
RealVector tmp = cachedU.getColumnVector(i);
double product=matrix.operate(cachedV.getColumnVector(i)).dotProduct(tmp);
if (product < 0) {
cachedU.setColumnVector(i, tmp.mapMultiply(-1));
}
}
}
/** {@inheritDoc} */
public RealMatrix getU() {
// return the cached matrix
return cachedU;
}
/** {@inheritDoc} */
public RealMatrix getUT() {
if (cachedUt == null) {
cachedUt = getU().transpose();
}
// return the cached matrix
return cachedUt;
}
/** {@inheritDoc} */
public RealMatrix getS() {
if (cachedS == null) {
// cache the matrix for subsequent calls
cachedS = MatrixUtils.createRealDiagonalMatrix(singularValues);
}
return cachedS;
}
/** {@inheritDoc} */
public double[] getSingularValues() {
return singularValues.clone();
}
/** {@inheritDoc} */
public RealMatrix getV() {
// return the cached matrix
return cachedV;
}
/** {@inheritDoc} */
public RealMatrix getVT() {
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
final int p = singularValues.length;
int dimension = 0;
while ((dimension < p) && (singularValues[dimension] >= minSingularValue)) {
++dimension;
}
if (dimension == 0) {
throw new NumberIsTooLargeException(LocalizedFormats.TOO_LARGE_CUTOFF_SINGULAR_VALUE,
minSingularValue, singularValues[0], true);
}
final double[][] data = new double[dimension][p];
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, p - 1);
RealMatrix jv = new Array2DRowRealMatrix(data, false);
return jv.transpose().multiply(jv);
}
/** {@inheritDoc} */
public double getNorm() {
return singularValues[0];
}
/** {@inheritDoc} */
public double getConditionNumber() {
return singularValues[0] / singularValues[singularValues.length - 1];
}
/** {@inheritDoc} */
public int getRank() {
final double threshold = FastMath.max(m, n) * FastMath.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() == Math.max(m, n));
}
/** Specialized solver. */
private static class Solver implements DecompositionSolver {
/** Pseudo-inverse of the initial matrix. */
private final RealMatrix pseudoInverse;
/** Singularity indicator. */
private boolean nonSingular;
/**
* Build a solver from decomposed matrix.
*
* @param singularValues Singular values.
* @param uT U<sup>T</sup> 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) {
double[][] suT = uT.getData();
for (int i = 0; i < singularValues.length; ++i) {
final double a;
if (singularValues[i] > 0) {
a = 1 / singularValues[i];
} else {
a = 0;
}
final double[] suTi = suT[i];
for (int j = 0; j < suTi.length; ++j) {
suTi[j] *= a;
}
}
pseudoInverse = v.multiply(new Array2DRowRealMatrix(suT, false));
this.nonSingular = nonSingular;
}
/**
* Solve the linear equation A × X = B in least square sense.
* <p>
* The m×n matrix A may not be square, the solution X is such that
* ||A × X - B|| is minimal.
* </p>
* @param b Right-hand side of the equation A × X = B
* @return a vector X that minimizes the two norm of A × X - B
* @throws org.apache.commons.math.exception.DimensionMismatchException
* if the matrices dimensions do not match.
*/
public double[] solve(final double[] b) {
return pseudoInverse.operate(b);
}
/**
* Solve the linear equation A × X = B in least square sense.
* <p>
* The m×n matrix A may not be square, the solution X is such that
* ||A × X - B|| is minimal.
* </p>
* @param b Right-hand side of the equation A × X = B
* @return a vector X that minimizes the two norm of A × X - B
* @throws org.apache.commons.math.exception.DimensionMismatchException
* if the matrices dimensions do not match.
*/
public RealVector solve(final RealVector b) {
return pseudoInverse.operate(b);
}
/**
* Solve the linear equation A × X = B in least square sense.
* <p>
* The m×n matrix A may not be square, the solution X is such that
* ||A × X - B|| is minimal.
* </p>
*
* @param b Right-hand side of the equation A × X = B
* @return a matrix X that minimizes the two norm of A × X - B
* @throws org.apache.commons.math.exception.DimensionMismatchException
* if the matrices dimensions do not match.
*/
public double[][] solve(final double[][] b) {
return pseudoInverse.multiply(MatrixUtils.createRealMatrix(b)).getData();
}
/**
* Solve the linear equation A × X = B in least square sense.
* <p>
* The m×n matrix A may not be square, the solution X is such that
* ||A × X - B|| is minimal.
* </p>
*
* @param b Right-hand side of the equation A × X = B
* @return a matrix X that minimizes the two norm of A × X - B
* @throws org.apache.commons.math.exception.DimensionMismatchException
* if the matrices dimensions do not match.
*/
public RealMatrix solve(final RealMatrix b) {
return pseudoInverse.multiply(b);
}
/**
* Check if the decomposed matrix is non-singular.
*
* @return {@code true} if the decomposed matrix is non-singular.
*/
public boolean isNonSingular() {
return nonSingular;
}
/**
* Get the pseudo-inverse of the decomposed matrix.
*
* @return the inverse matrix.
*/
public RealMatrix getInverse() {
return pseudoInverse;
}
}
}