1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 18 package org.apache.commons.math.linear; 19 20 import org.apache.commons.math.MathRuntimeException; 21 import org.apache.commons.math.MaxIterationsExceededException; 22 import org.apache.commons.math.exception.util.LocalizedFormats; 23 import org.apache.commons.math.util.MathUtils; 24 import org.apache.commons.math.util.FastMath; 25 26 /** 27 * Calculates the eigen decomposition of a real <strong>symmetric</strong> 28 * matrix. 29 * <p> 30 * The eigen decomposition of matrix A is a set of two matrices: V and D such 31 * that A = V D V<sup>T</sup>. A, V and D are all m × m matrices. 32 * </p> 33 * <p> 34 * As of 2.0, this class supports only <strong>symmetric</strong> matrices, and 35 * hence computes only real realEigenvalues. This implies the D matrix returned 36 * by {@link #getD()} is always diagonal and the imaginary values returned 37 * {@link #getImagEigenvalue(int)} and {@link #getImagEigenvalues()} are always 38 * null. 39 * </p> 40 * <p> 41 * When called with a {@link RealMatrix} argument, this implementation only uses 42 * the upper part of the matrix, the part below the diagonal is not accessed at 43 * all. 44 * </p> 45 * <p> 46 * This implementation is based on the paper by A. Drubrulle, R.S. Martin and 47 * J.H. Wilkinson 'The Implicit QL Algorithm' in Wilksinson and Reinsch (1971) 48 * Handbook for automatic computation, vol. 2, Linear algebra, Springer-Verlag, 49 * New-York 50 * </p> 51 * @version $Revision: 1002040 $ $Date: 2010-09-28 09:18:31 +0200 (mar. 28 sept. 2010) $ 52 * @since 2.0 53 */ 54 public class EigenDecompositionImpl implements EigenDecomposition { 55 56 /** Maximum number of iterations accepted in the implicit QL transformation */ 57 private byte maxIter = 30; 58 59 /** Main diagonal of the tridiagonal matrix. */ 60 private double[] main; 61 62 /** Secondary diagonal of the tridiagonal matrix. */ 63 private double[] secondary; 64 65 /** 66 * Transformer to tridiagonal (may be null if matrix is already 67 * tridiagonal). 68 */ 69 private TriDiagonalTransformer transformer; 70 71 /** Real part of the realEigenvalues. */ 72 private double[] realEigenvalues; 73 74 /** Imaginary part of the realEigenvalues. */ 75 private double[] imagEigenvalues; 76 77 /** Eigenvectors. */ 78 private ArrayRealVector[] eigenvectors; 79 80 /** Cached value of V. */ 81 private RealMatrix cachedV; 82 83 /** Cached value of D. */ 84 private RealMatrix cachedD; 85 86 /** Cached value of Vt. */ 87 private RealMatrix cachedVt; 88 89 /** 90 * Calculates the eigen decomposition of the given symmetric matrix. 91 * @param matrix The <strong>symmetric</strong> matrix to decompose. 92 * @param splitTolerance dummy parameter, present for backward compatibility only. 93 * @exception InvalidMatrixException (wrapping a 94 * {@link org.apache.commons.math.ConvergenceException} if algorithm 95 * fails to converge 96 */ EigenDecompositionImpl(final RealMatrix matrix,final double splitTolerance)97 public EigenDecompositionImpl(final RealMatrix matrix,final double splitTolerance) 98 throws InvalidMatrixException { 99 if (isSymmetric(matrix)) { 100 transformToTridiagonal(matrix); 101 findEigenVectors(transformer.getQ().getData()); 102 } else { 103 // as of 2.0, non-symmetric matrices (i.e. complex eigenvalues) are 104 // NOT supported 105 // see issue https://issues.apache.org/jira/browse/MATH-235 106 throw new InvalidMatrixException( 107 LocalizedFormats.ASSYMETRIC_EIGEN_NOT_SUPPORTED); 108 } 109 } 110 111 /** 112 * Calculates the eigen decomposition of the symmetric tridiagonal 113 * matrix. The Householder matrix is assumed to be the identity matrix. 114 * @param main Main diagonal of the symmetric triadiagonal form 115 * @param secondary Secondary of the tridiagonal form 116 * @param splitTolerance dummy parameter, present for backward compatibility only. 117 * @exception InvalidMatrixException (wrapping a 118 * {@link org.apache.commons.math.ConvergenceException} if algorithm 119 * fails to converge 120 */ EigenDecompositionImpl(final double[] main,final double[] secondary, final double splitTolerance)121 public EigenDecompositionImpl(final double[] main,final double[] secondary, 122 final double splitTolerance) 123 throws InvalidMatrixException { 124 this.main = main.clone(); 125 this.secondary = secondary.clone(); 126 transformer = null; 127 final int size=main.length; 128 double[][] z = new double[size][size]; 129 for (int i=0;i<size;i++) { 130 z[i][i]=1.0; 131 } 132 findEigenVectors(z); 133 } 134 135 /** 136 * Check if a matrix is symmetric. 137 * @param matrix 138 * matrix to check 139 * @return true if matrix is symmetric 140 */ isSymmetric(final RealMatrix matrix)141 private boolean isSymmetric(final RealMatrix matrix) { 142 final int rows = matrix.getRowDimension(); 143 final int columns = matrix.getColumnDimension(); 144 final double eps = 10 * rows * columns * MathUtils.EPSILON; 145 for (int i = 0; i < rows; ++i) { 146 for (int j = i + 1; j < columns; ++j) { 147 final double mij = matrix.getEntry(i, j); 148 final double mji = matrix.getEntry(j, i); 149 if (FastMath.abs(mij - mji) > 150 (FastMath.max(FastMath.abs(mij), FastMath.abs(mji)) * eps)) { 151 return false; 152 } 153 } 154 } 155 return true; 156 } 157 158 /** {@inheritDoc} */ getV()159 public RealMatrix getV() throws InvalidMatrixException { 160 161 if (cachedV == null) { 162 final int m = eigenvectors.length; 163 cachedV = MatrixUtils.createRealMatrix(m, m); 164 for (int k = 0; k < m; ++k) { 165 cachedV.setColumnVector(k, eigenvectors[k]); 166 } 167 } 168 // return the cached matrix 169 return cachedV; 170 171 } 172 173 /** {@inheritDoc} */ getD()174 public RealMatrix getD() throws InvalidMatrixException { 175 if (cachedD == null) { 176 // cache the matrix for subsequent calls 177 cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues); 178 } 179 return cachedD; 180 } 181 182 /** {@inheritDoc} */ getVT()183 public RealMatrix getVT() throws InvalidMatrixException { 184 185 if (cachedVt == null) { 186 final int m = eigenvectors.length; 187 cachedVt = MatrixUtils.createRealMatrix(m, m); 188 for (int k = 0; k < m; ++k) { 189 cachedVt.setRowVector(k, eigenvectors[k]); 190 } 191 192 } 193 194 // return the cached matrix 195 return cachedVt; 196 } 197 198 /** {@inheritDoc} */ getRealEigenvalues()199 public double[] getRealEigenvalues() throws InvalidMatrixException { 200 return realEigenvalues.clone(); 201 } 202 203 /** {@inheritDoc} */ getRealEigenvalue(final int i)204 public double getRealEigenvalue(final int i) throws InvalidMatrixException, 205 ArrayIndexOutOfBoundsException { 206 return realEigenvalues[i]; 207 } 208 209 /** {@inheritDoc} */ getImagEigenvalues()210 public double[] getImagEigenvalues() throws InvalidMatrixException { 211 return imagEigenvalues.clone(); 212 } 213 214 /** {@inheritDoc} */ getImagEigenvalue(final int i)215 public double getImagEigenvalue(final int i) throws InvalidMatrixException, 216 ArrayIndexOutOfBoundsException { 217 return imagEigenvalues[i]; 218 } 219 220 /** {@inheritDoc} */ getEigenvector(final int i)221 public RealVector getEigenvector(final int i) 222 throws InvalidMatrixException, ArrayIndexOutOfBoundsException { 223 return eigenvectors[i].copy(); 224 } 225 226 /** 227 * Return the determinant of the matrix 228 * @return determinant of the matrix 229 */ getDeterminant()230 public double getDeterminant() { 231 double determinant = 1; 232 for (double lambda : realEigenvalues) { 233 determinant *= lambda; 234 } 235 return determinant; 236 } 237 238 /** {@inheritDoc} */ getSolver()239 public DecompositionSolver getSolver() { 240 return new Solver(realEigenvalues, imagEigenvalues, eigenvectors); 241 } 242 243 /** Specialized solver. */ 244 private static class Solver implements DecompositionSolver { 245 246 /** Real part of the realEigenvalues. */ 247 private double[] realEigenvalues; 248 249 /** Imaginary part of the realEigenvalues. */ 250 private double[] imagEigenvalues; 251 252 /** Eigenvectors. */ 253 private final ArrayRealVector[] eigenvectors; 254 255 /** 256 * Build a solver from decomposed matrix. 257 * @param realEigenvalues 258 * real parts of the eigenvalues 259 * @param imagEigenvalues 260 * imaginary parts of the eigenvalues 261 * @param eigenvectors 262 * eigenvectors 263 */ Solver(final double[] realEigenvalues, final double[] imagEigenvalues, final ArrayRealVector[] eigenvectors)264 private Solver(final double[] realEigenvalues, 265 final double[] imagEigenvalues, 266 final ArrayRealVector[] eigenvectors) { 267 this.realEigenvalues = realEigenvalues; 268 this.imagEigenvalues = imagEigenvalues; 269 this.eigenvectors = eigenvectors; 270 } 271 272 /** 273 * Solve the linear equation A × X = B for symmetric matrices A. 274 * <p> 275 * This method only find exact linear solutions, i.e. solutions for 276 * which ||A × X - B|| is exactly 0. 277 * </p> 278 * @param b 279 * right-hand side of the equation A × X = B 280 * @return a vector X that minimizes the two norm of A × X - B 281 * @exception IllegalArgumentException 282 * if matrices dimensions don't match 283 * @exception InvalidMatrixException 284 * if decomposed matrix is singular 285 */ solve(final double[] b)286 public double[] solve(final double[] b) 287 throws IllegalArgumentException, InvalidMatrixException { 288 289 if (!isNonSingular()) { 290 throw new SingularMatrixException(); 291 } 292 293 final int m = realEigenvalues.length; 294 if (b.length != m) { 295 throw MathRuntimeException.createIllegalArgumentException( 296 LocalizedFormats.VECTOR_LENGTH_MISMATCH, 297 b.length, m); 298 } 299 300 final double[] bp = new double[m]; 301 for (int i = 0; i < m; ++i) { 302 final ArrayRealVector v = eigenvectors[i]; 303 final double[] vData = v.getDataRef(); 304 final double s = v.dotProduct(b) / realEigenvalues[i]; 305 for (int j = 0; j < m; ++j) { 306 bp[j] += s * vData[j]; 307 } 308 } 309 310 return bp; 311 312 } 313 314 /** 315 * Solve the linear equation A × X = B for symmetric matrices A. 316 * <p> 317 * This method only find exact linear solutions, i.e. solutions for 318 * which ||A × X - B|| is exactly 0. 319 * </p> 320 * @param b 321 * right-hand side of the equation A × X = B 322 * @return a vector X that minimizes the two norm of A × X - B 323 * @exception IllegalArgumentException 324 * if matrices dimensions don't match 325 * @exception InvalidMatrixException 326 * if decomposed matrix is singular 327 */ solve(final RealVector b)328 public RealVector solve(final RealVector b) 329 throws IllegalArgumentException, InvalidMatrixException { 330 331 if (!isNonSingular()) { 332 throw new SingularMatrixException(); 333 } 334 335 final int m = realEigenvalues.length; 336 if (b.getDimension() != m) { 337 throw MathRuntimeException.createIllegalArgumentException( 338 LocalizedFormats.VECTOR_LENGTH_MISMATCH, b 339 .getDimension(), m); 340 } 341 342 final double[] bp = new double[m]; 343 for (int i = 0; i < m; ++i) { 344 final ArrayRealVector v = eigenvectors[i]; 345 final double[] vData = v.getDataRef(); 346 final double s = v.dotProduct(b) / realEigenvalues[i]; 347 for (int j = 0; j < m; ++j) { 348 bp[j] += s * vData[j]; 349 } 350 } 351 352 return new ArrayRealVector(bp, false); 353 354 } 355 356 /** 357 * Solve the linear equation A × X = B for symmetric matrices A. 358 * <p> 359 * This method only find exact linear solutions, i.e. solutions for 360 * which ||A × X - B|| is exactly 0. 361 * </p> 362 * @param b 363 * right-hand side of the equation A × X = B 364 * @return a matrix X that minimizes the two norm of A × X - B 365 * @exception IllegalArgumentException 366 * if matrices dimensions don't match 367 * @exception InvalidMatrixException 368 * if decomposed matrix is singular 369 */ solve(final RealMatrix b)370 public RealMatrix solve(final RealMatrix b) 371 throws IllegalArgumentException, InvalidMatrixException { 372 373 if (!isNonSingular()) { 374 throw new SingularMatrixException(); 375 } 376 377 final int m = realEigenvalues.length; 378 if (b.getRowDimension() != m) { 379 throw MathRuntimeException 380 .createIllegalArgumentException( 381 LocalizedFormats.DIMENSIONS_MISMATCH_2x2, 382 b.getRowDimension(), b.getColumnDimension(), m, 383 "n"); 384 } 385 386 final int nColB = b.getColumnDimension(); 387 final double[][] bp = new double[m][nColB]; 388 for (int k = 0; k < nColB; ++k) { 389 for (int i = 0; i < m; ++i) { 390 final ArrayRealVector v = eigenvectors[i]; 391 final double[] vData = v.getDataRef(); 392 double s = 0; 393 for (int j = 0; j < m; ++j) { 394 s += v.getEntry(j) * b.getEntry(j, k); 395 } 396 s /= realEigenvalues[i]; 397 for (int j = 0; j < m; ++j) { 398 bp[j][k] += s * vData[j]; 399 } 400 } 401 } 402 403 return MatrixUtils.createRealMatrix(bp); 404 405 } 406 407 /** 408 * Check if the decomposed matrix is non-singular. 409 * @return true if the decomposed matrix is non-singular 410 */ isNonSingular()411 public boolean isNonSingular() { 412 for (int i = 0; i < realEigenvalues.length; ++i) { 413 if ((realEigenvalues[i] == 0) && (imagEigenvalues[i] == 0)) { 414 return false; 415 } 416 } 417 return true; 418 } 419 420 /** 421 * Get the inverse of the decomposed matrix. 422 * @return inverse matrix 423 * @throws InvalidMatrixException 424 * if decomposed matrix is singular 425 */ getInverse()426 public RealMatrix getInverse() throws InvalidMatrixException { 427 428 if (!isNonSingular()) { 429 throw new SingularMatrixException(); 430 } 431 432 final int m = realEigenvalues.length; 433 final double[][] invData = new double[m][m]; 434 435 for (int i = 0; i < m; ++i) { 436 final double[] invI = invData[i]; 437 for (int j = 0; j < m; ++j) { 438 double invIJ = 0; 439 for (int k = 0; k < m; ++k) { 440 final double[] vK = eigenvectors[k].getDataRef(); 441 invIJ += vK[i] * vK[j] / realEigenvalues[k]; 442 } 443 invI[j] = invIJ; 444 } 445 } 446 return MatrixUtils.createRealMatrix(invData); 447 448 } 449 450 } 451 452 /** 453 * Transform matrix to tridiagonal. 454 * @param matrix 455 * matrix to transform 456 */ transformToTridiagonal(final RealMatrix matrix)457 private void transformToTridiagonal(final RealMatrix matrix) { 458 459 // transform the matrix to tridiagonal 460 transformer = new TriDiagonalTransformer(matrix); 461 main = transformer.getMainDiagonalRef(); 462 secondary = transformer.getSecondaryDiagonalRef(); 463 464 } 465 466 /** 467 * Find eigenvalues and eigenvectors (Dubrulle et al., 1971) 468 * @param householderMatrix Householder matrix of the transformation 469 * to tri-diagonal form. 470 */ findEigenVectors(double[][] householderMatrix)471 private void findEigenVectors(double[][] householderMatrix) { 472 473 double[][]z = householderMatrix.clone(); 474 final int n = main.length; 475 realEigenvalues = new double[n]; 476 imagEigenvalues = new double[n]; 477 double[] e = new double[n]; 478 for (int i = 0; i < n - 1; i++) { 479 realEigenvalues[i] = main[i]; 480 e[i] = secondary[i]; 481 } 482 realEigenvalues[n - 1] = main[n - 1]; 483 e[n - 1] = 0.0; 484 485 // Determine the largest main and secondary value in absolute term. 486 double maxAbsoluteValue=0.0; 487 for (int i = 0; i < n; i++) { 488 if (FastMath.abs(realEigenvalues[i])>maxAbsoluteValue) { 489 maxAbsoluteValue=FastMath.abs(realEigenvalues[i]); 490 } 491 if (FastMath.abs(e[i])>maxAbsoluteValue) { 492 maxAbsoluteValue=FastMath.abs(e[i]); 493 } 494 } 495 // Make null any main and secondary value too small to be significant 496 if (maxAbsoluteValue!=0.0) { 497 for (int i=0; i < n; i++) { 498 if (FastMath.abs(realEigenvalues[i])<=MathUtils.EPSILON*maxAbsoluteValue) { 499 realEigenvalues[i]=0.0; 500 } 501 if (FastMath.abs(e[i])<=MathUtils.EPSILON*maxAbsoluteValue) { 502 e[i]=0.0; 503 } 504 } 505 } 506 507 for (int j = 0; j < n; j++) { 508 int its = 0; 509 int m; 510 do { 511 for (m = j; m < n - 1; m++) { 512 double delta = FastMath.abs(realEigenvalues[m]) + FastMath.abs(realEigenvalues[m + 1]); 513 if (FastMath.abs(e[m]) + delta == delta) { 514 break; 515 } 516 } 517 if (m != j) { 518 if (its == maxIter) 519 throw new InvalidMatrixException( 520 new MaxIterationsExceededException(maxIter)); 521 its++; 522 double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]); 523 double t = FastMath.sqrt(1 + q * q); 524 if (q < 0.0) { 525 q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q - t); 526 } else { 527 q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q + t); 528 } 529 double u = 0.0; 530 double s = 1.0; 531 double c = 1.0; 532 int i; 533 for (i = m - 1; i >= j; i--) { 534 double p = s * e[i]; 535 double h = c * e[i]; 536 if (FastMath.abs(p) >= FastMath.abs(q)) { 537 c = q / p; 538 t = FastMath.sqrt(c * c + 1.0); 539 e[i + 1] = p * t; 540 s = 1.0 / t; 541 c = c * s; 542 } else { 543 s = p / q; 544 t = FastMath.sqrt(s * s + 1.0); 545 e[i + 1] = q * t; 546 c = 1.0 / t; 547 s = s * c; 548 } 549 if (e[i + 1] == 0.0) { 550 realEigenvalues[i + 1] -= u; 551 e[m] = 0.0; 552 break; 553 } 554 q = realEigenvalues[i + 1] - u; 555 t = (realEigenvalues[i] - q) * s + 2.0 * c * h; 556 u = s * t; 557 realEigenvalues[i + 1] = q + u; 558 q = c * t - h; 559 for (int ia = 0; ia < n; ia++) { 560 p = z[ia][i + 1]; 561 z[ia][i + 1] = s * z[ia][i] + c * p; 562 z[ia][i] = c * z[ia][i] - s * p; 563 } 564 } 565 if (t == 0.0 && i >= j) 566 continue; 567 realEigenvalues[j] -= u; 568 e[j] = q; 569 e[m] = 0.0; 570 } 571 } while (m != j); 572 } 573 574 //Sort the eigen values (and vectors) in increase order 575 for (int i = 0; i < n; i++) { 576 int k = i; 577 double p = realEigenvalues[i]; 578 for (int j = i + 1; j < n; j++) { 579 if (realEigenvalues[j] > p) { 580 k = j; 581 p = realEigenvalues[j]; 582 } 583 } 584 if (k != i) { 585 realEigenvalues[k] = realEigenvalues[i]; 586 realEigenvalues[i] = p; 587 for (int j = 0; j < n; j++) { 588 p = z[j][i]; 589 z[j][i] = z[j][k]; 590 z[j][k] = p; 591 } 592 } 593 } 594 595 // Determine the largest eigen value in absolute term. 596 maxAbsoluteValue=0.0; 597 for (int i = 0; i < n; i++) { 598 if (FastMath.abs(realEigenvalues[i])>maxAbsoluteValue) { 599 maxAbsoluteValue=FastMath.abs(realEigenvalues[i]); 600 } 601 } 602 // Make null any eigen value too small to be significant 603 if (maxAbsoluteValue!=0.0) { 604 for (int i=0; i < n; i++) { 605 if (FastMath.abs(realEigenvalues[i])<MathUtils.EPSILON*maxAbsoluteValue) { 606 realEigenvalues[i]=0.0; 607 } 608 } 609 } 610 eigenvectors = new ArrayRealVector[n]; 611 double[] tmp = new double[n]; 612 for (int i = 0; i < n; i++) { 613 for (int j = 0; j < n; j++) { 614 tmp[j] = z[j][i]; 615 } 616 eigenvectors[i] = new ArrayRealVector(tmp); 617 } 618 } 619 } 620