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 import org.apache.commons.math.MaxIterationsExceededException; 022 import org.apache.commons.math.util.MathUtils; 023 024 /** 025 * Calculates the eigen decomposition of a real <strong>symmetric</strong> 026 * matrix. 027 * <p> 028 * The eigen decomposition of matrix A is a set of two matrices: V and D such 029 * that A = V D V<sup>T</sup>. A, V and D are all m × m matrices. 030 * </p> 031 * <p> 032 * As of 2.0, this class supports only <strong>symmetric</strong> matrices, and 033 * hence computes only real realEigenvalues. This implies the D matrix returned 034 * by {@link #getD()} is always diagonal and the imaginary values returned 035 * {@link #getImagEigenvalue(int)} and {@link #getImagEigenvalues()} are always 036 * null. 037 * </p> 038 * <p> 039 * When called with a {@link RealMatrix} argument, this implementation only uses 040 * the upper part of the matrix, the part below the diagonal is not accessed at 041 * all. 042 * </p> 043 * <p> 044 * This implementation is based on the paper by A. Drubrulle, R.S. Martin and 045 * J.H. Wilkinson 'The Implicit QL Algorithm' in Wilksinson and Reinsch (1971) 046 * Handbook for automatic computation, vol. 2, Linear algebra, Springer-Verlag, 047 * New-York 048 * </p> 049 * @version $Revision: 912413 $ $Date: 2010-02-21 16:46:12 -0500 (Sun, 21 Feb 2010) $ 050 * @since 2.0 051 */ 052 public class EigenDecompositionImpl implements EigenDecomposition { 053 054 /** Maximum number of iterations accepted in the implicit QL transformation */ 055 private byte maxIter = 30; 056 057 /** Main diagonal of the tridiagonal matrix. */ 058 private double[] main; 059 060 /** Secondary diagonal of the tridiagonal matrix. */ 061 private double[] secondary; 062 063 /** 064 * Transformer to tridiagonal (may be null if matrix is already 065 * tridiagonal). 066 */ 067 private TriDiagonalTransformer transformer; 068 069 /** Real part of the realEigenvalues. */ 070 private double[] realEigenvalues; 071 072 /** Imaginary part of the realEigenvalues. */ 073 private double[] imagEigenvalues; 074 075 /** Eigenvectors. */ 076 private ArrayRealVector[] eigenvectors; 077 078 /** Cached value of V. */ 079 private RealMatrix cachedV; 080 081 /** Cached value of D. */ 082 private RealMatrix cachedD; 083 084 /** Cached value of Vt. */ 085 private RealMatrix cachedVt; 086 087 /** 088 * Calculates the eigen decomposition of the given symmetric matrix. 089 * @param matrix The <strong>symmetric</strong> matrix to decompose. 090 * @param splitTolerance dummy parameter, present for backward compatibility only. 091 * @exception InvalidMatrixException (wrapping a 092 * {@link org.apache.commons.math.ConvergenceException} if algorithm 093 * fails to converge 094 */ 095 public EigenDecompositionImpl(final RealMatrix matrix,final double splitTolerance) 096 throws InvalidMatrixException { 097 if (isSymmetric(matrix)) { 098 transformToTridiagonal(matrix); 099 findEigenVectors(transformer.getQ().getData()); 100 } else { 101 // as of 2.0, non-symmetric matrices (i.e. complex eigenvalues) are 102 // NOT supported 103 // see issue https://issues.apache.org/jira/browse/MATH-235 104 throw new InvalidMatrixException( 105 "eigen decomposition of assymetric matrices not supported yet"); 106 } 107 } 108 109 /** 110 * Calculates the eigen decomposition of the symmetric tridiagonal 111 * matrix. The Householder matrix is assumed to be the identity matrix. 112 * @param main Main diagonal of the symmetric triadiagonal form 113 * @param secondary Secondary of the tridiagonal form 114 * @param splitTolerance dummy parameter, present for backward compatibility only. 115 * @exception InvalidMatrixException (wrapping a 116 * {@link org.apache.commons.math.ConvergenceException} if algorithm 117 * fails to converge 118 */ 119 public EigenDecompositionImpl(final double[] main,final double[] secondary, 120 final double splitTolerance) 121 throws InvalidMatrixException { 122 this.main = main.clone(); 123 this.secondary = secondary.clone(); 124 transformer = null; 125 final int size=main.length; 126 double[][] z = new double[size][size]; 127 for (int i=0;i<size;i++) { 128 z[i][i]=1.0; 129 } 130 findEigenVectors(z); 131 } 132 133 /** 134 * Check if a matrix is symmetric. 135 * @param matrix 136 * matrix to check 137 * @return true if matrix is symmetric 138 */ 139 private boolean isSymmetric(final RealMatrix matrix) { 140 final int rows = matrix.getRowDimension(); 141 final int columns = matrix.getColumnDimension(); 142 final double eps = 10 * rows * columns * MathUtils.EPSILON; 143 for (int i = 0; i < rows; ++i) { 144 for (int j = i + 1; j < columns; ++j) { 145 final double mij = matrix.getEntry(i, j); 146 final double mji = matrix.getEntry(j, i); 147 if (Math.abs(mij - mji) > (Math.max(Math.abs(mij), Math 148 .abs(mji)) * eps)) { 149 return false; 150 } 151 } 152 } 153 return true; 154 } 155 156 /** {@inheritDoc} */ 157 public RealMatrix getV() throws InvalidMatrixException { 158 159 if (cachedV == null) { 160 final int m = eigenvectors.length; 161 cachedV = MatrixUtils.createRealMatrix(m, m); 162 for (int k = 0; k < m; ++k) { 163 cachedV.setColumnVector(k, eigenvectors[k]); 164 } 165 } 166 // return the cached matrix 167 return cachedV; 168 169 } 170 171 /** {@inheritDoc} */ 172 public RealMatrix getD() throws InvalidMatrixException { 173 if (cachedD == null) { 174 // cache the matrix for subsequent calls 175 cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues); 176 } 177 return cachedD; 178 } 179 180 /** {@inheritDoc} */ 181 public RealMatrix getVT() throws InvalidMatrixException { 182 183 if (cachedVt == null) { 184 final int m = eigenvectors.length; 185 cachedVt = MatrixUtils.createRealMatrix(m, m); 186 for (int k = 0; k < m; ++k) { 187 cachedVt.setRowVector(k, eigenvectors[k]); 188 } 189 190 } 191 192 // return the cached matrix 193 return cachedVt; 194 } 195 196 /** {@inheritDoc} */ 197 public double[] getRealEigenvalues() throws InvalidMatrixException { 198 return realEigenvalues.clone(); 199 } 200 201 /** {@inheritDoc} */ 202 public double getRealEigenvalue(final int i) throws InvalidMatrixException, 203 ArrayIndexOutOfBoundsException { 204 return realEigenvalues[i]; 205 } 206 207 /** {@inheritDoc} */ 208 public double[] getImagEigenvalues() throws InvalidMatrixException { 209 return imagEigenvalues.clone(); 210 } 211 212 /** {@inheritDoc} */ 213 public double getImagEigenvalue(final int i) throws InvalidMatrixException, 214 ArrayIndexOutOfBoundsException { 215 return imagEigenvalues[i]; 216 } 217 218 /** {@inheritDoc} */ 219 public RealVector getEigenvector(final int i) 220 throws InvalidMatrixException, ArrayIndexOutOfBoundsException { 221 return eigenvectors[i].copy(); 222 } 223 224 /** 225 * Return the determinant of the matrix 226 * @return determinant of the matrix 227 */ 228 public double getDeterminant() { 229 double determinant = 1; 230 for (double lambda : realEigenvalues) { 231 determinant *= lambda; 232 } 233 return determinant; 234 } 235 236 /** {@inheritDoc} */ 237 public DecompositionSolver getSolver() { 238 return new Solver(realEigenvalues, imagEigenvalues, eigenvectors); 239 } 240 241 /** Specialized solver. */ 242 private static class Solver implements DecompositionSolver { 243 244 /** Real part of the realEigenvalues. */ 245 private double[] realEigenvalues; 246 247 /** Imaginary part of the realEigenvalues. */ 248 private double[] imagEigenvalues; 249 250 /** Eigenvectors. */ 251 private final ArrayRealVector[] eigenvectors; 252 253 /** 254 * Build a solver from decomposed matrix. 255 * @param realEigenvalues 256 * real parts of the eigenvalues 257 * @param imagEigenvalues 258 * imaginary parts of the eigenvalues 259 * @param eigenvectors 260 * eigenvectors 261 */ 262 private Solver(final double[] realEigenvalues, 263 final double[] imagEigenvalues, 264 final ArrayRealVector[] eigenvectors) { 265 this.realEigenvalues = realEigenvalues; 266 this.imagEigenvalues = imagEigenvalues; 267 this.eigenvectors = eigenvectors; 268 } 269 270 /** 271 * Solve the linear equation A × X = B for symmetric matrices A. 272 * <p> 273 * This method only find exact linear solutions, i.e. solutions for 274 * which ||A × X - B|| is exactly 0. 275 * </p> 276 * @param b 277 * right-hand side of the equation A × X = B 278 * @return a vector X that minimizes the two norm of A × X - B 279 * @exception IllegalArgumentException 280 * if matrices dimensions don't match 281 * @exception InvalidMatrixException 282 * if decomposed matrix is singular 283 */ 284 public double[] solve(final double[] b) 285 throws IllegalArgumentException, InvalidMatrixException { 286 287 if (!isNonSingular()) { 288 throw new SingularMatrixException(); 289 } 290 291 final int m = realEigenvalues.length; 292 if (b.length != m) { 293 throw MathRuntimeException.createIllegalArgumentException( 294 "vector length mismatch: got {0} but expected {1}", 295 b.length, m); 296 } 297 298 final double[] bp = new double[m]; 299 for (int i = 0; i < m; ++i) { 300 final ArrayRealVector v = eigenvectors[i]; 301 final double[] vData = v.getDataRef(); 302 final double s = v.dotProduct(b) / realEigenvalues[i]; 303 for (int j = 0; j < m; ++j) { 304 bp[j] += s * vData[j]; 305 } 306 } 307 308 return bp; 309 310 } 311 312 /** 313 * Solve the linear equation A × X = B for symmetric matrices A. 314 * <p> 315 * This method only find exact linear solutions, i.e. solutions for 316 * which ||A × X - B|| is exactly 0. 317 * </p> 318 * @param b 319 * right-hand side of the equation A × X = B 320 * @return a vector X that minimizes the two norm of A × X - B 321 * @exception IllegalArgumentException 322 * if matrices dimensions don't match 323 * @exception InvalidMatrixException 324 * if decomposed matrix is singular 325 */ 326 public RealVector solve(final RealVector b) 327 throws IllegalArgumentException, InvalidMatrixException { 328 329 if (!isNonSingular()) { 330 throw new SingularMatrixException(); 331 } 332 333 final int m = realEigenvalues.length; 334 if (b.getDimension() != m) { 335 throw MathRuntimeException.createIllegalArgumentException( 336 "vector length mismatch: got {0} but expected {1}", b 337 .getDimension(), m); 338 } 339 340 final double[] bp = new double[m]; 341 for (int i = 0; i < m; ++i) { 342 final ArrayRealVector v = eigenvectors[i]; 343 final double[] vData = v.getDataRef(); 344 final double s = v.dotProduct(b) / realEigenvalues[i]; 345 for (int j = 0; j < m; ++j) { 346 bp[j] += s * vData[j]; 347 } 348 } 349 350 return new ArrayRealVector(bp, false); 351 352 } 353 354 /** 355 * Solve the linear equation A × X = B for symmetric matrices A. 356 * <p> 357 * This method only find exact linear solutions, i.e. solutions for 358 * which ||A × X - B|| is exactly 0. 359 * </p> 360 * @param b 361 * right-hand side of the equation A × X = B 362 * @return a matrix X that minimizes the two norm of A × X - B 363 * @exception IllegalArgumentException 364 * if matrices dimensions don't match 365 * @exception InvalidMatrixException 366 * if decomposed matrix is singular 367 */ 368 public RealMatrix solve(final RealMatrix b) 369 throws IllegalArgumentException, InvalidMatrixException { 370 371 if (!isNonSingular()) { 372 throw new SingularMatrixException(); 373 } 374 375 final int m = realEigenvalues.length; 376 if (b.getRowDimension() != m) { 377 throw MathRuntimeException 378 .createIllegalArgumentException( 379 "dimensions mismatch: got {0}x{1} but expected {2}x{3}", 380 b.getRowDimension(), b.getColumnDimension(), m, 381 "n"); 382 } 383 384 final int nColB = b.getColumnDimension(); 385 final double[][] bp = new double[m][nColB]; 386 for (int k = 0; k < nColB; ++k) { 387 for (int i = 0; i < m; ++i) { 388 final ArrayRealVector v = eigenvectors[i]; 389 final double[] vData = v.getDataRef(); 390 double s = 0; 391 for (int j = 0; j < m; ++j) { 392 s += v.getEntry(j) * b.getEntry(j, k); 393 } 394 s /= realEigenvalues[i]; 395 for (int j = 0; j < m; ++j) { 396 bp[j][k] += s * vData[j]; 397 } 398 } 399 } 400 401 return MatrixUtils.createRealMatrix(bp); 402 403 } 404 405 /** 406 * Check if the decomposed matrix is non-singular. 407 * @return true if the decomposed matrix is non-singular 408 */ 409 public boolean isNonSingular() { 410 for (int i = 0; i < realEigenvalues.length; ++i) { 411 if ((realEigenvalues[i] == 0) && (imagEigenvalues[i] == 0)) { 412 return false; 413 } 414 } 415 return true; 416 } 417 418 /** 419 * Get the inverse of the decomposed matrix. 420 * @return inverse matrix 421 * @throws InvalidMatrixException 422 * if decomposed matrix is singular 423 */ 424 public RealMatrix getInverse() throws InvalidMatrixException { 425 426 if (!isNonSingular()) { 427 throw new SingularMatrixException(); 428 } 429 430 final int m = realEigenvalues.length; 431 final double[][] invData = new double[m][m]; 432 433 for (int i = 0; i < m; ++i) { 434 final double[] invI = invData[i]; 435 for (int j = 0; j < m; ++j) { 436 double invIJ = 0; 437 for (int k = 0; k < m; ++k) { 438 final double[] vK = eigenvectors[k].getDataRef(); 439 invIJ += vK[i] * vK[j] / realEigenvalues[k]; 440 } 441 invI[j] = invIJ; 442 } 443 } 444 return MatrixUtils.createRealMatrix(invData); 445 446 } 447 448 } 449 450 /** 451 * Transform matrix to tridiagonal. 452 * @param matrix 453 * matrix to transform 454 */ 455 private void transformToTridiagonal(final RealMatrix matrix) { 456 457 // transform the matrix to tridiagonal 458 transformer = new TriDiagonalTransformer(matrix); 459 main = transformer.getMainDiagonalRef(); 460 secondary = transformer.getSecondaryDiagonalRef(); 461 462 } 463 464 /** 465 * Find eigenvalues and eigenvectors (Dubrulle et al., 1971) 466 * @param householderMatrix Householder matrix of the transformation 467 * to tri-diagonal form. 468 */ 469 private void findEigenVectors(double[][] householderMatrix) { 470 471 double[][]z = householderMatrix.clone(); 472 final int n = main.length; 473 realEigenvalues = new double[n]; 474 imagEigenvalues = new double[n]; 475 double[] e = new double[n]; 476 for (int i = 0; i < n - 1; i++) { 477 realEigenvalues[i] = main[i]; 478 e[i] = secondary[i]; 479 } 480 realEigenvalues[n - 1] = main[n - 1]; 481 e[n - 1] = 0.0; 482 483 // Determine the largest main and secondary value in absolute term. 484 double maxAbsoluteValue=0.0; 485 for (int i = 0; i < n; i++) { 486 if (Math.abs(realEigenvalues[i])>maxAbsoluteValue) { 487 maxAbsoluteValue=Math.abs(realEigenvalues[i]); 488 } 489 if (Math.abs(e[i])>maxAbsoluteValue) { 490 maxAbsoluteValue=Math.abs(e[i]); 491 } 492 } 493 // Make null any main and secondary value too small to be significant 494 if (maxAbsoluteValue!=0.0) { 495 for (int i=0; i < n; i++) { 496 if (Math.abs(realEigenvalues[i])<=MathUtils.EPSILON*maxAbsoluteValue) { 497 realEigenvalues[i]=0.0; 498 } 499 if (Math.abs(e[i])<=MathUtils.EPSILON*maxAbsoluteValue) { 500 e[i]=0.0; 501 } 502 } 503 } 504 505 for (int j = 0; j < n; j++) { 506 int its = 0; 507 int m; 508 do { 509 for (m = j; m < n - 1; m++) { 510 double delta = Math.abs(realEigenvalues[m]) + Math.abs(realEigenvalues[m + 1]); 511 if (Math.abs(e[m]) + delta == delta) { 512 break; 513 } 514 } 515 if (m != j) { 516 if (its == maxIter) 517 throw new InvalidMatrixException( 518 new MaxIterationsExceededException(maxIter)); 519 its++; 520 double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]); 521 double t = Math.sqrt(1 + q * q); 522 if (q < 0.0) { 523 q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q - t); 524 } else { 525 q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q + t); 526 } 527 double u = 0.0; 528 double s = 1.0; 529 double c = 1.0; 530 int i; 531 for (i = m - 1; i >= j; i--) { 532 double p = s * e[i]; 533 double h = c * e[i]; 534 if (Math.abs(p) >= Math.abs(q)) { 535 c = q / p; 536 t = Math.sqrt(c * c + 1.0); 537 e[i + 1] = p * t; 538 s = 1.0 / t; 539 c = c * s; 540 } else { 541 s = p / q; 542 t = Math.sqrt(s * s + 1.0); 543 e[i + 1] = q * t; 544 c = 1.0 / t; 545 s = s * c; 546 } 547 if (e[i + 1] == 0.0) { 548 realEigenvalues[i + 1] -= u; 549 e[m] = 0.0; 550 break; 551 } 552 q = realEigenvalues[i + 1] - u; 553 t = (realEigenvalues[i] - q) * s + 2.0 * c * h; 554 u = s * t; 555 realEigenvalues[i + 1] = q + u; 556 q = c * t - h; 557 for (int ia = 0; ia < n; ia++) { 558 p = z[ia][i + 1]; 559 z[ia][i + 1] = s * z[ia][i] + c * p; 560 z[ia][i] = c * z[ia][i] - s * p; 561 } 562 } 563 if (e[i + 1] == 0.0 && i >= j) 564 continue; 565 realEigenvalues[j] -= u; 566 e[j] = q; 567 e[m] = 0.0; 568 } 569 } while (m != j); 570 } 571 572 //Sort the eigen values (and vectors) in increase order 573 for (int i = 0; i < n; i++) { 574 int k = i; 575 double p = realEigenvalues[i]; 576 for (int j = i + 1; j < n; j++) { 577 if (realEigenvalues[j] > p) { 578 k = j; 579 p = realEigenvalues[j]; 580 } 581 } 582 if (k != i) { 583 realEigenvalues[k] = realEigenvalues[i]; 584 realEigenvalues[i] = p; 585 for (int j = 0; j < n; j++) { 586 p = z[j][i]; 587 z[j][i] = z[j][k]; 588 z[j][k] = p; 589 } 590 } 591 } 592 593 // Determine the largest eigen value in absolute term. 594 maxAbsoluteValue=0.0; 595 for (int i = 0; i < n; i++) { 596 if (Math.abs(realEigenvalues[i])>maxAbsoluteValue) { 597 maxAbsoluteValue=Math.abs(realEigenvalues[i]); 598 } 599 } 600 // Make null any eigen value too small to be significant 601 if (maxAbsoluteValue!=0.0) { 602 for (int i=0; i < n; i++) { 603 if (Math.abs(realEigenvalues[i])<MathUtils.EPSILON*maxAbsoluteValue) { 604 realEigenvalues[i]=0.0; 605 } 606 } 607 } 608 eigenvectors = new ArrayRealVector[n]; 609 double[] tmp = new double[n]; 610 for (int i = 0; i < n; i++) { 611 for (int j = 0; j < n; j++) { 612 tmp[j] = z[j][i]; 613 } 614 eigenvectors[i] = new ArrayRealVector(tmp); 615 } 616 } 617 }