| /* |
| * 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.MathRuntimeException; |
| 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 $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $ |
| * @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 |
| * The matrix to decompose. |
| * @exception InvalidMatrixException |
| * (wrapping a |
| * {@link org.apache.commons.math.ConvergenceException} if |
| * algorithm fails to converge |
| */ |
| public SingularValueDecompositionImpl(final RealMatrix matrix) |
| throws InvalidMatrixException { |
| |
| 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.0); |
| singularValues = eigenDecomposition.getRealEigenvalues(); |
| cachedV = eigenDecomposition.getV(); |
| // compute eigen decomposition of A*A^T |
| eigenDecomposition = new EigenDecompositionImpl( |
| new Array2DRowRealMatrix(matAAT),1.0); |
| 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.0); |
| singularValues = eigenDecomposition.getRealEigenvalues(); |
| cachedU = eigenDecomposition.getV(); |
| |
| // compute eigen decomposition of A^T*A |
| eigenDecomposition = new EigenDecompositionImpl( |
| new Array2DRowRealMatrix(matATA),1.0); |
| 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.0)); |
| } |
| } |
| } |
| |
| /** {@inheritDoc} */ |
| public RealMatrix getU() throws InvalidMatrixException { |
| // 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 { |
| // 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 |
| final int p = singularValues.length; |
| int dimension = 0; |
| while ((dimension < p) && (singularValues[dimension] >= minSingularValue)) { |
| ++dimension; |
| } |
| |
| if (dimension == 0) { |
| throw MathRuntimeException.createIllegalArgumentException( |
| LocalizedFormats.TOO_LARGE_CUTOFF_SINGULAR_VALUE, |
| minSingularValue, singularValues[0]); |
| } |
| |
| 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() 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 = 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 |
| * singularValues |
| * @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.0 / singularValues[i]; |
| } else { |
| a=0.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 |
| * @exception IllegalArgumentException |
| * if matrices dimensions don't match |
| */ |
| public double[] solve(final double[] b) throws IllegalArgumentException { |
| 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 |
| * @exception IllegalArgumentException |
| * if matrices dimensions don't match |
| */ |
| public RealVector solve(final RealVector b) |
| throws IllegalArgumentException { |
| 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 |
| * @exception IllegalArgumentException |
| * if matrices dimensions don't match |
| */ |
| public RealMatrix solve(final RealMatrix b) |
| throws IllegalArgumentException { |
| return pseudoInverse.multiply(b); |
| } |
| |
| /** |
| * 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 |
| */ |
| public RealMatrix getInverse() { |
| return pseudoInverse; |
| } |
| |
| } |
| |
| } |
