1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 package org.apache.commons.math3.linear;
19
20 import org.apache.commons.math3.complex.Complex;
21 import org.apache.commons.math3.exception.MathArithmeticException;
22 import org.apache.commons.math3.exception.MathUnsupportedOperationException;
23 import org.apache.commons.math3.exception.MaxCountExceededException;
24 import org.apache.commons.math3.exception.DimensionMismatchException;
25 import org.apache.commons.math3.exception.util.LocalizedFormats;
26 import org.apache.commons.math3.util.Precision;
27 import org.apache.commons.math3.util.FastMath;
28
29 /**
30 * Calculates the eigen decomposition of a real matrix.
31 * <p>The eigen decomposition of matrix A is a set of two matrices:
32 * V and D such that A = V × D × V<sup>T</sup>.
33 * A, V and D are all m × m matrices.</p>
34 * <p>This class is similar in spirit to the <code>EigenvalueDecomposition</code>
35 * class from the <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a>
36 * library, with the following changes:</p>
37 * <ul>
38 * <li>a {@link #getVT() getVt} method has been added,</li>
39 * <li>two {@link #getRealEigenvalue(int) getRealEigenvalue} and {@link #getImagEigenvalue(int)
40 * getImagEigenvalue} methods to pick up a single eigenvalue have been added,</li>
41 * <li>a {@link #getEigenvector(int) getEigenvector} method to pick up a single
42 * eigenvector has been added,</li>
43 * <li>a {@link #getDeterminant() getDeterminant} method has been added.</li>
44 * <li>a {@link #getSolver() getSolver} method has been added.</li>
45 * </ul>
46 * <p>
47 * As of 3.1, this class supports general real matrices (both symmetric and non-symmetric):
48 * </p>
49 * <p>
50 * If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is diagonal and the eigenvector
51 * matrix V is orthogonal, i.e. A = V.multiply(D.multiply(V.transpose())) and
52 * V.multiply(V.transpose()) equals the identity matrix.
53 * </p>
54 * <p>
55 * If A is not symmetric, then the eigenvalue matrix D is block diagonal with the real eigenvalues
56 * in 1-by-1 blocks and any complex eigenvalues, lambda + i*mu, in 2-by-2 blocks:
57 * <pre>
58 * [lambda, mu ]
59 * [ -mu, lambda]
60 * </pre>
61 * The columns of V represent the eigenvectors in the sense that A*V = V*D,
62 * i.e. A.multiply(V) equals V.multiply(D).
63 * The matrix V may be badly conditioned, or even singular, so the validity of the equation
64 * A = V*D*inverse(V) depends upon the condition of V.
65 * </p>
66 * <p>
67 * This implementation is based on the paper by A. Drubrulle, R.S. Martin and
68 * J.H. Wilkinson "The Implicit QL Algorithm" in Wilksinson and Reinsch (1971)
69 * Handbook for automatic computation, vol. 2, Linear algebra, Springer-Verlag,
70 * New-York
71 * </p>
72 * @see <a href="http://mathworld.wolfram.com/EigenDecomposition.html">MathWorld</a>
73 * @see <a href="http://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix">Wikipedia</a>
74 * @version $Id: EigenDecomposition.java 1452595 2013-03-04 23:29:39Z tn $
75 * @since 2.0 (changed to concrete class in 3.0)
76 */
77 public class EigenDecomposition {
78 /** Internally used epsilon criteria. */
79 private static final double EPSILON = 1e-12;
80 /** Maximum number of iterations accepted in the implicit QL transformation */
81 private byte maxIter = 30;
82 /** Main diagonal of the tridiagonal matrix. */
83 private double[] main;
84 /** Secondary diagonal of the tridiagonal matrix. */
85 private double[] secondary;
86 /**
87 * Transformer to tridiagonal (may be null if matrix is already
88 * tridiagonal).
89 */
90 private TriDiagonalTransformer transformer;
91 /** Real part of the realEigenvalues. */
92 private double[] realEigenvalues;
93 /** Imaginary part of the realEigenvalues. */
94 private double[] imagEigenvalues;
95 /** Eigenvectors. */
96 private ArrayRealVector[] eigenvectors;
97 /** Cached value of V. */
98 private RealMatrix cachedV;
99 /** 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 positive
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 }