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