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 018package org.apache.commons.math3.linear; 019 020import org.apache.commons.math3.complex.Complex; 021import org.apache.commons.math3.exception.MathArithmeticException; 022import org.apache.commons.math3.exception.MathUnsupportedOperationException; 023import org.apache.commons.math3.exception.MaxCountExceededException; 024import org.apache.commons.math3.exception.DimensionMismatchException; 025import org.apache.commons.math3.exception.util.LocalizedFormats; 026import org.apache.commons.math3.util.Precision; 027import org.apache.commons.math3.util.FastMath; 028 029/** 030 * Calculates the eigen decomposition of a real matrix. 031 * <p>The eigen decomposition of matrix A is a set of two matrices: 032 * V and D such that A = V × D × V<sup>T</sup>. 033 * A, V and D are all m × m matrices.</p> 034 * <p>This class is similar in spirit to the <code>EigenvalueDecomposition</code> 035 * class from the <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> 036 * library, with the following changes:</p> 037 * <ul> 038 * <li>a {@link #getVT() getVt} method has been added,</li> 039 * <li>two {@link #getRealEigenvalue(int) getRealEigenvalue} and {@link #getImagEigenvalue(int) 040 * getImagEigenvalue} methods to pick up a single eigenvalue have been added,</li> 041 * <li>a {@link #getEigenvector(int) getEigenvector} method to pick up a single 042 * eigenvector has been added,</li> 043 * <li>a {@link #getDeterminant() getDeterminant} method has been added.</li> 044 * <li>a {@link #getSolver() getSolver} method has been added.</li> 045 * </ul> 046 * <p> 047 * As of 3.1, this class supports general real matrices (both symmetric and non-symmetric): 048 * </p> 049 * <p> 050 * If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is diagonal and the eigenvector 051 * matrix V is orthogonal, i.e. A = V.multiply(D.multiply(V.transpose())) and 052 * V.multiply(V.transpose()) equals the identity matrix. 053 * </p> 054 * <p> 055 * If A is not symmetric, then the eigenvalue matrix D is block diagonal with the real eigenvalues 056 * in 1-by-1 blocks and any complex eigenvalues, lambda + i*mu, in 2-by-2 blocks: 057 * <pre> 058 * [lambda, mu ] 059 * [ -mu, lambda] 060 * </pre> 061 * The columns of V represent the eigenvectors in the sense that A*V = V*D, 062 * i.e. A.multiply(V) equals V.multiply(D). 063 * The matrix V may be badly conditioned, or even singular, so the validity of the equation 064 * A = V*D*inverse(V) depends upon the condition of V. 065 * </p> 066 * <p> 067 * This implementation is based on the paper by A. Drubrulle, R.S. Martin and 068 * J.H. Wilkinson "The Implicit QL Algorithm" in Wilksinson and Reinsch (1971) 069 * Handbook for automatic computation, vol. 2, Linear algebra, Springer-Verlag, 070 * New-York 071 * </p> 072 * @see <a href="http://mathworld.wolfram.com/EigenDecomposition.html">MathWorld</a> 073 * @see <a href="http://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix">Wikipedia</a> 074 * @since 2.0 (changed to concrete class in 3.0) 075 */ 076public class EigenDecomposition { 077 /** Internally used epsilon criteria. */ 078 private static final double EPSILON = 1e-12; 079 /** Maximum number of iterations accepted in the implicit QL transformation */ 080 private byte maxIter = 30; 081 /** Main diagonal of the tridiagonal matrix. */ 082 private double[] main; 083 /** Secondary diagonal of the tridiagonal matrix. */ 084 private double[] secondary; 085 /** 086 * Transformer to tridiagonal (may be null if matrix is already 087 * tridiagonal). 088 */ 089 private TriDiagonalTransformer transformer; 090 /** Real part of the realEigenvalues. */ 091 private double[] realEigenvalues; 092 /** Imaginary part of the realEigenvalues. */ 093 private double[] imagEigenvalues; 094 /** Eigenvectors. */ 095 private ArrayRealVector[] eigenvectors; 096 /** Cached value of V. */ 097 private RealMatrix cachedV; 098 /** Cached value of D. */ 099 private RealMatrix cachedD; 100 /** Cached value of Vt. */ 101 private RealMatrix cachedVt; 102 /** Whether the matrix is symmetric. */ 103 private final boolean isSymmetric; 104 105 /** 106 * Calculates the eigen decomposition of the given real matrix. 107 * <p> 108 * Supports decomposition of a general matrix since 3.1. 109 * 110 * @param matrix Matrix to decompose. 111 * @throws MaxCountExceededException if the algorithm fails to converge. 112 * @throws MathArithmeticException if the decomposition of a general matrix 113 * results in a matrix with zero norm 114 * @since 3.1 115 */ 116 public EigenDecomposition(final RealMatrix matrix) 117 throws MathArithmeticException { 118 final double symTol = 10 * matrix.getRowDimension() * matrix.getColumnDimension() * Precision.EPSILON; 119 isSymmetric = MatrixUtils.isSymmetric(matrix, symTol); 120 if (isSymmetric) { 121 transformToTridiagonal(matrix); 122 findEigenVectors(transformer.getQ().getData()); 123 } else { 124 final SchurTransformer t = transformToSchur(matrix); 125 findEigenVectorsFromSchur(t); 126 } 127 } 128 129 /** 130 * Calculates the eigen decomposition of the given real matrix. 131 * 132 * @param matrix Matrix to decompose. 133 * @param splitTolerance Dummy parameter (present for backward 134 * compatibility only). 135 * @throws MathArithmeticException if the decomposition of a general matrix 136 * results in a matrix with zero norm 137 * @throws MaxCountExceededException if the algorithm fails to converge. 138 * @deprecated in 3.1 (to be removed in 4.0) due to unused parameter 139 */ 140 @Deprecated 141 public EigenDecomposition(final RealMatrix matrix, 142 final double splitTolerance) 143 throws MathArithmeticException { 144 this(matrix); 145 } 146 147 /** 148 * Calculates the eigen decomposition of the symmetric tridiagonal 149 * matrix. The Householder matrix is assumed to be the identity matrix. 150 * 151 * @param main Main diagonal of the symmetric tridiagonal form. 152 * @param secondary Secondary of the tridiagonal form. 153 * @throws MaxCountExceededException if the algorithm fails to converge. 154 * @since 3.1 155 */ 156 public EigenDecomposition(final double[] main, final double[] secondary) { 157 isSymmetric = true; 158 this.main = main.clone(); 159 this.secondary = secondary.clone(); 160 transformer = null; 161 final int size = main.length; 162 final double[][] z = new double[size][size]; 163 for (int i = 0; i < size; i++) { 164 z[i][i] = 1.0; 165 } 166 findEigenVectors(z); 167 } 168 169 /** 170 * Calculates the eigen decomposition of the symmetric tridiagonal 171 * matrix. The Householder matrix is assumed to be the identity matrix. 172 * 173 * @param main Main diagonal of the symmetric tridiagonal form. 174 * @param secondary Secondary of the tridiagonal form. 175 * @param splitTolerance Dummy parameter (present for backward 176 * compatibility only). 177 * @throws MaxCountExceededException if the algorithm fails to converge. 178 * @deprecated in 3.1 (to be removed in 4.0) due to unused parameter 179 */ 180 @Deprecated 181 public EigenDecomposition(final double[] main, final double[] secondary, 182 final double splitTolerance) { 183 this(main, secondary); 184 } 185 186 /** 187 * Gets the matrix V of the decomposition. 188 * V is an orthogonal matrix, i.e. its transpose is also its inverse. 189 * The columns of V are the eigenvectors of the original matrix. 190 * No assumption is made about the orientation of the system axes formed 191 * by the columns of V (e.g. in a 3-dimension space, V can form a left- 192 * or right-handed system). 193 * 194 * @return the V matrix. 195 */ 196 public RealMatrix getV() { 197 198 if (cachedV == null) { 199 final int m = eigenvectors.length; 200 cachedV = MatrixUtils.createRealMatrix(m, m); 201 for (int k = 0; k < m; ++k) { 202 cachedV.setColumnVector(k, eigenvectors[k]); 203 } 204 } 205 // return the cached matrix 206 return cachedV; 207 } 208 209 /** 210 * Gets the block diagonal matrix D of the decomposition. 211 * D is a block diagonal matrix. 212 * Real eigenvalues are on the diagonal while complex values are on 213 * 2x2 blocks { {real +imaginary}, {-imaginary, real} }. 214 * 215 * @return the D matrix. 216 * 217 * @see #getRealEigenvalues() 218 * @see #getImagEigenvalues() 219 */ 220 public RealMatrix getD() { 221 222 if (cachedD == null) { 223 // cache the matrix for subsequent calls 224 cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues); 225 226 for (int i = 0; i < imagEigenvalues.length; i++) { 227 if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) > 0) { 228 cachedD.setEntry(i, i+1, imagEigenvalues[i]); 229 } else if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) { 230 cachedD.setEntry(i, i-1, imagEigenvalues[i]); 231 } 232 } 233 } 234 return cachedD; 235 } 236 237 /** 238 * Gets the transpose of the matrix V of the decomposition. 239 * V is an orthogonal matrix, i.e. its transpose is also its inverse. 240 * The columns of V are the eigenvectors of the original matrix. 241 * No assumption is made about the orientation of the system axes formed 242 * by the columns of V (e.g. in a 3-dimension space, V can form a left- 243 * or right-handed system). 244 * 245 * @return the transpose of the V matrix. 246 */ 247 public RealMatrix getVT() { 248 249 if (cachedVt == null) { 250 final int m = eigenvectors.length; 251 cachedVt = MatrixUtils.createRealMatrix(m, m); 252 for (int k = 0; k < m; ++k) { 253 cachedVt.setRowVector(k, eigenvectors[k]); 254 } 255 } 256 257 // return the cached matrix 258 return cachedVt; 259 } 260 261 /** 262 * Returns whether the calculated eigen values are complex or real. 263 * <p>The method performs a zero check for each element of the 264 * {@link #getImagEigenvalues()} array and returns {@code true} if any 265 * element is not equal to zero. 266 * 267 * @return {@code true} if the eigen values are complex, {@code false} otherwise 268 * @since 3.1 269 */ 270 public boolean hasComplexEigenvalues() { 271 for (int i = 0; i < imagEigenvalues.length; i++) { 272 if (!Precision.equals(imagEigenvalues[i], 0.0, EPSILON)) { 273 return true; 274 } 275 } 276 return false; 277 } 278 279 /** 280 * Gets a copy of the real parts of the eigenvalues of the original matrix. 281 * 282 * @return a copy of the real parts of the eigenvalues of the original matrix. 283 * 284 * @see #getD() 285 * @see #getRealEigenvalue(int) 286 * @see #getImagEigenvalues() 287 */ 288 public double[] getRealEigenvalues() { 289 return realEigenvalues.clone(); 290 } 291 292 /** 293 * Returns the real part of the i<sup>th</sup> eigenvalue of the original 294 * matrix. 295 * 296 * @param i index of the eigenvalue (counting from 0) 297 * @return real part of the i<sup>th</sup> eigenvalue of the original 298 * matrix. 299 * 300 * @see #getD() 301 * @see #getRealEigenvalues() 302 * @see #getImagEigenvalue(int) 303 */ 304 public double getRealEigenvalue(final int i) { 305 return realEigenvalues[i]; 306 } 307 308 /** 309 * Gets a copy of the imaginary parts of the eigenvalues of the original 310 * matrix. 311 * 312 * @return a copy of the imaginary parts of the eigenvalues of the original 313 * matrix. 314 * 315 * @see #getD() 316 * @see #getImagEigenvalue(int) 317 * @see #getRealEigenvalues() 318 */ 319 public double[] getImagEigenvalues() { 320 return imagEigenvalues.clone(); 321 } 322 323 /** 324 * Gets the imaginary part of the i<sup>th</sup> eigenvalue of the original 325 * matrix. 326 * 327 * @param i Index of the eigenvalue (counting from 0). 328 * @return the imaginary part of the i<sup>th</sup> eigenvalue of the original 329 * matrix. 330 * 331 * @see #getD() 332 * @see #getImagEigenvalues() 333 * @see #getRealEigenvalue(int) 334 */ 335 public double getImagEigenvalue(final int i) { 336 return imagEigenvalues[i]; 337 } 338 339 /** 340 * Gets a copy of the i<sup>th</sup> eigenvector of the original matrix. 341 * 342 * @param i Index of the eigenvector (counting from 0). 343 * @return a copy of the i<sup>th</sup> eigenvector of the original matrix. 344 * @see #getD() 345 */ 346 public RealVector getEigenvector(final int i) { 347 return eigenvectors[i].copy(); 348 } 349 350 /** 351 * Computes the determinant of the matrix. 352 * 353 * @return the determinant of the matrix. 354 */ 355 public double getDeterminant() { 356 double determinant = 1; 357 for (double lambda : realEigenvalues) { 358 determinant *= lambda; 359 } 360 return determinant; 361 } 362 363 /** 364 * Computes the square-root of the matrix. 365 * This implementation assumes that the matrix is symmetric and positive 366 * definite. 367 * 368 * @return the square-root of the matrix. 369 * @throws MathUnsupportedOperationException if the matrix is not 370 * symmetric or not positive definite. 371 * @since 3.1 372 */ 373 public RealMatrix getSquareRoot() { 374 if (!isSymmetric) { 375 throw new MathUnsupportedOperationException(); 376 } 377 378 final double[] sqrtEigenValues = new double[realEigenvalues.length]; 379 for (int i = 0; i < realEigenvalues.length; i++) { 380 final double eigen = realEigenvalues[i]; 381 if (eigen <= 0) { 382 throw new MathUnsupportedOperationException(); 383 } 384 sqrtEigenValues[i] = FastMath.sqrt(eigen); 385 } 386 final RealMatrix sqrtEigen = MatrixUtils.createRealDiagonalMatrix(sqrtEigenValues); 387 final RealMatrix v = getV(); 388 final RealMatrix vT = getVT(); 389 390 return v.multiply(sqrtEigen).multiply(vT); 391 } 392 393 /** 394 * Gets a solver for finding the A × X = B solution in exact 395 * linear sense. 396 * <p> 397 * Since 3.1, eigen decomposition of a general matrix is supported, 398 * but the {@link DecompositionSolver} only supports real eigenvalues. 399 * 400 * @return a solver 401 * @throws MathUnsupportedOperationException if the decomposition resulted in 402 * complex eigenvalues 403 */ 404 public DecompositionSolver getSolver() { 405 if (hasComplexEigenvalues()) { 406 throw new MathUnsupportedOperationException(); 407 } 408 return new Solver(realEigenvalues, imagEigenvalues, eigenvectors); 409 } 410 411 /** Specialized solver. */ 412 private static class Solver implements DecompositionSolver { 413 /** Real part of the realEigenvalues. */ 414 private double[] realEigenvalues; 415 /** Imaginary part of the realEigenvalues. */ 416 private double[] imagEigenvalues; 417 /** Eigenvectors. */ 418 private final ArrayRealVector[] eigenvectors; 419 420 /** 421 * Builds a solver from decomposed matrix. 422 * 423 * @param realEigenvalues Real parts of the eigenvalues. 424 * @param imagEigenvalues Imaginary parts of the eigenvalues. 425 * @param eigenvectors Eigenvectors. 426 */ 427 private Solver(final double[] realEigenvalues, 428 final double[] imagEigenvalues, 429 final ArrayRealVector[] eigenvectors) { 430 this.realEigenvalues = realEigenvalues; 431 this.imagEigenvalues = imagEigenvalues; 432 this.eigenvectors = eigenvectors; 433 } 434 435 /** 436 * Solves the linear equation A × X = B for symmetric matrices A. 437 * <p> 438 * This method only finds exact linear solutions, i.e. solutions for 439 * which ||A × X - B|| is exactly 0. 440 * </p> 441 * 442 * @param b Right-hand side of the equation A × X = B. 443 * @return a Vector X that minimizes the two norm of A × X - B. 444 * 445 * @throws DimensionMismatchException if the matrices dimensions do not match. 446 * @throws SingularMatrixException if the decomposed matrix is singular. 447 */ 448 public RealVector solve(final RealVector b) { 449 if (!isNonSingular()) { 450 throw new SingularMatrixException(); 451 } 452 453 final int m = realEigenvalues.length; 454 if (b.getDimension() != m) { 455 throw new DimensionMismatchException(b.getDimension(), m); 456 } 457 458 final double[] bp = new double[m]; 459 for (int i = 0; i < m; ++i) { 460 final ArrayRealVector v = eigenvectors[i]; 461 final double[] vData = v.getDataRef(); 462 final double s = v.dotProduct(b) / realEigenvalues[i]; 463 for (int j = 0; j < m; ++j) { 464 bp[j] += s * vData[j]; 465 } 466 } 467 468 return new ArrayRealVector(bp, false); 469 } 470 471 /** {@inheritDoc} */ 472 public RealMatrix solve(RealMatrix b) { 473 474 if (!isNonSingular()) { 475 throw new SingularMatrixException(); 476 } 477 478 final int m = realEigenvalues.length; 479 if (b.getRowDimension() != m) { 480 throw new DimensionMismatchException(b.getRowDimension(), m); 481 } 482 483 final int nColB = b.getColumnDimension(); 484 final double[][] bp = new double[m][nColB]; 485 final double[] tmpCol = new double[m]; 486 for (int k = 0; k < nColB; ++k) { 487 for (int i = 0; i < m; ++i) { 488 tmpCol[i] = b.getEntry(i, k); 489 bp[i][k] = 0; 490 } 491 for (int i = 0; i < m; ++i) { 492 final ArrayRealVector v = eigenvectors[i]; 493 final double[] vData = v.getDataRef(); 494 double s = 0; 495 for (int j = 0; j < m; ++j) { 496 s += v.getEntry(j) * tmpCol[j]; 497 } 498 s /= realEigenvalues[i]; 499 for (int j = 0; j < m; ++j) { 500 bp[j][k] += s * vData[j]; 501 } 502 } 503 } 504 505 return new Array2DRowRealMatrix(bp, false); 506 507 } 508 509 /** 510 * Checks whether the decomposed matrix is non-singular. 511 * 512 * @return true if the decomposed matrix is non-singular. 513 */ 514 public boolean isNonSingular() { 515 double largestEigenvalueNorm = 0.0; 516 // Looping over all values (in case they are not sorted in decreasing 517 // order of their norm). 518 for (int i = 0; i < realEigenvalues.length; ++i) { 519 largestEigenvalueNorm = FastMath.max(largestEigenvalueNorm, eigenvalueNorm(i)); 520 } 521 // Corner case: zero matrix, all exactly 0 eigenvalues 522 if (largestEigenvalueNorm == 0.0) { 523 return false; 524 } 525 for (int i = 0; i < realEigenvalues.length; ++i) { 526 // Looking for eigenvalues that are 0, where we consider anything much much smaller 527 // than the largest eigenvalue to be effectively 0. 528 if (Precision.equals(eigenvalueNorm(i) / largestEigenvalueNorm, 0, EPSILON)) { 529 return false; 530 } 531 } 532 return true; 533 } 534 535 /** 536 * @param i which eigenvalue to find the norm of 537 * @return the norm of ith (complex) eigenvalue. 538 */ 539 private double eigenvalueNorm(int i) { 540 final double re = realEigenvalues[i]; 541 final double im = imagEigenvalues[i]; 542 return FastMath.sqrt(re * re + im * im); 543 } 544 545 /** 546 * Get the inverse of the decomposed matrix. 547 * 548 * @return the inverse matrix. 549 * @throws SingularMatrixException if the decomposed matrix is singular. 550 */ 551 public RealMatrix getInverse() { 552 if (!isNonSingular()) { 553 throw new SingularMatrixException(); 554 } 555 556 final int m = realEigenvalues.length; 557 final double[][] invData = new double[m][m]; 558 559 for (int i = 0; i < m; ++i) { 560 final double[] invI = invData[i]; 561 for (int j = 0; j < m; ++j) { 562 double invIJ = 0; 563 for (int k = 0; k < m; ++k) { 564 final double[] vK = eigenvectors[k].getDataRef(); 565 invIJ += vK[i] * vK[j] / realEigenvalues[k]; 566 } 567 invI[j] = invIJ; 568 } 569 } 570 return MatrixUtils.createRealMatrix(invData); 571 } 572 } 573 574 /** 575 * Transforms the matrix to tridiagonal form. 576 * 577 * @param matrix Matrix to transform. 578 */ 579 private void transformToTridiagonal(final RealMatrix matrix) { 580 // transform the matrix to tridiagonal 581 transformer = new TriDiagonalTransformer(matrix); 582 main = transformer.getMainDiagonalRef(); 583 secondary = transformer.getSecondaryDiagonalRef(); 584 } 585 586 /** 587 * Find eigenvalues and eigenvectors (Dubrulle et al., 1971) 588 * 589 * @param householderMatrix Householder matrix of the transformation 590 * to tridiagonal form. 591 */ 592 private void findEigenVectors(final double[][] householderMatrix) { 593 final double[][]z = householderMatrix.clone(); 594 final int n = main.length; 595 realEigenvalues = new double[n]; 596 imagEigenvalues = new double[n]; 597 final double[] e = new double[n]; 598 for (int i = 0; i < n - 1; i++) { 599 realEigenvalues[i] = main[i]; 600 e[i] = secondary[i]; 601 } 602 realEigenvalues[n - 1] = main[n - 1]; 603 e[n - 1] = 0; 604 605 // Determine the largest main and secondary value in absolute term. 606 double maxAbsoluteValue = 0; 607 for (int i = 0; i < n; i++) { 608 if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) { 609 maxAbsoluteValue = FastMath.abs(realEigenvalues[i]); 610 } 611 if (FastMath.abs(e[i]) > maxAbsoluteValue) { 612 maxAbsoluteValue = FastMath.abs(e[i]); 613 } 614 } 615 // Make null any main and secondary value too small to be significant 616 if (maxAbsoluteValue != 0) { 617 for (int i=0; i < n; i++) { 618 if (FastMath.abs(realEigenvalues[i]) <= Precision.EPSILON * maxAbsoluteValue) { 619 realEigenvalues[i] = 0; 620 } 621 if (FastMath.abs(e[i]) <= Precision.EPSILON * maxAbsoluteValue) { 622 e[i]=0; 623 } 624 } 625 } 626 627 for (int j = 0; j < n; j++) { 628 int its = 0; 629 int m; 630 do { 631 for (m = j; m < n - 1; m++) { 632 double delta = FastMath.abs(realEigenvalues[m]) + 633 FastMath.abs(realEigenvalues[m + 1]); 634 if (FastMath.abs(e[m]) + delta == delta) { 635 break; 636 } 637 } 638 if (m != j) { 639 if (its == maxIter) { 640 throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED, 641 maxIter); 642 } 643 its++; 644 double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]); 645 double t = FastMath.sqrt(1 + q * q); 646 if (q < 0.0) { 647 q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q - t); 648 } else { 649 q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q + t); 650 } 651 double u = 0.0; 652 double s = 1.0; 653 double c = 1.0; 654 int i; 655 for (i = m - 1; i >= j; i--) { 656 double p = s * e[i]; 657 double h = c * e[i]; 658 if (FastMath.abs(p) >= FastMath.abs(q)) { 659 c = q / p; 660 t = FastMath.sqrt(c * c + 1.0); 661 e[i + 1] = p * t; 662 s = 1.0 / t; 663 c *= s; 664 } else { 665 s = p / q; 666 t = FastMath.sqrt(s * s + 1.0); 667 e[i + 1] = q * t; 668 c = 1.0 / t; 669 s *= c; 670 } 671 if (e[i + 1] == 0.0) { 672 realEigenvalues[i + 1] -= u; 673 e[m] = 0.0; 674 break; 675 } 676 q = realEigenvalues[i + 1] - u; 677 t = (realEigenvalues[i] - q) * s + 2.0 * c * h; 678 u = s * t; 679 realEigenvalues[i + 1] = q + u; 680 q = c * t - h; 681 for (int ia = 0; ia < n; ia++) { 682 p = z[ia][i + 1]; 683 z[ia][i + 1] = s * z[ia][i] + c * p; 684 z[ia][i] = c * z[ia][i] - s * p; 685 } 686 } 687 if (t == 0.0 && i >= j) { 688 continue; 689 } 690 realEigenvalues[j] -= u; 691 e[j] = q; 692 e[m] = 0.0; 693 } 694 } while (m != j); 695 } 696 697 //Sort the eigen values (and vectors) in increase order 698 for (int i = 0; i < n; i++) { 699 int k = i; 700 double p = realEigenvalues[i]; 701 for (int j = i + 1; j < n; j++) { 702 if (realEigenvalues[j] > p) { 703 k = j; 704 p = realEigenvalues[j]; 705 } 706 } 707 if (k != i) { 708 realEigenvalues[k] = realEigenvalues[i]; 709 realEigenvalues[i] = p; 710 for (int j = 0; j < n; j++) { 711 p = z[j][i]; 712 z[j][i] = z[j][k]; 713 z[j][k] = p; 714 } 715 } 716 } 717 718 // Determine the largest eigen value in absolute term. 719 maxAbsoluteValue = 0; 720 for (int i = 0; i < n; i++) { 721 if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) { 722 maxAbsoluteValue=FastMath.abs(realEigenvalues[i]); 723 } 724 } 725 // Make null any eigen value too small to be significant 726 if (maxAbsoluteValue != 0.0) { 727 for (int i=0; i < n; i++) { 728 if (FastMath.abs(realEigenvalues[i]) < Precision.EPSILON * maxAbsoluteValue) { 729 realEigenvalues[i] = 0; 730 } 731 } 732 } 733 eigenvectors = new ArrayRealVector[n]; 734 final double[] tmp = new double[n]; 735 for (int i = 0; i < n; i++) { 736 for (int j = 0; j < n; j++) { 737 tmp[j] = z[j][i]; 738 } 739 eigenvectors[i] = new ArrayRealVector(tmp); 740 } 741 } 742 743 /** 744 * Transforms the matrix to Schur form and calculates the eigenvalues. 745 * 746 * @param matrix Matrix to transform. 747 * @return the {@link SchurTransformer Shur transform} for this matrix 748 */ 749 private SchurTransformer transformToSchur(final RealMatrix matrix) { 750 final SchurTransformer schurTransform = new SchurTransformer(matrix); 751 final double[][] matT = schurTransform.getT().getData(); 752 753 realEigenvalues = new double[matT.length]; 754 imagEigenvalues = new double[matT.length]; 755 756 for (int i = 0; i < realEigenvalues.length; i++) { 757 if (i == (realEigenvalues.length - 1) || 758 Precision.equals(matT[i + 1][i], 0.0, EPSILON)) { 759 realEigenvalues[i] = matT[i][i]; 760 } else { 761 final double x = matT[i + 1][i + 1]; 762 final double p = 0.5 * (matT[i][i] - x); 763 final double z = FastMath.sqrt(FastMath.abs(p * p + matT[i + 1][i] * matT[i][i + 1])); 764 realEigenvalues[i] = x + p; 765 imagEigenvalues[i] = z; 766 realEigenvalues[i + 1] = x + p; 767 imagEigenvalues[i + 1] = -z; 768 i++; 769 } 770 } 771 return schurTransform; 772 } 773 774 /** 775 * Performs a division of two complex numbers. 776 * 777 * @param xr real part of the first number 778 * @param xi imaginary part of the first number 779 * @param yr real part of the second number 780 * @param yi imaginary part of the second number 781 * @return result of the complex division 782 */ 783 private Complex cdiv(final double xr, final double xi, 784 final double yr, final double yi) { 785 return new Complex(xr, xi).divide(new Complex(yr, yi)); 786 } 787 788 /** 789 * Find eigenvectors from a matrix transformed to Schur form. 790 * 791 * @param schur the schur transformation of the matrix 792 * @throws MathArithmeticException if the Schur form has a norm of zero 793 */ 794 private void findEigenVectorsFromSchur(final SchurTransformer schur) 795 throws MathArithmeticException { 796 final double[][] matrixT = schur.getT().getData(); 797 final double[][] matrixP = schur.getP().getData(); 798 799 final int n = matrixT.length; 800 801 // compute matrix norm 802 double norm = 0.0; 803 for (int i = 0; i < n; i++) { 804 for (int j = FastMath.max(i - 1, 0); j < n; j++) { 805 norm += FastMath.abs(matrixT[i][j]); 806 } 807 } 808 809 // we can not handle a matrix with zero norm 810 if (Precision.equals(norm, 0.0, EPSILON)) { 811 throw new MathArithmeticException(LocalizedFormats.ZERO_NORM); 812 } 813 814 // Backsubstitute to find vectors of upper triangular form 815 816 double r = 0.0; 817 double s = 0.0; 818 double z = 0.0; 819 820 for (int idx = n - 1; idx >= 0; idx--) { 821 double p = realEigenvalues[idx]; 822 double q = imagEigenvalues[idx]; 823 824 if (Precision.equals(q, 0.0)) { 825 // Real vector 826 int l = idx; 827 matrixT[idx][idx] = 1.0; 828 for (int i = idx - 1; i >= 0; i--) { 829 double w = matrixT[i][i] - p; 830 r = 0.0; 831 for (int j = l; j <= idx; j++) { 832 r += matrixT[i][j] * matrixT[j][idx]; 833 } 834 if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) { 835 z = w; 836 s = r; 837 } else { 838 l = i; 839 if (Precision.equals(imagEigenvalues[i], 0.0)) { 840 if (w != 0.0) { 841 matrixT[i][idx] = -r / w; 842 } else { 843 matrixT[i][idx] = -r / (Precision.EPSILON * norm); 844 } 845 } else { 846 // Solve real equations 847 double x = matrixT[i][i + 1]; 848 double y = matrixT[i + 1][i]; 849 q = (realEigenvalues[i] - p) * (realEigenvalues[i] - p) + 850 imagEigenvalues[i] * imagEigenvalues[i]; 851 double t = (x * s - z * r) / q; 852 matrixT[i][idx] = t; 853 if (FastMath.abs(x) > FastMath.abs(z)) { 854 matrixT[i + 1][idx] = (-r - w * t) / x; 855 } else { 856 matrixT[i + 1][idx] = (-s - y * t) / z; 857 } 858 } 859 860 // Overflow control 861 double t = FastMath.abs(matrixT[i][idx]); 862 if ((Precision.EPSILON * t) * t > 1) { 863 for (int j = i; j <= idx; j++) { 864 matrixT[j][idx] /= t; 865 } 866 } 867 } 868 } 869 } else if (q < 0.0) { 870 // Complex vector 871 int l = idx - 1; 872 873 // Last vector component imaginary so matrix is triangular 874 if (FastMath.abs(matrixT[idx][idx - 1]) > FastMath.abs(matrixT[idx - 1][idx])) { 875 matrixT[idx - 1][idx - 1] = q / matrixT[idx][idx - 1]; 876 matrixT[idx - 1][idx] = -(matrixT[idx][idx] - p) / matrixT[idx][idx - 1]; 877 } else { 878 final Complex result = cdiv(0.0, -matrixT[idx - 1][idx], 879 matrixT[idx - 1][idx - 1] - p, q); 880 matrixT[idx - 1][idx - 1] = result.getReal(); 881 matrixT[idx - 1][idx] = result.getImaginary(); 882 } 883 884 matrixT[idx][idx - 1] = 0.0; 885 matrixT[idx][idx] = 1.0; 886 887 for (int i = idx - 2; i >= 0; i--) { 888 double ra = 0.0; 889 double sa = 0.0; 890 for (int j = l; j <= idx; j++) { 891 ra += matrixT[i][j] * matrixT[j][idx - 1]; 892 sa += matrixT[i][j] * matrixT[j][idx]; 893 } 894 double w = matrixT[i][i] - p; 895 896 if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) { 897 z = w; 898 r = ra; 899 s = sa; 900 } else { 901 l = i; 902 if (Precision.equals(imagEigenvalues[i], 0.0)) { 903 final Complex c = cdiv(-ra, -sa, w, q); 904 matrixT[i][idx - 1] = c.getReal(); 905 matrixT[i][idx] = c.getImaginary(); 906 } else { 907 // Solve complex equations 908 double x = matrixT[i][i + 1]; 909 double y = matrixT[i + 1][i]; 910 double vr = (realEigenvalues[i] - p) * (realEigenvalues[i] - p) + 911 imagEigenvalues[i] * imagEigenvalues[i] - q * q; 912 final double vi = (realEigenvalues[i] - p) * 2.0 * q; 913 if (Precision.equals(vr, 0.0) && Precision.equals(vi, 0.0)) { 914 vr = Precision.EPSILON * norm * 915 (FastMath.abs(w) + FastMath.abs(q) + FastMath.abs(x) + 916 FastMath.abs(y) + FastMath.abs(z)); 917 } 918 final Complex c = cdiv(x * r - z * ra + q * sa, 919 x * s - z * sa - q * ra, vr, vi); 920 matrixT[i][idx - 1] = c.getReal(); 921 matrixT[i][idx] = c.getImaginary(); 922 923 if (FastMath.abs(x) > (FastMath.abs(z) + FastMath.abs(q))) { 924 matrixT[i + 1][idx - 1] = (-ra - w * matrixT[i][idx - 1] + 925 q * matrixT[i][idx]) / x; 926 matrixT[i + 1][idx] = (-sa - w * matrixT[i][idx] - 927 q * matrixT[i][idx - 1]) / x; 928 } else { 929 final Complex c2 = cdiv(-r - y * matrixT[i][idx - 1], 930 -s - y * matrixT[i][idx], z, q); 931 matrixT[i + 1][idx - 1] = c2.getReal(); 932 matrixT[i + 1][idx] = c2.getImaginary(); 933 } 934 } 935 936 // Overflow control 937 double t = FastMath.max(FastMath.abs(matrixT[i][idx - 1]), 938 FastMath.abs(matrixT[i][idx])); 939 if ((Precision.EPSILON * t) * t > 1) { 940 for (int j = i; j <= idx; j++) { 941 matrixT[j][idx - 1] /= t; 942 matrixT[j][idx] /= t; 943 } 944 } 945 } 946 } 947 } 948 } 949 950 // Back transformation to get eigenvectors of original matrix 951 for (int j = n - 1; j >= 0; j--) { 952 for (int i = 0; i <= n - 1; i++) { 953 z = 0.0; 954 for (int k = 0; k <= FastMath.min(j, n - 1); k++) { 955 z += matrixP[i][k] * matrixT[k][j]; 956 } 957 matrixP[i][j] = z; 958 } 959 } 960 961 eigenvectors = new ArrayRealVector[n]; 962 final double[] tmp = new double[n]; 963 for (int i = 0; i < n; i++) { 964 for (int j = 0; j < n; j++) { 965 tmp[j] = matrixP[j][i]; 966 } 967 eigenvectors[i] = new ArrayRealVector(tmp); 968 } 969 } 970}