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.ConvergenceException; 021 import org.apache.commons.math.MathRuntimeException; 022 import org.apache.commons.math.util.MathUtils; 023 024 /** 025 * Calculates the Singular Value Decomposition of a matrix. 026 * <p>The Singular Value Decomposition of matrix A is a set of three matrices: 027 * U, Σ and V such that A = U × Σ × V<sup>T</sup>. 028 * Let A be an m × n matrix, then U is an m × m orthogonal matrix, 029 * Σ is a m × n diagonal matrix with positive diagonal elements, 030 * and V is an n × n orthogonal matrix.</p> 031 * 032 * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $ 033 * @since 2.0 034 */ 035 public class SingularValueDecompositionImpl implements SingularValueDecomposition { 036 037 /** Number of rows of the initial matrix. */ 038 private int m; 039 040 /** Number of columns of the initial matrix. */ 041 private int n; 042 043 /** Transformer to bidiagonal. */ 044 private BiDiagonalTransformer transformer; 045 046 /** Main diagonal of the bidiagonal matrix. */ 047 private double[] mainBidiagonal; 048 049 /** Secondary diagonal of the bidiagonal matrix. */ 050 private double[] secondaryBidiagonal; 051 052 /** Main diagonal of the tridiagonal matrix. */ 053 private double[] mainTridiagonal; 054 055 /** Secondary diagonal of the tridiagonal matrix. */ 056 private double[] secondaryTridiagonal; 057 058 /** Eigen decomposition of the tridiagonal matrix. */ 059 private EigenDecomposition eigenDecomposition; 060 061 /** Singular values. */ 062 private double[] singularValues; 063 064 /** Cached value of U. */ 065 private RealMatrix cachedU; 066 067 /** Cached value of U<sup>T</sup>. */ 068 private RealMatrix cachedUt; 069 070 /** Cached value of S. */ 071 private RealMatrix cachedS; 072 073 /** Cached value of V. */ 074 private RealMatrix cachedV; 075 076 /** Cached value of V<sup>T</sup>. */ 077 private RealMatrix cachedVt; 078 079 /** 080 * Calculates the Singular Value Decomposition of the given matrix. 081 * @param matrix The matrix to decompose. 082 * @exception InvalidMatrixException (wrapping a {@link ConvergenceException} 083 * if algorithm fails to converge 084 */ 085 public SingularValueDecompositionImpl(RealMatrix matrix) 086 throws InvalidMatrixException { 087 088 m = matrix.getRowDimension(); 089 n = matrix.getColumnDimension(); 090 091 cachedU = null; 092 cachedS = null; 093 cachedV = null; 094 cachedVt = null; 095 096 // transform the matrix to bidiagonal 097 transformer = new BiDiagonalTransformer(matrix); 098 mainBidiagonal = transformer.getMainDiagonalRef(); 099 secondaryBidiagonal = transformer.getSecondaryDiagonalRef(); 100 101 // compute Bt.B (if upper diagonal) or B.Bt (if lower diagonal) 102 mainTridiagonal = new double[mainBidiagonal.length]; 103 secondaryTridiagonal = new double[mainBidiagonal.length - 1]; 104 double a = mainBidiagonal[0]; 105 mainTridiagonal[0] = a * a; 106 for (int i = 1; i < mainBidiagonal.length; ++i) { 107 final double b = secondaryBidiagonal[i - 1]; 108 secondaryTridiagonal[i - 1] = a * b; 109 a = mainBidiagonal[i]; 110 mainTridiagonal[i] = a * a + b * b; 111 } 112 113 // compute singular values 114 eigenDecomposition = 115 new EigenDecompositionImpl(mainTridiagonal, secondaryTridiagonal, 116 MathUtils.SAFE_MIN); 117 singularValues = eigenDecomposition.getRealEigenvalues(); 118 for (int i = 0; i < singularValues.length; ++i) { 119 singularValues[i] = Math.sqrt(singularValues[i]); 120 } 121 122 } 123 124 /** {@inheritDoc} */ 125 public RealMatrix getU() 126 throws InvalidMatrixException { 127 128 if (cachedU == null) { 129 130 if (m >= n) { 131 // the tridiagonal matrix is Bt.B, where B is upper bidiagonal 132 final double[][] eData = eigenDecomposition.getV().getData(); 133 final double[][] iData = new double[m][]; 134 double[] ei1 = eData[0]; 135 iData[0] = ei1; 136 for (int i = 0; i < n - 1; ++i) { 137 // compute Bt.E.S^(-1) where E is the eigenvectors matrix 138 // we reuse the array from matrix E to store the result 139 final double[] ei0 = ei1; 140 ei1 = eData[i + 1]; 141 iData[i + 1] = ei1; 142 for (int j = 0; j < n; ++j) { 143 ei0[j] = (mainBidiagonal[i] * ei0[j] + 144 secondaryBidiagonal[i] * ei1[j]) / singularValues[j]; 145 } 146 } 147 // last row 148 final double lastMain = mainBidiagonal[n - 1]; 149 for (int j = 0; j < n; ++j) { 150 ei1[j] *= lastMain / singularValues[j]; 151 } 152 for (int i = n; i < m; ++i) { 153 iData[i] = new double[n]; 154 } 155 cachedU = 156 transformer.getU().multiply(MatrixUtils.createRealMatrix(iData)); 157 } else { 158 // the tridiagonal matrix is B.Bt, where B is lower bidiagonal 159 cachedU = transformer.getU().multiply(eigenDecomposition.getV()); 160 } 161 162 } 163 164 // return the cached matrix 165 return cachedU; 166 167 } 168 169 /** {@inheritDoc} */ 170 public RealMatrix getUT() 171 throws InvalidMatrixException { 172 173 if (cachedUt == null) { 174 cachedUt = getU().transpose(); 175 } 176 177 // return the cached matrix 178 return cachedUt; 179 180 } 181 182 /** {@inheritDoc} */ 183 public RealMatrix getS() 184 throws InvalidMatrixException { 185 186 if (cachedS == null) { 187 188 // cache the matrix for subsequent calls 189 cachedS = MatrixUtils.createRealDiagonalMatrix(singularValues); 190 191 } 192 return cachedS; 193 } 194 195 /** {@inheritDoc} */ 196 public double[] getSingularValues() 197 throws InvalidMatrixException { 198 return singularValues.clone(); 199 } 200 201 /** {@inheritDoc} */ 202 public RealMatrix getV() 203 throws InvalidMatrixException { 204 205 if (cachedV == null) { 206 207 if (m >= n) { 208 // the tridiagonal matrix is Bt.B, where B is upper bidiagonal 209 cachedV = transformer.getV().multiply(eigenDecomposition.getV()); 210 } else { 211 // the tridiagonal matrix is B.Bt, where B is lower bidiagonal 212 final double[][] eData = eigenDecomposition.getV().getData(); 213 final double[][] iData = new double[n][]; 214 double[] ei1 = eData[0]; 215 iData[0] = ei1; 216 for (int i = 0; i < m - 1; ++i) { 217 // compute Bt.E.S^(-1) where E is the eigenvectors matrix 218 // we reuse the array from matrix E to store the result 219 final double[] ei0 = ei1; 220 ei1 = eData[i + 1]; 221 iData[i + 1] = ei1; 222 for (int j = 0; j < m; ++j) { 223 ei0[j] = (mainBidiagonal[i] * ei0[j] + 224 secondaryBidiagonal[i] * ei1[j]) / singularValues[j]; 225 } 226 } 227 // last row 228 final double lastMain = mainBidiagonal[m - 1]; 229 for (int j = 0; j < m; ++j) { 230 ei1[j] *= lastMain / singularValues[j]; 231 } 232 for (int i = m; i < n; ++i) { 233 iData[i] = new double[m]; 234 } 235 cachedV = 236 transformer.getV().multiply(MatrixUtils.createRealMatrix(iData)); 237 } 238 239 } 240 241 // return the cached matrix 242 return cachedV; 243 244 } 245 246 /** {@inheritDoc} */ 247 public RealMatrix getVT() 248 throws InvalidMatrixException { 249 250 if (cachedVt == null) { 251 cachedVt = getV().transpose(); 252 } 253 254 // return the cached matrix 255 return cachedVt; 256 257 } 258 259 /** {@inheritDoc} */ 260 public RealMatrix getCovariance(final double minSingularValue) { 261 262 // get the number of singular values to consider 263 int dimension = 0; 264 while ((dimension < n) && (singularValues[dimension] >= minSingularValue)) { 265 ++dimension; 266 } 267 268 if (dimension == 0) { 269 throw MathRuntimeException.createIllegalArgumentException( 270 "cutoff singular value is {0}, should be at most {1}", 271 minSingularValue, singularValues[0]); 272 } 273 274 final double[][] data = new double[dimension][n]; 275 getVT().walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() { 276 /** {@inheritDoc} */ 277 @Override 278 public void visit(final int row, final int column, final double value) { 279 data[row][column] = value / singularValues[row]; 280 } 281 }, 0, dimension - 1, 0, n - 1); 282 283 RealMatrix jv = new Array2DRowRealMatrix(data, false); 284 return jv.transpose().multiply(jv); 285 286 } 287 288 /** {@inheritDoc} */ 289 public double getNorm() 290 throws InvalidMatrixException { 291 return singularValues[0]; 292 } 293 294 /** {@inheritDoc} */ 295 public double getConditionNumber() 296 throws InvalidMatrixException { 297 return singularValues[0] / singularValues[singularValues.length - 1]; 298 } 299 300 /** {@inheritDoc} */ 301 public int getRank() 302 throws IllegalStateException { 303 304 final double threshold = Math.max(m, n) * Math.ulp(singularValues[0]); 305 306 for (int i = singularValues.length - 1; i >= 0; --i) { 307 if (singularValues[i] > threshold) { 308 return i + 1; 309 } 310 } 311 return 0; 312 313 } 314 315 /** {@inheritDoc} */ 316 public DecompositionSolver getSolver() { 317 return new Solver(singularValues, getUT(), getV(), 318 getRank() == singularValues.length); 319 } 320 321 /** Specialized solver. */ 322 private static class Solver implements DecompositionSolver { 323 324 /** Singular values. */ 325 private final double[] singularValues; 326 327 /** U<sup>T</sup> matrix of the decomposition. */ 328 private final RealMatrix uT; 329 330 /** V matrix of the decomposition. */ 331 private final RealMatrix v; 332 333 /** Singularity indicator. */ 334 private boolean nonSingular; 335 336 /** 337 * Build a solver from decomposed matrix. 338 * @param singularValues singularValues 339 * @param uT U<sup>T</sup> matrix of the decomposition 340 * @param v V matrix of the decomposition 341 * @param nonSingular singularity indicator 342 */ 343 private Solver(final double[] singularValues, final RealMatrix uT, final RealMatrix v, 344 final boolean nonSingular) { 345 this.singularValues = singularValues; 346 this.uT = uT; 347 this.v = v; 348 this.nonSingular = nonSingular; 349 } 350 351 /** Solve the linear equation A × X = B in least square sense. 352 * <p>The m×n matrix A may not be square, the solution X is 353 * such that ||A × X - B|| is minimal.</p> 354 * @param b right-hand side of the equation A × X = B 355 * @return a vector X that minimizes the two norm of A × X - B 356 * @exception IllegalArgumentException if matrices dimensions don't match 357 * @exception InvalidMatrixException if decomposed matrix is singular 358 */ 359 public double[] solve(final double[] b) 360 throws IllegalArgumentException, InvalidMatrixException { 361 362 if (b.length != singularValues.length) { 363 throw MathRuntimeException.createIllegalArgumentException( 364 "vector length mismatch: got {0} but expected {1}", 365 b.length, singularValues.length); 366 } 367 368 final double[] w = uT.operate(b); 369 for (int i = 0; i < singularValues.length; ++i) { 370 final double si = singularValues[i]; 371 if (si == 0) { 372 throw new SingularMatrixException(); 373 } 374 w[i] /= si; 375 } 376 return v.operate(w); 377 378 } 379 380 /** Solve the linear equation A × X = B in least square sense. 381 * <p>The m×n matrix A may not be square, the solution X is 382 * such that ||A × X - B|| is minimal.</p> 383 * @param b right-hand side of the equation A × X = B 384 * @return a vector X that minimizes the two norm of A × X - B 385 * @exception IllegalArgumentException if matrices dimensions don't match 386 * @exception InvalidMatrixException if decomposed matrix is singular 387 */ 388 public RealVector solve(final RealVector b) 389 throws IllegalArgumentException, InvalidMatrixException { 390 391 if (b.getDimension() != singularValues.length) { 392 throw MathRuntimeException.createIllegalArgumentException( 393 "vector length mismatch: got {0} but expected {1}", 394 b.getDimension(), singularValues.length); 395 } 396 397 final RealVector w = uT.operate(b); 398 for (int i = 0; i < singularValues.length; ++i) { 399 final double si = singularValues[i]; 400 if (si == 0) { 401 throw new SingularMatrixException(); 402 } 403 w.setEntry(i, w.getEntry(i) / si); 404 } 405 return v.operate(w); 406 407 } 408 409 /** Solve the linear equation A × X = B in least square sense. 410 * <p>The m×n matrix A may not be square, the solution X is 411 * such that ||A × X - B|| is minimal.</p> 412 * @param b right-hand side of the equation A × X = B 413 * @return a matrix X that minimizes the two norm of A × X - B 414 * @exception IllegalArgumentException if matrices dimensions don't match 415 * @exception InvalidMatrixException if decomposed matrix is singular 416 */ 417 public RealMatrix solve(final RealMatrix b) 418 throws IllegalArgumentException, InvalidMatrixException { 419 420 if (b.getRowDimension() != singularValues.length) { 421 throw MathRuntimeException.createIllegalArgumentException( 422 "dimensions mismatch: got {0}x{1} but expected {2}x{3}", 423 b.getRowDimension(), b.getColumnDimension(), 424 singularValues.length, "n"); 425 } 426 427 final RealMatrix w = uT.multiply(b); 428 for (int i = 0; i < singularValues.length; ++i) { 429 final double si = singularValues[i]; 430 if (si == 0) { 431 throw new SingularMatrixException(); 432 } 433 final double inv = 1.0 / si; 434 for (int j = 0; j < b.getColumnDimension(); ++j) { 435 w.multiplyEntry(i, j, inv); 436 } 437 } 438 return v.multiply(w); 439 440 } 441 442 /** 443 * Check if the decomposed matrix is non-singular. 444 * @return true if the decomposed matrix is non-singular 445 */ 446 public boolean isNonSingular() { 447 return nonSingular; 448 } 449 450 /** Get the pseudo-inverse of the decomposed matrix. 451 * @return inverse matrix 452 * @throws InvalidMatrixException if decomposed matrix is singular 453 */ 454 public RealMatrix getInverse() 455 throws InvalidMatrixException { 456 457 if (!isNonSingular()) { 458 throw new SingularMatrixException(); 459 } 460 461 return solve(MatrixUtils.createRealIdentityMatrix(singularValues.length)); 462 463 } 464 465 } 466 467 }