001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018 package org.apache.commons.math.linear; 019 020 import org.apache.commons.math.MathRuntimeException; 021 022 /** 023 * Calculates the compact Singular Value Decomposition of a matrix. 024 * <p> 025 * The Singular Value Decomposition of matrix A is a set of three matrices: U, 026 * Σ and V such that A = U × Σ × V<sup>T</sup>. Let A be 027 * a m × n matrix, then U is a m × p orthogonal matrix, Σ is a 028 * p × p diagonal matrix with positive or null elements, V is a p × 029 * n orthogonal matrix (hence V<sup>T</sup> is also orthogonal) where 030 * p=min(m,n). 031 * </p> 032 * @version $Revision: 912413 $ $Date: 2010-02-21 16:46:12 -0500 (Sun, 21 Feb 2010) $ 033 * @since 2.0 034 */ 035 public class SingularValueDecompositionImpl implements 036 SingularValueDecomposition { 037 038 /** Number of rows of the initial matrix. */ 039 private int m; 040 041 /** Number of columns of the initial matrix. */ 042 private int n; 043 044 /** Eigen decomposition of the tridiagonal matrix. */ 045 private EigenDecomposition eigenDecomposition; 046 047 /** Singular values. */ 048 private double[] singularValues; 049 050 /** Cached value of U. */ 051 private RealMatrix cachedU; 052 053 /** Cached value of U<sup>T</sup>. */ 054 private RealMatrix cachedUt; 055 056 /** Cached value of S. */ 057 private RealMatrix cachedS; 058 059 /** Cached value of V. */ 060 private RealMatrix cachedV; 061 062 /** Cached value of V<sup>T</sup>. */ 063 private RealMatrix cachedVt; 064 065 /** 066 * Calculates the compact Singular Value Decomposition of the given matrix. 067 * @param matrix 068 * The matrix to decompose. 069 * @exception InvalidMatrixException 070 * (wrapping a 071 * {@link org.apache.commons.math.ConvergenceException} if 072 * algorithm fails to converge 073 */ 074 public SingularValueDecompositionImpl(final RealMatrix matrix) 075 throws InvalidMatrixException { 076 077 m = matrix.getRowDimension(); 078 n = matrix.getColumnDimension(); 079 080 cachedU = null; 081 cachedS = null; 082 cachedV = null; 083 cachedVt = null; 084 085 double[][] localcopy = matrix.getData(); 086 double[][] matATA = new double[n][n]; 087 // 088 // create A^T*A 089 // 090 for (int i = 0; i < n; i++) { 091 for (int j = i; j < n; j++) { 092 matATA[i][j] = 0.0; 093 for (int k = 0; k < m; k++) { 094 matATA[i][j] += localcopy[k][i] * localcopy[k][j]; 095 } 096 matATA[j][i]=matATA[i][j]; 097 } 098 } 099 100 double[][] matAAT = new double[m][m]; 101 // 102 // create A*A^T 103 // 104 for (int i = 0; i < m; i++) { 105 for (int j = i; j < m; j++) { 106 matAAT[i][j] = 0.0; 107 for (int k = 0; k < n; k++) { 108 matAAT[i][j] += localcopy[i][k] * localcopy[j][k]; 109 } 110 matAAT[j][i]=matAAT[i][j]; 111 } 112 } 113 int p; 114 if (m>=n) { 115 p=n; 116 // compute eigen decomposition of A^T*A 117 eigenDecomposition = new EigenDecompositionImpl( 118 new Array2DRowRealMatrix(matATA),1.0); 119 singularValues = eigenDecomposition.getRealEigenvalues(); 120 cachedV = eigenDecomposition.getV(); 121 122 // compute eigen decomposition of A*A^T 123 eigenDecomposition = new EigenDecompositionImpl( 124 new Array2DRowRealMatrix(matAAT),1.0); 125 cachedU = eigenDecomposition.getV().getSubMatrix(0, m - 1, 0, p - 1); 126 } else { 127 p=m; 128 // compute eigen decomposition of A*A^T 129 eigenDecomposition = new EigenDecompositionImpl( 130 new Array2DRowRealMatrix(matAAT),1.0); 131 singularValues = eigenDecomposition.getRealEigenvalues(); 132 cachedU = eigenDecomposition.getV(); 133 134 // compute eigen decomposition of A^T*A 135 eigenDecomposition = new EigenDecompositionImpl( 136 new Array2DRowRealMatrix(matATA),1.0); 137 cachedV = eigenDecomposition.getV().getSubMatrix(0,n-1,0,p-1); 138 } 139 for (int i = 0; i < p; i++) { 140 singularValues[i] = Math.sqrt(Math.abs(singularValues[i])); 141 } 142 // Up to this point, U and V are computed independently of each other. 143 // There still an sign indetermination of each column of, say, U. 144 // The sign is set such that A.V_i=sigma_i.U_i (i<=p) 145 // The right sign corresponds to a positive dot product of A.V_i and U_i 146 for (int i = 0; i < p; i++) { 147 RealVector tmp = cachedU.getColumnVector(i); 148 double product=matrix.operate(cachedV.getColumnVector(i)).dotProduct(tmp); 149 if (product<0) { 150 cachedU.setColumnVector(i, tmp.mapMultiply(-1.0)); 151 } 152 } 153 } 154 155 /** {@inheritDoc} */ 156 public RealMatrix getU() throws InvalidMatrixException { 157 // return the cached matrix 158 return cachedU; 159 160 } 161 162 /** {@inheritDoc} */ 163 public RealMatrix getUT() throws InvalidMatrixException { 164 165 if (cachedUt == null) { 166 cachedUt = getU().transpose(); 167 } 168 169 // return the cached matrix 170 return cachedUt; 171 172 } 173 174 /** {@inheritDoc} */ 175 public RealMatrix getS() throws InvalidMatrixException { 176 177 if (cachedS == null) { 178 179 // cache the matrix for subsequent calls 180 cachedS = MatrixUtils.createRealDiagonalMatrix(singularValues); 181 182 } 183 return cachedS; 184 } 185 186 /** {@inheritDoc} */ 187 public double[] getSingularValues() throws InvalidMatrixException { 188 return singularValues.clone(); 189 } 190 191 /** {@inheritDoc} */ 192 public RealMatrix getV() throws InvalidMatrixException { 193 // return the cached matrix 194 return cachedV; 195 196 } 197 198 /** {@inheritDoc} */ 199 public RealMatrix getVT() throws InvalidMatrixException { 200 201 if (cachedVt == null) { 202 cachedVt = getV().transpose(); 203 } 204 205 // return the cached matrix 206 return cachedVt; 207 208 } 209 210 /** {@inheritDoc} */ 211 public RealMatrix getCovariance(final double minSingularValue) { 212 213 // get the number of singular values to consider 214 final int p = singularValues.length; 215 int dimension = 0; 216 while ((dimension < p) && (singularValues[dimension] >= minSingularValue)) { 217 ++dimension; 218 } 219 220 if (dimension == 0) { 221 throw MathRuntimeException.createIllegalArgumentException( 222 "cutoff singular value is {0}, should be at most {1}", 223 minSingularValue, singularValues[0]); 224 } 225 226 final double[][] data = new double[dimension][p]; 227 getVT().walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() { 228 /** {@inheritDoc} */ 229 @Override 230 public void visit(final int row, final int column, 231 final double value) { 232 data[row][column] = value / singularValues[row]; 233 } 234 }, 0, dimension - 1, 0, p - 1); 235 236 RealMatrix jv = new Array2DRowRealMatrix(data, false); 237 return jv.transpose().multiply(jv); 238 239 } 240 241 /** {@inheritDoc} */ 242 public double getNorm() throws InvalidMatrixException { 243 return singularValues[0]; 244 } 245 246 /** {@inheritDoc} */ 247 public double getConditionNumber() throws InvalidMatrixException { 248 return singularValues[0] / singularValues[singularValues.length - 1]; 249 } 250 251 /** {@inheritDoc} */ 252 public int getRank() throws IllegalStateException { 253 254 final double threshold = Math.max(m, n) * Math.ulp(singularValues[0]); 255 256 for (int i = singularValues.length - 1; i >= 0; --i) { 257 if (singularValues[i] > threshold) { 258 return i + 1; 259 } 260 } 261 return 0; 262 263 } 264 265 /** {@inheritDoc} */ 266 public DecompositionSolver getSolver() { 267 return new Solver(singularValues, getUT(), getV(), getRank() == Math 268 .max(m, n)); 269 } 270 271 /** Specialized solver. */ 272 private static class Solver implements DecompositionSolver { 273 274 /** Pseudo-inverse of the initial matrix. */ 275 private final RealMatrix pseudoInverse; 276 277 /** Singularity indicator. */ 278 private boolean nonSingular; 279 280 /** 281 * Build a solver from decomposed matrix. 282 * @param singularValues 283 * singularValues 284 * @param uT 285 * U<sup>T</sup> matrix of the decomposition 286 * @param v 287 * V matrix of the decomposition 288 * @param nonSingular 289 * singularity indicator 290 */ 291 private Solver(final double[] singularValues, final RealMatrix uT, 292 final RealMatrix v, final boolean nonSingular) { 293 double[][] suT = uT.getData(); 294 for (int i = 0; i < singularValues.length; ++i) { 295 final double a; 296 if (singularValues[i]>0) { 297 a=1.0 / singularValues[i]; 298 } else { 299 a=0.0; 300 } 301 final double[] suTi = suT[i]; 302 for (int j = 0; j < suTi.length; ++j) { 303 suTi[j] *= a; 304 } 305 } 306 pseudoInverse = v.multiply(new Array2DRowRealMatrix(suT, false)); 307 this.nonSingular = nonSingular; 308 } 309 310 /** 311 * Solve the linear equation A × X = B in least square sense. 312 * <p> 313 * The m×n matrix A may not be square, the solution X is such that 314 * ||A × X - B|| is minimal. 315 * </p> 316 * @param b 317 * right-hand side of the equation A × X = B 318 * @return a vector X that minimizes the two norm of A × X - B 319 * @exception IllegalArgumentException 320 * if matrices dimensions don't match 321 */ 322 public double[] solve(final double[] b) throws IllegalArgumentException { 323 return pseudoInverse.operate(b); 324 } 325 326 /** 327 * Solve the linear equation A × X = B in least square sense. 328 * <p> 329 * The m×n matrix A may not be square, the solution X is such that 330 * ||A × X - B|| is minimal. 331 * </p> 332 * @param b 333 * right-hand side of the equation A × X = B 334 * @return a vector X that minimizes the two norm of A × X - B 335 * @exception IllegalArgumentException 336 * if matrices dimensions don't match 337 */ 338 public RealVector solve(final RealVector b) 339 throws IllegalArgumentException { 340 return pseudoInverse.operate(b); 341 } 342 343 /** 344 * Solve the linear equation A × X = B in least square sense. 345 * <p> 346 * The m×n matrix A may not be square, the solution X is such that 347 * ||A × X - B|| is minimal. 348 * </p> 349 * @param b 350 * right-hand side of the equation A × X = B 351 * @return a matrix X that minimizes the two norm of A × X - B 352 * @exception IllegalArgumentException 353 * if matrices dimensions don't match 354 */ 355 public RealMatrix solve(final RealMatrix b) 356 throws IllegalArgumentException { 357 return pseudoInverse.multiply(b); 358 } 359 360 /** 361 * Check if the decomposed matrix is non-singular. 362 * @return true if the decomposed matrix is non-singular 363 */ 364 public boolean isNonSingular() { 365 return nonSingular; 366 } 367 368 /** 369 * Get the pseudo-inverse of the decomposed matrix. 370 * @return inverse matrix 371 */ 372 public RealMatrix getInverse() { 373 return pseudoInverse; 374 } 375 376 } 377 378 }