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