EigenDecomposition.java

  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.  */

  18. package org.apache.commons.math4.legacy.linear;

  20. import org.apache.commons.numbers.complex.Complex;
  21. import org.apache.commons.numbers.core.Precision;
  22. import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
  23. import org.apache.commons.math4.legacy.exception.MathArithmeticException;
  24. import org.apache.commons.math4.legacy.exception.MathUnsupportedOperationException;
  25. import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
  26. import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
  27. import org.apache.commons.math4.core.jdkmath.JdkMath;

  29. /**
  30.  * Calculates the eigen decomposition of a real matrix.
  31.  * <p>
  32.  * The eigen decomposition of matrix A is a set of two matrices:
  33.  * V and D such that A = V &times; D &times; V<sup>T</sup>.
  34.  * A, V and D are all m &times; m matrices.
  35.  * <p>
  36.  * This class is similar in spirit to the {@code EigenvalueDecomposition}
  37.  * class from the <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a>
  38.  * library, with the following changes:
  39.  * <ul>
  40.  *   <li>a {@link #getVT() getVt} method has been added,</li>
  41.  *   <li>two {@link #getRealEigenvalue(int) getRealEigenvalue} and
  42.  *       {@link #getImagEigenvalue(int) getImagEigenvalue} methods to pick up a
  43.  *       single eigenvalue have been added,</li>
  44.  *   <li>a {@link #getEigenvector(int) getEigenvector} method to pick up a
  45.  *       single eigenvector has been added,</li>
  46.  *   <li>a {@link #getDeterminant() getDeterminant} method has been added.</li>
  47.  *   <li>a {@link #getSolver() getSolver} method has been added.</li>
  48.  * </ul>
  49.  * <p>
  50.  * As of 3.1, this class supports general real matrices (both symmetric and non-symmetric):
  51.  * <p>
  52.  * If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is diagonal
  53.  * and the eigenvector matrix V is orthogonal, i.e.
  54.  * {@code A = V.multiply(D.multiply(V.transpose()))} and
  55.  * {@code V.multiply(V.transpose())} equals the identity matrix.
  56.  * </p>
  57.  * <p>
  58.  * If A is not symmetric, then the eigenvalue matrix D is block diagonal with the real
  59.  * eigenvalues in 1-by-1 blocks and any complex eigenvalues, lambda + i*mu, in 2-by-2
  60.  * blocks:
  61.  * <pre>
  62.  *    [lambda, mu    ]
  63.  *    [   -mu, lambda]
  64.  * </pre>
  65.  * The columns of V represent the eigenvectors in the sense that {@code A*V = V*D},
  66.  * i.e. A.multiply(V) equals V.multiply(D).
  67.  * The matrix V may be badly conditioned, or even singular, so the validity of the
  68.  * equation {@code A = V*D*inverse(V)} depends upon the condition of V.
  69.  * <p>
  70.  * This implementation is based on the paper by A. Drubrulle, R.S. Martin and
  71.  * J.H. Wilkinson "The Implicit QL Algorithm" in Wilksinson and Reinsch (1971)
  72.  * Handbook for automatic computation, vol. 2, Linear algebra, Springer-Verlag,
  73.  * New-York.
  74.  *
  75.  * @see <a href="http://mathworld.wolfram.com/EigenDecomposition.html">MathWorld</a>
  76.  * @see <a href="http://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix">Wikipedia</a>
  77.  * @since 2.0 (changed to concrete class in 3.0)
  78.  */
  79. public class EigenDecomposition {
  80.     /** Internally used epsilon criteria. */
  81.     private static final double EPSILON = 1e-12;
  82.     /** Maximum number of iterations accepted in the implicit QL transformation. */
  83.     private static final byte MAX_ITER = 30;
  84.     /** Main diagonal of the tridiagonal matrix. */
  85.     private double[] main;
  86.     /** Secondary diagonal of the tridiagonal matrix. */
  87.     private double[] secondary;
  88.     /**
  89.      * Transformer to tridiagonal (may be null if matrix is already
  90.      * tridiagonal).
  91.      */
  92.     private TriDiagonalTransformer transformer;
  93.     /** Real part of the realEigenvalues. */
  94.     private double[] realEigenvalues;
  95.     /** Imaginary part of the realEigenvalues. */
  96.     private double[] imagEigenvalues;
  97.     /** Eigenvectors. */
  98.     private ArrayRealVector[] eigenvectors;
  99.     /** 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;

  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.     }

  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.     }

  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() {

  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.     }

  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() {

  190.         if (cachedD == null) {
  191.             // cache the matrix for subsequent calls
  192.             cachedD = MatrixUtils.createRealMatrixWithDiagonal(realEigenvalues);

  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.     }

  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() {

  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.         }

  225.         // return the cached matrix
  226.         return cachedVt;
  227.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.         }

  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();

  358.         return v.multiply(sqrtEigen).multiply(vT);
  359.     }

  361.     /**
  362.      * Gets a solver for finding the A &times; 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.     }

  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;

  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.         }

  403.         /**
  404.          * Solves the linear equation A &times; X = B for symmetric matrices A.
  405.          * <p>
  406.          * This method only finds exact linear solutions, i.e. solutions for
  407.          * which ||A &times; X - B|| is exactly 0.
  408.          * </p>
  409.          *
  410.          * @param b Right-hand side of the equation A &times; X = B.
  411.          * @return a Vector X that minimizes the two norm of A &times; 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.             }

  422.             final int m = realEigenvalues.length;
  423.             if (b.getDimension() != m) {
  424.                 throw new DimensionMismatchException(b.getDimension(), m);
  425.             }

  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.             }

  437.             return new ArrayRealVector(bp, false);
  438.         }

  440.         /** {@inheritDoc} */
  441.         @Override
  442.         public RealMatrix solve(RealMatrix b) {

  444.             if (!isNonSingular()) {
  445.                 throw new SingularMatrixException();
  446.             }

  448.             final int m = realEigenvalues.length;
  449.             if (b.getRowDimension() != m) {
  450.                 throw new DimensionMismatchException(b.getRowDimension(), m);
  451.             }

  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.             }

  475.             return new Array2DRowRealMatrix(bp, false);
  476.         }

  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.         }

  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.         }

  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.             }

  527.             final int m = realEigenvalues.length;
  528.             final double[][] invData = new double[m][m];

  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.     }

  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.     }

  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;

  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.         }

  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.         }

  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.         }

  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.     }

  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();

  724.         realEigenvalues = new double[matT.length];
  725.         imagEigenvalues = new double[matT.length];

  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.     }

  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.     }

  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();

  770.         final int n = matrixT.length;

  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.         }

  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.         }

  785.         // Backsubstitute to find vectors of upper triangular form

  787.         double r = 0.0;
  788.         double s = 0.0;
  789.         double z = 0.0;

  791.         for (int idx = n - 1; idx >= 0; idx--) {
  792.             double p = realEigenvalues[idx];
  793.             double q = imagEigenvalues[idx];

  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.                         }

  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;

  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.                 }

  855.                 matrixT[idx][idx - 1] = 0.0;
  856.                 matrixT[idx][idx]     = 1.0;

  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;

  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();

  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.                         }

  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.         }

  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.         }

  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. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT