View Javadoc
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 &times; D &times; V<sup>T</sup>.
33   * A, V and D are all m &times; 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 1538368 2013-11-03 13:57:37Z erans $
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 &times; 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 &times; X = B for symmetric matrices A.
438          * <p>
439          * This method only finds exact linear solutions, i.e. solutions for
440          * which ||A &times; X - B|| is exactly 0.
441          * </p>
442          *
443          * @param b Right-hand side of the equation A &times; X = B.
444          * @return a Vector X that minimizes the two norm of A &times; 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             double largestEigenvalueNorm = 0.0;
517             // Looping over all values (in case they are not sorted in decreasing
518             // order of their norm).
519             for (int i = 0; i < realEigenvalues.length; ++i) {
520                 largestEigenvalueNorm = FastMath.max(largestEigenvalueNorm, eigenvalueNorm(i));
521             }
522             // Corner case: zero matrix, all exactly 0 eigenvalues
523             if (largestEigenvalueNorm == 0.0) {
524                 return false;
525             }
526             for (int i = 0; i < realEigenvalues.length; ++i) {
527                 // Looking for eigenvalues that are 0, where we consider anything much much smaller
528                 // than the largest eigenvalue to be effectively 0.
529                 if (Precision.equals(eigenvalueNorm(i) / largestEigenvalueNorm, 0, EPSILON)) {
530                     return false;
531                 }
532             }
533             return true;
534         }
535 
536         /**
537          * @param i which eigenvalue to find the norm of
538          * @return the norm of ith (complex) eigenvalue.
539          */
540         private double eigenvalueNorm(int i) {
541             final double re = realEigenvalues[i];
542             final double im = imagEigenvalues[i];
543             return FastMath.sqrt(re * re + im * im);
544         }
545 
546         /**
547          * Get the inverse of the decomposed matrix.
548          *
549          * @return the inverse matrix.
550          * @throws SingularMatrixException if the decomposed matrix is singular.
551          */
552         public RealMatrix getInverse() {
553             if (!isNonSingular()) {
554                 throw new SingularMatrixException();
555             }
556 
557             final int m = realEigenvalues.length;
558             final double[][] invData = new double[m][m];
559 
560             for (int i = 0; i < m; ++i) {
561                 final double[] invI = invData[i];
562                 for (int j = 0; j < m; ++j) {
563                     double invIJ = 0;
564                     for (int k = 0; k < m; ++k) {
565                         final double[] vK = eigenvectors[k].getDataRef();
566                         invIJ += vK[i] * vK[j] / realEigenvalues[k];
567                     }
568                     invI[j] = invIJ;
569                 }
570             }
571             return MatrixUtils.createRealMatrix(invData);
572         }
573     }
574 
575     /**
576      * Transforms the matrix to tridiagonal form.
577      *
578      * @param matrix Matrix to transform.
579      */
580     private void transformToTridiagonal(final RealMatrix matrix) {
581         // transform the matrix to tridiagonal
582         transformer = new TriDiagonalTransformer(matrix);
583         main = transformer.getMainDiagonalRef();
584         secondary = transformer.getSecondaryDiagonalRef();
585     }
586 
587     /**
588      * Find eigenvalues and eigenvectors (Dubrulle et al., 1971)
589      *
590      * @param householderMatrix Householder matrix of the transformation
591      * to tridiagonal form.
592      */
593     private void findEigenVectors(final double[][] householderMatrix) {
594         final double[][]z = householderMatrix.clone();
595         final int n = main.length;
596         realEigenvalues = new double[n];
597         imagEigenvalues = new double[n];
598         final double[] e = new double[n];
599         for (int i = 0; i < n - 1; i++) {
600             realEigenvalues[i] = main[i];
601             e[i] = secondary[i];
602         }
603         realEigenvalues[n - 1] = main[n - 1];
604         e[n - 1] = 0;
605 
606         // Determine the largest main and secondary value in absolute term.
607         double maxAbsoluteValue = 0;
608         for (int i = 0; i < n; i++) {
609             if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) {
610                 maxAbsoluteValue = FastMath.abs(realEigenvalues[i]);
611             }
612             if (FastMath.abs(e[i]) > maxAbsoluteValue) {
613                 maxAbsoluteValue = FastMath.abs(e[i]);
614             }
615         }
616         // Make null any main and secondary value too small to be significant
617         if (maxAbsoluteValue != 0) {
618             for (int i=0; i < n; i++) {
619                 if (FastMath.abs(realEigenvalues[i]) <= Precision.EPSILON * maxAbsoluteValue) {
620                     realEigenvalues[i] = 0;
621                 }
622                 if (FastMath.abs(e[i]) <= Precision.EPSILON * maxAbsoluteValue) {
623                     e[i]=0;
624                 }
625             }
626         }
627 
628         for (int j = 0; j < n; j++) {
629             int its = 0;
630             int m;
631             do {
632                 for (m = j; m < n - 1; m++) {
633                     double delta = FastMath.abs(realEigenvalues[m]) +
634                         FastMath.abs(realEigenvalues[m + 1]);
635                     if (FastMath.abs(e[m]) + delta == delta) {
636                         break;
637                     }
638                 }
639                 if (m != j) {
640                     if (its == maxIter) {
641                         throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED,
642                                                             maxIter);
643                     }
644                     its++;
645                     double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]);
646                     double t = FastMath.sqrt(1 + q * q);
647                     if (q < 0.0) {
648                         q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q - t);
649                     } else {
650                         q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q + t);
651                     }
652                     double u = 0.0;
653                     double s = 1.0;
654                     double c = 1.0;
655                     int i;
656                     for (i = m - 1; i >= j; i--) {
657                         double p = s * e[i];
658                         double h = c * e[i];
659                         if (FastMath.abs(p) >= FastMath.abs(q)) {
660                             c = q / p;
661                             t = FastMath.sqrt(c * c + 1.0);
662                             e[i + 1] = p * t;
663                             s = 1.0 / t;
664                             c *= s;
665                         } else {
666                             s = p / q;
667                             t = FastMath.sqrt(s * s + 1.0);
668                             e[i + 1] = q * t;
669                             c = 1.0 / t;
670                             s *= c;
671                         }
672                         if (e[i + 1] == 0.0) {
673                             realEigenvalues[i + 1] -= u;
674                             e[m] = 0.0;
675                             break;
676                         }
677                         q = realEigenvalues[i + 1] - u;
678                         t = (realEigenvalues[i] - q) * s + 2.0 * c * h;
679                         u = s * t;
680                         realEigenvalues[i + 1] = q + u;
681                         q = c * t - h;
682                         for (int ia = 0; ia < n; ia++) {
683                             p = z[ia][i + 1];
684                             z[ia][i + 1] = s * z[ia][i] + c * p;
685                             z[ia][i] = c * z[ia][i] - s * p;
686                         }
687                     }
688                     if (t == 0.0 && i >= j) {
689                         continue;
690                     }
691                     realEigenvalues[j] -= u;
692                     e[j] = q;
693                     e[m] = 0.0;
694                 }
695             } while (m != j);
696         }
697 
698         //Sort the eigen values (and vectors) in increase order
699         for (int i = 0; i < n; i++) {
700             int k = i;
701             double p = realEigenvalues[i];
702             for (int j = i + 1; j < n; j++) {
703                 if (realEigenvalues[j] > p) {
704                     k = j;
705                     p = realEigenvalues[j];
706                 }
707             }
708             if (k != i) {
709                 realEigenvalues[k] = realEigenvalues[i];
710                 realEigenvalues[i] = p;
711                 for (int j = 0; j < n; j++) {
712                     p = z[j][i];
713                     z[j][i] = z[j][k];
714                     z[j][k] = p;
715                 }
716             }
717         }
718 
719         // Determine the largest eigen value in absolute term.
720         maxAbsoluteValue = 0;
721         for (int i = 0; i < n; i++) {
722             if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) {
723                 maxAbsoluteValue=FastMath.abs(realEigenvalues[i]);
724             }
725         }
726         // Make null any eigen value too small to be significant
727         if (maxAbsoluteValue != 0.0) {
728             for (int i=0; i < n; i++) {
729                 if (FastMath.abs(realEigenvalues[i]) < Precision.EPSILON * maxAbsoluteValue) {
730                     realEigenvalues[i] = 0;
731                 }
732             }
733         }
734         eigenvectors = new ArrayRealVector[n];
735         final double[] tmp = new double[n];
736         for (int i = 0; i < n; i++) {
737             for (int j = 0; j < n; j++) {
738                 tmp[j] = z[j][i];
739             }
740             eigenvectors[i] = new ArrayRealVector(tmp);
741         }
742     }
743 
744     /**
745      * Transforms the matrix to Schur form and calculates the eigenvalues.
746      *
747      * @param matrix Matrix to transform.
748      * @return the {@link SchurTransformer Shur transform} for this matrix
749      */
750     private SchurTransformer transformToSchur(final RealMatrix matrix) {
751         final SchurTransformer schurTransform = new SchurTransformer(matrix);
752         final double[][] matT = schurTransform.getT().getData();
753 
754         realEigenvalues = new double[matT.length];
755         imagEigenvalues = new double[matT.length];
756 
757         for (int i = 0; i < realEigenvalues.length; i++) {
758             if (i == (realEigenvalues.length - 1) ||
759                 Precision.equals(matT[i + 1][i], 0.0, EPSILON)) {
760                 realEigenvalues[i] = matT[i][i];
761             } else {
762                 final double x = matT[i + 1][i + 1];
763                 final double p = 0.5 * (matT[i][i] - x);
764                 final double z = FastMath.sqrt(FastMath.abs(p * p + matT[i + 1][i] * matT[i][i + 1]));
765                 realEigenvalues[i] = x + p;
766                 imagEigenvalues[i] = z;
767                 realEigenvalues[i + 1] = x + p;
768                 imagEigenvalues[i + 1] = -z;
769                 i++;
770             }
771         }
772         return schurTransform;
773     }
774 
775     /**
776      * Performs a division of two complex numbers.
777      *
778      * @param xr real part of the first number
779      * @param xi imaginary part of the first number
780      * @param yr real part of the second number
781      * @param yi imaginary part of the second number
782      * @return result of the complex division
783      */
784     private Complex cdiv(final double xr, final double xi,
785                          final double yr, final double yi) {
786         return new Complex(xr, xi).divide(new Complex(yr, yi));
787     }
788 
789     /**
790      * Find eigenvectors from a matrix transformed to Schur form.
791      *
792      * @param schur the schur transformation of the matrix
793      * @throws MathArithmeticException if the Schur form has a norm of zero
794      */
795     private void findEigenVectorsFromSchur(final SchurTransformer schur)
796         throws MathArithmeticException {
797         final double[][] matrixT = schur.getT().getData();
798         final double[][] matrixP = schur.getP().getData();
799 
800         final int n = matrixT.length;
801 
802         // compute matrix norm
803         double norm = 0.0;
804         for (int i = 0; i < n; i++) {
805            for (int j = FastMath.max(i - 1, 0); j < n; j++) {
806                norm += FastMath.abs(matrixT[i][j]);
807            }
808         }
809 
810         // we can not handle a matrix with zero norm
811         if (Precision.equals(norm, 0.0, EPSILON)) {
812            throw new MathArithmeticException(LocalizedFormats.ZERO_NORM);
813         }
814 
815         // Backsubstitute to find vectors of upper triangular form
816 
817         double r = 0.0;
818         double s = 0.0;
819         double z = 0.0;
820 
821         for (int idx = n - 1; idx >= 0; idx--) {
822             double p = realEigenvalues[idx];
823             double q = imagEigenvalues[idx];
824 
825             if (Precision.equals(q, 0.0)) {
826                 // Real vector
827                 int l = idx;
828                 matrixT[idx][idx] = 1.0;
829                 for (int i = idx - 1; i >= 0; i--) {
830                     double w = matrixT[i][i] - p;
831                     r = 0.0;
832                     for (int j = l; j <= idx; j++) {
833                         r += matrixT[i][j] * matrixT[j][idx];
834                     }
835                     if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) {
836                         z = w;
837                         s = r;
838                     } else {
839                         l = i;
840                         if (Precision.equals(imagEigenvalues[i], 0.0)) {
841                             if (w != 0.0) {
842                                 matrixT[i][idx] = -r / w;
843                             } else {
844                                 matrixT[i][idx] = -r / (Precision.EPSILON * norm);
845                             }
846                         } else {
847                             // Solve real equations
848                             double x = matrixT[i][i + 1];
849                             double y = matrixT[i + 1][i];
850                             q = (realEigenvalues[i] - p) * (realEigenvalues[i] - p) +
851                                 imagEigenvalues[i] * imagEigenvalues[i];
852                             double t = (x * s - z * r) / q;
853                             matrixT[i][idx] = t;
854                             if (FastMath.abs(x) > FastMath.abs(z)) {
855                                 matrixT[i + 1][idx] = (-r - w * t) / x;
856                             } else {
857                                 matrixT[i + 1][idx] = (-s - y * t) / z;
858                             }
859                         }
860 
861                         // Overflow control
862                         double t = FastMath.abs(matrixT[i][idx]);
863                         if ((Precision.EPSILON * t) * t > 1) {
864                             for (int j = i; j <= idx; j++) {
865                                 matrixT[j][idx] /= t;
866                             }
867                         }
868                     }
869                 }
870             } else if (q < 0.0) {
871                 // Complex vector
872                 int l = idx - 1;
873 
874                 // Last vector component imaginary so matrix is triangular
875                 if (FastMath.abs(matrixT[idx][idx - 1]) > FastMath.abs(matrixT[idx - 1][idx])) {
876                     matrixT[idx - 1][idx - 1] = q / matrixT[idx][idx - 1];
877                     matrixT[idx - 1][idx]     = -(matrixT[idx][idx] - p) / matrixT[idx][idx - 1];
878                 } else {
879                     final Complex result = cdiv(0.0, -matrixT[idx - 1][idx],
880                                                 matrixT[idx - 1][idx - 1] - p, q);
881                     matrixT[idx - 1][idx - 1] = result.getReal();
882                     matrixT[idx - 1][idx]     = result.getImaginary();
883                 }
884 
885                 matrixT[idx][idx - 1] = 0.0;
886                 matrixT[idx][idx]     = 1.0;
887 
888                 for (int i = idx - 2; i >= 0; i--) {
889                     double ra = 0.0;
890                     double sa = 0.0;
891                     for (int j = l; j <= idx; j++) {
892                         ra += matrixT[i][j] * matrixT[j][idx - 1];
893                         sa += matrixT[i][j] * matrixT[j][idx];
894                     }
895                     double w = matrixT[i][i] - p;
896 
897                     if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) {
898                         z = w;
899                         r = ra;
900                         s = sa;
901                     } else {
902                         l = i;
903                         if (Precision.equals(imagEigenvalues[i], 0.0)) {
904                             final Complex c = cdiv(-ra, -sa, w, q);
905                             matrixT[i][idx - 1] = c.getReal();
906                             matrixT[i][idx] = c.getImaginary();
907                         } else {
908                             // Solve complex equations
909                             double x = matrixT[i][i + 1];
910                             double y = matrixT[i + 1][i];
911                             double vr = (realEigenvalues[i] - p) * (realEigenvalues[i] - p) +
912                                         imagEigenvalues[i] * imagEigenvalues[i] - q * q;
913                             final double vi = (realEigenvalues[i] - p) * 2.0 * q;
914                             if (Precision.equals(vr, 0.0) && Precision.equals(vi, 0.0)) {
915                                 vr = Precision.EPSILON * norm *
916                                      (FastMath.abs(w) + FastMath.abs(q) + FastMath.abs(x) +
917                                       FastMath.abs(y) + FastMath.abs(z));
918                             }
919                             final Complex c     = cdiv(x * r - z * ra + q * sa,
920                                                        x * s - z * sa - q * ra, vr, vi);
921                             matrixT[i][idx - 1] = c.getReal();
922                             matrixT[i][idx]     = c.getImaginary();
923 
924                             if (FastMath.abs(x) > (FastMath.abs(z) + FastMath.abs(q))) {
925                                 matrixT[i + 1][idx - 1] = (-ra - w * matrixT[i][idx - 1] +
926                                                            q * matrixT[i][idx]) / x;
927                                 matrixT[i + 1][idx]     = (-sa - w * matrixT[i][idx] -
928                                                            q * matrixT[i][idx - 1]) / x;
929                             } else {
930                                 final Complex c2        = cdiv(-r - y * matrixT[i][idx - 1],
931                                                                -s - y * matrixT[i][idx], z, q);
932                                 matrixT[i + 1][idx - 1] = c2.getReal();
933                                 matrixT[i + 1][idx]     = c2.getImaginary();
934                             }
935                         }
936 
937                         // Overflow control
938                         double t = FastMath.max(FastMath.abs(matrixT[i][idx - 1]),
939                                                 FastMath.abs(matrixT[i][idx]));
940                         if ((Precision.EPSILON * t) * t > 1) {
941                             for (int j = i; j <= idx; j++) {
942                                 matrixT[j][idx - 1] /= t;
943                                 matrixT[j][idx] /= t;
944                             }
945                         }
946                     }
947                 }
948             }
949         }
950 
951         // Back transformation to get eigenvectors of original matrix
952         for (int j = n - 1; j >= 0; j--) {
953             for (int i = 0; i <= n - 1; i++) {
954                 z = 0.0;
955                 for (int k = 0; k <= FastMath.min(j, n - 1); k++) {
956                     z += matrixP[i][k] * matrixT[k][j];
957                 }
958                 matrixP[i][j] = z;
959             }
960         }
961 
962         eigenvectors = new ArrayRealVector[n];
963         final double[] tmp = new double[n];
964         for (int i = 0; i < n; i++) {
965             for (int j = 0; j < n; j++) {
966                 tmp[j] = matrixP[j][i];
967             }
968             eigenvectors[i] = new ArrayRealVector(tmp);
969         }
970     }
971 }