EigenDecompositionNonSymmetric.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.  *      https://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. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */

  23. package org.hipparchus.linear;

  25. import java.util.ArrayList;
  26. import java.util.List;

  28. import org.hipparchus.complex.Complex;
  29. import org.hipparchus.complex.ComplexField;
  30. import org.hipparchus.exception.LocalizedCoreFormats;
  31. import org.hipparchus.exception.MathIllegalStateException;
  32. import org.hipparchus.exception.MathRuntimeException;
  33. import org.hipparchus.util.FastMath;
  34. import org.hipparchus.util.Precision;

  36. /**
  37.  * Calculates the eigen decomposition of a non-symmetric real matrix.
  38.  * <p>
  39.  * The eigen decomposition of matrix A is a set of two matrices:
  40.  * \(V\) and \(D\) such that \(A V = V D\) where $\(A\),
  41.  * \(V\) and \(D\) are all \(m \times m\) matrices.
  42.  * <p>
  43.  * This class is similar in spirit to the {@code EigenvalueDecomposition}
  44.  * class from the <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a>
  45.  * library, with the following changes:
  46.  * <ul>
  47.  *   <li>a {@link #getVInv() getVInv} method has been added,</li>
  48.  *   <li>z {@link #getEigenvalue(int) getEigenvalue} method to pick up a
  49.  *       single eigenvalue has been added,</li>
  50.  *   <li>a {@link #getEigenvector(int) getEigenvector} method to pick up a
  51.  *       single eigenvector has been added,</li>
  52.  *   <li>a {@link #getDeterminant() getDeterminant} method has been added.</li>
  53.  * </ul>
  54.  * <p>
  55.  * This class supports non-symmetric matrices, which have complex eigenvalues.
  56.  * Support for symmetric matrices is provided by {@link EigenDecompositionSymmetric}.
  57.  * </p>
  58.  * <p>
  59.  * As \(A\) is not symmetric, then the eigenvalue matrix \(D\) is block diagonal with
  60.  * the real eigenvalues in 1-by-1 blocks and any complex eigenvalues, \(\lambda \pm i \mu\),
  61.  * in 2-by-2 blocks:
  62.  * </p>
  63.  * <p>
  64.  * \[
  65.  *   \begin{bmatrix}
  66.  *    \lambda &amp; \mu\\
  67.  *    -\mu    &amp; \lambda
  68.  *   \end{bmatrix}
  69.  * \]
  70.  * </p>
  71.  * <p>
  72.  * The columns of \(V\) represent the eigenvectors in the sense that \(A V = V D\),
  73.  * i.e. {@code A.multiply(V)} equals {@code V.multiply(D)}.
  74.  * The matrix \(V\) may be badly conditioned, or even singular, so the validity of the
  75.  * equation \(A = V D V^{-1}\) depends upon the condition of \(V\).
  76.  * </p>
  77.  * <p>
  78.  * This implementation is based on the paper by A. Drubrulle, R.S. Martin and
  79.  * J.H. Wilkinson "The Implicit QL Algorithm" in Wilksinson and Reinsch (1971)
  80.  * Handbook for automatic computation, vol. 2, Linear algebra, Springer-Verlag,
  81.  * New-York.
  82.  *
  83.  * @see <a href="http://mathworld.wolfram.com/EigenDecomposition.html">MathWorld</a>
  84.  * @see <a href="http://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix">Wikipedia</a>
  85.  * @since 3.0
  86.  */
  87. public class EigenDecompositionNonSymmetric {
  88.     /** Default epsilon value to use for internal epsilon **/
  89.     public static final double DEFAULT_EPSILON = 1e-12;
  90.     /** Internally used epsilon criteria. */
  91.     private final double epsilon;
  92.     /** eigenvalues. */
  93.     private Complex[] eigenvalues;
  94.     /** Eigenvectors. */
  95.     private List<FieldVector<Complex>> eigenvectors;
  96.     /** Cached value of \(V\). */
  97.     private RealMatrix cachedV;
  98.     /** Cached value of \(D\). */
  99.     private RealMatrix cachedD;
  100.     /** Cached value of \(V^{-1}\). */
  101.     private RealMatrix cachedVInv;

  103.     /** Calculates the eigen decomposition of the given real matrix.
  104.      * @param matrix Matrix to decompose.
  105.      * @throws MathIllegalStateException if the algorithm fails to converge.
  106.      * @throws MathRuntimeException if the decomposition of a general matrix
  107.      * results in a matrix with zero norm
  108.      */
  109.     public EigenDecompositionNonSymmetric(final RealMatrix matrix) {
  110.         this(matrix, DEFAULT_EPSILON);
  111.     }

  113.     /** Calculates the eigen decomposition of the given real matrix.
  114.      * @param matrix Matrix to decompose.
  115.      * @param epsilon Epsilon used for internal tests (e.g. is singular, eigenvalue ratio, etc.)
  116.      * @throws MathIllegalStateException if the algorithm fails to converge.
  117.      * @throws MathRuntimeException if the decomposition of a general matrix
  118.      * results in a matrix with zero norm
  119.      */
  120.     public EigenDecompositionNonSymmetric(final RealMatrix matrix, double epsilon)
  121.         throws MathRuntimeException {
  122.         this.epsilon = epsilon;
  123.         findEigenVectorsFromSchur(transformToSchur(matrix));
  124.     }

  126.     /**
  127.      * Gets the matrix V of the decomposition.
  128.      * V is a matrix whose columns hold either the real or the
  129.      * imaginary part of eigenvectors.
  130.      *
  131.      * @return the V matrix.
  132.      */
  133.     public RealMatrix getV() {

  135.         if (cachedV == null) {
  136.             final int m = eigenvectors.size();
  137.             cachedV = MatrixUtils.createRealMatrix(m, m);
  138.             for (int k = 0; k < m; ++k) {
  139.                 final FieldVector<Complex> ek = eigenvectors.get(k);
  140.                 if (eigenvalues[k].getImaginaryPart() >= 0) {
  141.                     // either it is a real eigenvalue, or it is the first of two conjugate eigenvalues,
  142.                     // we pick up the real part of the eigenvector
  143.                     for (int l = 0; l < m; ++l) {
  144.                         cachedV.setEntry(l, k, ek.getEntry(l).getRealPart());
  145.                     }
  146.                 } else {
  147.                     // second of two conjugate eigenvalues,
  148.                     // we pick up the imaginary part of the eigenvector
  149.                     for (int l = 0; l < m; ++l) {
  150.                         cachedV.setEntry(l, k, -ek.getEntry(l).getImaginaryPart());
  151.                     }
  152.                 }
  153.             }
  154.         }

  156.         // return the cached matrix
  157.         return cachedV;

  159.     }

  161.     /**
  162.      * Gets the block diagonal matrix D of the decomposition.
  163.      * D is a block diagonal matrix.
  164.      * Real eigenvalues are on the diagonal while complex values are on
  165.      * 2x2 blocks { {real +imaginary}, {-imaginary, real} }.
  166.      *
  167.      * @return the D matrix.
  168.      */
  169.     public RealMatrix getD() {

  171.         if (cachedD == null) {
  172.             // cache the matrix for subsequent calls
  173.             cachedD = MatrixUtils.createRealMatrix(eigenvalues.length, eigenvalues.length);
  174.             for (int i = 0; i < eigenvalues.length; ++i) {
  175.                 cachedD.setEntry(i, i, eigenvalues[i].getRealPart());
  176.                 if (Precision.compareTo(eigenvalues[i].getImaginaryPart(), 0.0, epsilon) > 0) {
  177.                     cachedD.setEntry(i, i + 1, eigenvalues[i].getImaginaryPart());
  178.                 } else if (Precision.compareTo(eigenvalues[i].getImaginaryPart(), 0.0, epsilon) < 0) {
  179.                     cachedD.setEntry(i, i - 1, eigenvalues[i].getImaginaryPart());
  180.                 }
  181.             }
  182.         }
  183.         return cachedD;
  184.     }

  186.     /**
  187.      * Get's the value for epsilon which is used for internal tests (e.g. is singular, eigenvalue ratio, etc.)
  188.      *
  189.      * @return the epsilon value.
  190.      */
  191.     public double getEpsilon() {
  192.         return epsilon;
  193.     }

  195.     /**
  196.      * Gets the inverse of the matrix V of the decomposition.
  197.      *
  198.      * @return the inverse of the V matrix.
  199.      */
  200.     public RealMatrix getVInv() {

  202.         if (cachedVInv == null) {
  203.             cachedVInv = MatrixUtils.inverse(getV(), epsilon);
  204.         }

  206.         // return the cached matrix
  207.         return cachedVInv;
  208.     }

  210.     /**
  211.      * Gets a copy of the eigenvalues of the original matrix.
  212.      *
  213.      * @return a copy of the eigenvalues of the original matrix.
  214.      *
  215.      * @see #getD()
  216.      * @see #getEigenvalue(int)
  217.      */
  218.     public Complex[] getEigenvalues() {
  219.         return eigenvalues.clone();
  220.     }

  222.     /**
  223.      * Returns the i<sup>th</sup> eigenvalue of the original matrix.
  224.      *
  225.      * @param i index of the eigenvalue (counting from 0)
  226.      * @return i<sup>th</sup> eigenvalue of the original matrix.
  227.      *
  228.      * @see #getD()
  229.      * @see #getEigenvalues()
  230.      */
  231.     public Complex getEigenvalue(final int i) {
  232.         return eigenvalues[i];
  233.     }

  235.     /**
  236.      * Gets a copy of the i<sup>th</sup> eigenvector of the original matrix.
  237.      * <p>
  238.      * Note that if the the i<sup>th</sup> is complex this method will throw
  239.      * an exception.
  240.      * </p>
  241.      * @param i Index of the eigenvector (counting from 0).
  242.      * @return a copy of the i<sup>th</sup> eigenvector of the original matrix.
  243.      * @see #getD()
  244.      */
  245.     public FieldVector<Complex> getEigenvector(final int i) {
  246.         return eigenvectors.get(i).copy();
  247.     }

  249.     /**
  250.      * Computes the determinant of the matrix.
  251.      *
  252.      * @return the determinant of the matrix.
  253.      */
  254.     public Complex getDeterminant() {
  255.         Complex determinant = Complex.ONE;
  256.         for (Complex eigenvalue : eigenvalues) {
  257.             determinant = determinant.multiply(eigenvalue);
  258.         }
  259.         return determinant;
  260.     }

  262.     /**
  263.      * Transforms the matrix to Schur form and calculates the eigenvalues.
  264.      *
  265.      * @param matrix Matrix to transform.
  266.      * @return the {@link SchurTransformer Schur transform} for this matrix
  267.      */
  268.     private SchurTransformer transformToSchur(final RealMatrix matrix) {
  269.         final SchurTransformer schurTransform = new SchurTransformer(matrix, epsilon);
  270.         final double[][] matT = schurTransform.getT().getData();
  271.         final double norm = matrix.getNorm1();

  273.         eigenvalues = new Complex[matT.length];

  275.         int i = 0;
  276.         while (i < eigenvalues.length) {
  277.             if (i == (eigenvalues.length - 1) ||
  278.                 Precision.equals(matT[i + 1][i], 0.0, norm * epsilon)) {
  279.                 eigenvalues[i] = new Complex(matT[i][i]);
  280.                 i++;
  281.             } else {
  282.                 final double x   = matT[i + 1][i + 1];
  283.                 final double p   = 0.5 * (matT[i][i] - x);
  284.                 final double z   = FastMath.sqrt(FastMath.abs(p * p + matT[i + 1][i] * matT[i][i + 1]));
  285.                 eigenvalues[i++] = new Complex(x + p, +z);
  286.                 eigenvalues[i++] = new Complex(x + p, -z);
  287.             }
  288.         }
  289.         return schurTransform;
  290.     }

  292.     /**
  293.      * Performs a division of two complex numbers.
  294.      *
  295.      * @param xr real part of the first number
  296.      * @param xi imaginary part of the first number
  297.      * @param yr real part of the second number
  298.      * @param yi imaginary part of the second number
  299.      * @return result of the complex division
  300.      */
  301.     private Complex cdiv(final double xr, final double xi,
  302.                          final double yr, final double yi) {
  303.         return new Complex(xr, xi).divide(new Complex(yr, yi));
  304.     }

  306.     /**
  307.      * Find eigenvectors from a matrix transformed to Schur form.
  308.      *
  309.      * @param schur the schur transformation of the matrix
  310.      * @throws MathRuntimeException if the Schur form has a norm of zero
  311.      */
  312.     private void findEigenVectorsFromSchur(final SchurTransformer schur)
  313.         throws MathRuntimeException {
  314.         final double[][] matrixT = schur.getT().getData();
  315.         final double[][] matrixP = schur.getP().getData();

  317.         final int n = matrixT.length;

  319.         // compute matrix norm
  320.         double norm = 0.0;
  321.         for (int i = 0; i < n; i++) {
  322.            for (int j = FastMath.max(i - 1, 0); j < n; j++) {
  323.                norm += FastMath.abs(matrixT[i][j]);
  324.            }
  325.         }

  327.         // we can not handle a matrix with zero norm
  328.         if (norm == 0.0) {
  329.            throw new MathRuntimeException(LocalizedCoreFormats.ZERO_NORM);
  330.         }

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

  334.         double r = 0.0;
  335.         double s = 0.0;
  336.         double z = 0.0;

  338.         for (int idx = n - 1; idx >= 0; idx--) {
  339.             double p = eigenvalues[idx].getRealPart();
  340.             double q = eigenvalues[idx].getImaginaryPart();

  342.             if (Precision.equals(q, 0.0)) {
  343.                 // Real vector
  344.                 int l = idx;
  345.                 matrixT[idx][idx] = 1.0;
  346.                 for (int i = idx - 1; i >= 0; i--) {
  347.                     double w = matrixT[i][i] - p;
  348.                     r = 0.0;
  349.                     for (int j = l; j <= idx; j++) {
  350.                         r += matrixT[i][j] * matrixT[j][idx];
  351.                     }
  352.                     if (Precision.compareTo(eigenvalues[i].getImaginaryPart(), 0.0, epsilon) < 0) {
  353.                         z = w;
  354.                         s = r;
  355.                     } else {
  356.                         l = i;
  357.                         if (Precision.equals(eigenvalues[i].getImaginaryPart(), 0.0)) {
  358.                             if (w != 0.0) {
  359.                                 matrixT[i][idx] = -r / w;
  360.                             } else {
  361.                                 matrixT[i][idx] = -r / (Precision.EPSILON * norm);
  362.                             }
  363.                         } else {
  364.                             // Solve real equations
  365.                             double x = matrixT[i][i + 1];
  366.                             double y = matrixT[i + 1][i];
  367.                             q = (eigenvalues[i].getRealPart() - p) * (eigenvalues[i].getRealPart() - p) +
  368.                                 eigenvalues[i].getImaginaryPart() * eigenvalues[i].getImaginaryPart();
  369.                             double t = (x * s - z * r) / q;
  370.                             matrixT[i][idx] = t;
  371.                             if (FastMath.abs(x) > FastMath.abs(z)) {
  372.                                 matrixT[i + 1][idx] = (-r - w * t) / x;
  373.                             } else {
  374.                                 matrixT[i + 1][idx] = (-s - y * t) / z;
  375.                             }
  376.                         }

  378.                         // Overflow control
  379.                         double t = FastMath.abs(matrixT[i][idx]);
  380.                         if ((Precision.EPSILON * t) * t > 1) {
  381.                             for (int j = i; j <= idx; j++) {
  382.                                 matrixT[j][idx] /= t;
  383.                             }
  384.                         }
  385.                     }
  386.                 }
  387.             } else if (q < 0.0) {
  388.                 // Complex vector
  389.                 int l = idx - 1;

  391.                 // Last vector component imaginary so matrix is triangular
  392.                 if (FastMath.abs(matrixT[idx][idx - 1]) > FastMath.abs(matrixT[idx - 1][idx])) {
  393.                     matrixT[idx - 1][idx - 1] = q / matrixT[idx][idx - 1];
  394.                     matrixT[idx - 1][idx]     = -(matrixT[idx][idx] - p) / matrixT[idx][idx - 1];
  395.                 } else {
  396.                     final Complex result = cdiv(0.0, -matrixT[idx - 1][idx],
  397.                                                 matrixT[idx - 1][idx - 1] - p, q);
  398.                     matrixT[idx - 1][idx - 1] = result.getReal();
  399.                     matrixT[idx - 1][idx]     = result.getImaginary();
  400.                 }

  402.                 matrixT[idx][idx - 1] = 0.0;
  403.                 matrixT[idx][idx]     = 1.0;

  405.                 for (int i = idx - 2; i >= 0; i--) {
  406.                     double ra = 0.0;
  407.                     double sa = 0.0;
  408.                     for (int j = l; j <= idx; j++) {
  409.                         ra += matrixT[i][j] * matrixT[j][idx - 1];
  410.                         sa += matrixT[i][j] * matrixT[j][idx];
  411.                     }
  412.                     double w = matrixT[i][i] - p;

  414.                     if (Precision.compareTo(eigenvalues[i].getImaginaryPart(), 0.0, epsilon) < 0) {
  415.                         z = w;
  416.                         r = ra;
  417.                         s = sa;
  418.                     } else {
  419.                         l = i;
  420.                         if (Precision.equals(eigenvalues[i].getImaginaryPart(), 0.0)) {
  421.                             final Complex c = cdiv(-ra, -sa, w, q);
  422.                             matrixT[i][idx - 1] = c.getReal();
  423.                             matrixT[i][idx] = c.getImaginary();
  424.                         } else {
  425.                             // Solve complex equations
  426.                             double x = matrixT[i][i + 1];
  427.                             double y = matrixT[i + 1][i];
  428.                             double vr = (eigenvalues[i].getRealPart() - p) * (eigenvalues[i].getRealPart() - p) +
  429.                                         eigenvalues[i].getImaginaryPart() * eigenvalues[i].getImaginaryPart() - q * q;
  430.                             final double vi = (eigenvalues[i].getRealPart() - p) * 2.0 * q;
  431.                             if (Precision.equals(vr, 0.0) && Precision.equals(vi, 0.0)) {
  432.                                 vr = Precision.EPSILON * norm *
  433.                                      (FastMath.abs(w) + FastMath.abs(q) + FastMath.abs(x) +
  434.                                       FastMath.abs(y) + FastMath.abs(z));
  435.                             }
  436.                             final Complex c     = cdiv(x * r - z * ra + q * sa,
  437.                                                        x * s - z * sa - q * ra, vr, vi);
  438.                             matrixT[i][idx - 1] = c.getReal();
  439.                             matrixT[i][idx]     = c.getImaginary();

  441.                             if (FastMath.abs(x) > (FastMath.abs(z) + FastMath.abs(q))) {
  442.                                 matrixT[i + 1][idx - 1] = (-ra - w * matrixT[i][idx - 1] +
  443.                                                            q * matrixT[i][idx]) / x;
  444.                                 matrixT[i + 1][idx]     = (-sa - w * matrixT[i][idx] -
  445.                                                            q * matrixT[i][idx - 1]) / x;
  446.                             } else {
  447.                                 final Complex c2        = cdiv(-r - y * matrixT[i][idx - 1],
  448.                                                                -s - y * matrixT[i][idx], z, q);
  449.                                 matrixT[i + 1][idx - 1] = c2.getReal();
  450.                                 matrixT[i + 1][idx]     = c2.getImaginary();
  451.                             }
  452.                         }

  454.                         // Overflow control
  455.                         double t = FastMath.max(FastMath.abs(matrixT[i][idx - 1]),
  456.                                                 FastMath.abs(matrixT[i][idx]));
  457.                         if ((Precision.EPSILON * t) * t > 1) {
  458.                             for (int j = i; j <= idx; j++) {
  459.                                 matrixT[j][idx - 1] /= t;
  460.                                 matrixT[j][idx] /= t;
  461.                             }
  462.                         }
  463.                     }
  464.                 }
  465.             }
  466.         }

  468.         // Back transformation to get eigenvectors of original matrix
  469.         for (int j = n - 1; j >= 0; j--) {
  470.             for (int i = 0; i <= n - 1; i++) {
  471.                 z = 0.0;
  472.                 for (int k = 0; k <= FastMath.min(j, n - 1); k++) {
  473.                     z += matrixP[i][k] * matrixT[k][j];
  474.                 }
  475.                 matrixP[i][j] = z;
  476.             }
  477.         }

  479.         eigenvectors = new ArrayList<>(n);
  480.         for (int i = 0; i < n; i++) {
  481.             FieldVector<Complex> ei = new ArrayFieldVector<>(ComplexField.getInstance(), n);
  482.             for (int j = 0; j < n; j++) {
  483.                 if (Precision.compareTo(eigenvalues[i].getImaginaryPart(), 0.0, epsilon) > 0) {
  484.                     ei.setEntry(j, new Complex(matrixP[j][i], +matrixP[j][i + 1]));
  485.                 } else if (Precision.compareTo(eigenvalues[i].getImaginaryPart(), 0.0, epsilon) < 0) {
  486.                     ei.setEntry(j, new Complex(matrixP[j][i - 1], -matrixP[j][i]));
  487.                 } else {
  488.                     ei.setEntry(j, new Complex(matrixP[j][i]));
  489.                 }
  490.             }
  491.             eigenvectors.add(ei);
  492.         }
  493.     }
  494. }
Created with JaCoCo 0.8.12.202403310830