ComplexEigenDecomposition.java

  1. /*
  2.  * Licensed to the Hipparchus project 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 Hipparchus project 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.  */
  17. package org.hipparchus.linear;

  19. import java.lang.reflect.Array;

  21. import org.hipparchus.complex.Complex;
  22. import org.hipparchus.complex.ComplexField;
  23. import org.hipparchus.exception.LocalizedCoreFormats;
  24. import org.hipparchus.exception.MathIllegalArgumentException;
  25. import org.hipparchus.exception.MathRuntimeException;
  26. import org.hipparchus.util.FastMath;
  27. import org.hipparchus.util.Precision;

  29. /**
  30.  * Given a matrix A, it computes a complex eigen decomposition AV = VD.
  31.  *
  32.  * <p>
  33.  * Complex Eigen Decomposition differs from the {@link EigenDecompositionSymmetric} since it
  34.  * computes the eigen vectors as complex eigen vectors (if applicable).
  35.  * </p>
  36.  *
  37.  * <p>
  38.  * Beware that in the complex case, you do not always have \(V \times V^{T} = I\) or even a
  39.  * diagonal matrix, even if the eigenvectors that form the columns of the V
  40.  * matrix are independent. On example is the square matrix
  41.  * \[
  42.  * A = \left(\begin{matrix}
  43.  * 3 &amp; -2\\
  44.  * 4 &amp; -1
  45.  * \end{matrix}\right)
  46.  * \]
  47.  * which has two conjugate eigenvalues \(\lambda_1=1+2i\) and \(\lambda_2=1-2i\)
  48.  * with associated eigenvectors \(v_1^T = (1, 1-i)\) and \(v_2^T = (1, 1+i)\).
  49.  * \[
  50.  * V\timesV^T = \left(\begin{matrix}
  51.  * 2 &amp; 2\\
  52.  * 2 &amp; 0
  53.  * \end{matrix}\right)
  54.  * \]
  55.  * which is not the identity matrix. Therefore, despite \(A \times V = V \times D\),
  56.  * \(A \ne V \times D \time V^T\), which would hold for real eigendecomposition.
  57.  * </p>
  58.  * <p>
  59.  * Also note that for consistency with Wolfram langage
  60.  * <a href="https://reference.wolfram.com/language/ref/Eigenvectors.html">eigenvectors</a>,
  61.  * we add zero vectors when the geometric multiplicity of the eigenvalue is smaller than
  62.  * its algebraic multiplicity (hence the regular eigenvector matrix should be non-square).
  63.  * With these additional null vectors, the eigenvectors matrix becomes square. This happens
  64.  * for example with the square matrix
  65.  * \[
  66.  * A = \left(\begin{matrix}
  67.  *  1 &amp; 0 &amp; 0\\
  68.  * -2 &amp; 1 &amp; 0\\
  69.  *  0 &amp; 0 &amp; 1
  70.  * \end{matrix}\right)
  71.  * \]
  72.  * Its characteristic polynomial is \((1-\lambda)^3\), hence is has one eigen value \(\lambda=1\)
  73.  * with algebraic multiplicity 3. However, this eigenvalue leads to only two eigenvectors
  74.  * \(v_1=(0, 1, 0)\) and \(v_2=(0, 0, 1)\), hence its geometric multiplicity is only 2, not 3.
  75.  * So we add a third zero vector \(v_3=(0, 0, 0)\), in the same way Wolfram language does.
  76.  * </p>
  77.  *
  78.  * Compute complex eigen values from the Schur transform. Compute complex eigen
  79.  * vectors based on eigen values and the inverse iteration method.
  80.  *
  81.  * see: <a href="https://en.wikipedia.org/wiki/Inverse_iteration">Inverse iteration</a>
  82.  * <a href="https://en.wikiversity.org/wiki/Shifted_inverse_iteration">Shifted inverse iteration</a>
  83.  * <a href="https://www.robots.ox.ac.uk/~sjrob/Teaching/EngComp/ecl4.pdf">Computation of matrix eigenvalues and eigenvectors</a>
  84.  *
  85.  */
  86. public class ComplexEigenDecomposition {

  88.     /** Default threshold below which eigenvectors are considered equal. */
  89.     public static final double DEFAULT_EIGENVECTORS_EQUALITY = 1.0e-5;
  90.     /** Default value to use for internal epsilon. */
  91.     public static final double DEFAULT_EPSILON = 1e-12;
  92.     /** Internally used epsilon criteria for final AV=VD check. */
  93.     public static final double DEFAULT_EPSILON_AV_VD_CHECK = 1e-6;
  94.     /** Maximum number of inverse iterations. */
  95.     private static final int MAX_ITER = 10;
  96.     /** complex eigenvalues. */
  97.     private Complex[] eigenvalues;
  98.     /** Eigenvectors. */
  99.     private FieldVector<Complex>[] eigenvectors;
  100.     /** Cached value of V. */
  101.     private FieldMatrix<Complex> V;
  102.     /** Cached value of D. */
  103.     private FieldMatrix<Complex> D;
  104.     /** Internally used threshold below which eigenvectors are considered equal. */
  105.     private final double eigenVectorsEquality;
  106.     /** Internally used epsilon criteria. */
  107.     private final double epsilon;
  108.     /** Internally used epsilon criteria for final AV=VD check. */
  109.     private final double epsilonAVVDCheck;

  111.     /**
  112.      * Constructor for decomposition.
  113.      * <p>
  114.      * This constructor uses the default values {@link #DEFAULT_EIGENVECTORS_EQUALITY},
  115.      * {@link #DEFAULT_EPSILON} and {@link #DEFAULT_EPSILON_AV_VD_CHECK}
  116.      * </p>
  117.      * @param matrix
  118.      *            real matrix.
  119.      */
  120.     public ComplexEigenDecomposition(final RealMatrix matrix) {
  121.         this(matrix, DEFAULT_EIGENVECTORS_EQUALITY,
  122.              DEFAULT_EPSILON, DEFAULT_EPSILON_AV_VD_CHECK);
  123.     }

  125.     /**
  126.      * Constructor for decomposition.
  127.      * <p>
  128.      * The {@code eigenVectorsEquality} threshold is used to ensure the L∞-normalized
  129.      * eigenvectors found using inverse iteration are different from each other.
  130.      * if \(min(|e_i-e_j|,|e_i+e_j|)\) is smaller than this threshold, the algorithm
  131.      * considers it has found again an already known vector, so it drops it and attempts
  132.      * a new inverse iteration with a different start vector. This value should be
  133.      * much larger than {@code epsilon} which is used for convergence
  134.      * </p>
  135.      * @param matrix real matrix.
  136.      * @param eigenVectorsEquality threshold below which eigenvectors are considered equal
  137.      * @param epsilon Epsilon used for internal tests (e.g. is singular, eigenvalue ratio, etc.)
  138.      * @param epsilonAVVDCheck Epsilon criteria for final AV=VD check
  139.      * @since 1.8
  140.      */
  141.     public ComplexEigenDecomposition(final RealMatrix matrix, final double eigenVectorsEquality,
  142.                                      final double epsilon, final double epsilonAVVDCheck) {

  144.         if (!matrix.isSquare()) {
  145.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SQUARE_MATRIX,
  146.                                                    matrix.getRowDimension(), matrix.getColumnDimension());
  147.         }
  148.         this.eigenVectorsEquality = eigenVectorsEquality;
  149.         this.epsilon              = epsilon;
  150.         this.epsilonAVVDCheck     = epsilonAVVDCheck;

  152.         // computing the eigen values
  153.         findEigenValues(matrix);
  154.         // computing the eigen vectors
  155.         findEigenVectors(convertToFieldComplex(matrix));

  157.         // V
  158.         final int m = eigenvectors.length;
  159.         V = MatrixUtils.createFieldMatrix(ComplexField.getInstance(), m, m);
  160.         for (int k = 0; k < m; ++k) {
  161.             V.setColumnVector(k, eigenvectors[k]);
  162.         }

  164.         // D
  165.         D = MatrixUtils.createFieldDiagonalMatrix(eigenvalues);

  167.         checkDefinition(matrix);
  168.     }

  170.     /**
  171.      * Getter of the eigen values.
  172.      *
  173.      * @return eigen values.
  174.      */
  175.     public Complex[] getEigenvalues() {
  176.         return eigenvalues.clone();
  177.     }

  179.     /**
  180.      * Getter of the eigen vectors.
  181.      *
  182.      * @param i
  183.      *            which eigen vector.
  184.      * @return eigen vector.
  185.      */
  186.     public FieldVector<Complex> getEigenvector(final int i) {
  187.         return eigenvectors[i].copy();
  188.     }

  190.     /** Reset eigenvalues and eigen vectors from matrices.
  191.      * <p>
  192.      * This method is intended to be called by sub-classes (mainly {@link OrderedComplexEigenDecomposition})
  193.      * that reorder the matrices elements. It rebuild the eigenvalues and eigen vectors arrays
  194.      * from the D and V matrices.
  195.      * </p>
  196.      * @since 2.1
  197.      */
  198.     protected void matricesToEigenArrays() {
  199.         for (int i = 0; i < eigenvalues.length; ++i) {
  200.             eigenvalues[i] = D.getEntry(i, i);
  201.         }
  202.         for (int i = 0; i < eigenvectors.length; ++i) {
  203.             for (int j = 0; j < eigenvectors[i].getDimension(); ++j) {
  204.                 eigenvectors[i].setEntry(j, V.getEntry(j, i));
  205.             }
  206.         }
  207.     }

  209.     /**
  210.      * Confirm if there are complex eigen values.
  211.      *
  212.      * @return true if there are complex eigen values.
  213.      */
  214.     public boolean hasComplexEigenvalues() {
  215.         for (Complex eigenvalue : eigenvalues) {
  216.             if (!Precision.equals(eigenvalue.getImaginary(), 0.0, epsilon)) {
  217.                 return true;
  218.             }
  219.         }
  220.         return false;
  221.     }

  223.     /**
  224.      * Computes the determinant.
  225.      *
  226.      * @return the determinant.
  227.      */
  228.     public double getDeterminant() {
  229.         Complex determinant = new Complex(1, 0);
  230.         for (Complex lambda : eigenvalues) {
  231.             determinant = determinant.multiply(lambda);
  232.         }
  233.         return determinant.getReal();
  234.     }

  236.     /**
  237.      * Getter V.
  238.      *
  239.      * @return V.
  240.      */
  241.     public FieldMatrix<Complex> getV() {
  242.         return V;
  243.     }

  245.     /**
  246.      * Getter D.
  247.      *
  248.      * @return D.
  249.      */
  250.     public FieldMatrix<Complex> getD() {
  251.         return D;
  252.     }

  254.     /**
  255.      * Getter VT.
  256.      *
  257.      * @return VT.
  258.      */
  259.     public FieldMatrix<Complex> getVT() {
  260.         return V.transpose();
  261.     }

  263.     /**
  264.      * Compute eigen values using the Schur transform.
  265.      *
  266.      * @param matrix
  267.      *            real matrix to compute eigen values.
  268.      */
  269.     protected void findEigenValues(final RealMatrix matrix) {
  270.         final SchurTransformer schurTransform = new SchurTransformer(matrix);
  271.         final double[][] matT = schurTransform.getT().getData();

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

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

  288.     }

  290.     /**
  291.      * Compute the eigen vectors using the inverse power method.
  292.      *
  293.      * @param matrix
  294.      *            real matrix to compute eigen vectors.
  295.      */
  296.     @SuppressWarnings("unchecked")
  297.     protected void findEigenVectors(final FieldMatrix<Complex> matrix) {
  298.         // number of eigen values/vectors
  299.         int n = eigenvalues.length;

  301.         // eigen vectors
  302.         eigenvectors = (FieldVector<Complex>[]) Array.newInstance(FieldVector.class, n);

  304.         // computing eigen vector based on eigen values and inverse iteration
  305.         for (int i = 0; i < eigenvalues.length; i++) {

  307.             // shifted non-singular matrix matrix A-(λ+ε)I that is close to the singular matrix A-λI
  308.             Complex mu = eigenvalues[i].add(epsilon);
  309.             final FieldMatrix<Complex> shifted = matrix.copy();
  310.             for (int k = 0; k < matrix.getColumnDimension(); ++k) {
  311.                 shifted.setEntry(k, k, shifted.getEntry(k, k).subtract(mu));
  312.             }

  314.             // solver for linear system (A - (λ+ε)I) Bₖ₊₁ = Bₖ
  315.             FieldDecompositionSolver<Complex> solver = new FieldQRDecomposition<>(shifted).getSolver();

  317.             // loop over possible start vectors
  318.             for (int p = 0; eigenvectors[i] == null && p < matrix.getColumnDimension(); ++p) {

  320.                 // find a vector to start iterations
  321.                 FieldVector<Complex> b = findStart(p);

  323.                 if (getNorm(b).norm() > Precision.SAFE_MIN) {
  324.                     // start vector is a good candidate for inverse iteration

  326.                     // perform inverse iteration
  327.                     double delta = Double.POSITIVE_INFINITY;
  328.                     for (int k = 0; delta > epsilon && k < MAX_ITER; k++) {

  330.                         // solve (A - (λ+ε)) Bₖ₊₁ = Bₖ
  331.                         final FieldVector<Complex> bNext = solver.solve(b);

  333.                         // normalize according to L∞ norm
  334.                         normalize(bNext);

  336.                         // compute convergence criterion, comparing Bₖ and both ±Bₖ₊₁
  337.                         // as iterations sometimes flip between two opposite vectors
  338.                         delta = separation(b, bNext);

  340.                         // prepare next iteration
  341.                         b = bNext;

  343.                     }

  345.                     // check we have not found again an already known vector
  346.                     for (int j = 0; b != null && j < i; ++j) {
  347.                         if (separation(eigenvectors[j], b) <= eigenVectorsEquality) {
  348.                             // the selected start vector leads us to found a known vector again,
  349.                             // we must try another start
  350.                             b = null;
  351.                         }
  352.                     }
  353.                     eigenvectors[i] = b;

  355.                 }
  356.             }

  358.             if (eigenvectors[i] == null) {
  359.                 // for consistency with Wolfram langage
  360.                 // https://reference.wolfram.com/language/ref/Eigenvectors.html
  361.                 // we add zero vectors when the geometric multiplicity of the eigenvalue
  362.                 // is smaller than its algebraic multiplicity (hence the regular eigenvector
  363.                 // matrix should be non-square). With these additional null vectors, the
  364.                 // eigenvectors matrix becomes square
  365.                 eigenvectors[i] = MatrixUtils.createFieldVector(ComplexField.getInstance(), n);
  366.             }

  368.         }
  369.     }

  371.     /** Find a start vector orthogonal to all already found normalized eigenvectors.
  372.      * @param index index of the vector
  373.      * @return start vector
  374.      */
  375.     private FieldVector<Complex> findStart(final int index) {

  377.         // create vector
  378.         final FieldVector<Complex> start =
  379.                         MatrixUtils.createFieldVector(ComplexField.getInstance(),
  380.                                                       eigenvalues.length);

  382.         // initialize with a canonical vector
  383.         start.setEntry(index, Complex.ONE);

  385.         return start;

  387.     }

  389.     /**
  390.      * Compute the L∞ norm of the a given vector.
  391.      *
  392.      * @param vector
  393.      *            vector.
  394.      * @return L∞ norm.
  395.      */
  396.     private Complex getNorm(FieldVector<Complex> vector) {
  397.         double  normR = 0;
  398.         Complex normC = Complex.ZERO;
  399.         for (int i = 0; i < vector.getDimension(); i++) {
  400.             final Complex ci = vector.getEntry(i);
  401.             final double  ni = FastMath.hypot(ci.getReal(), ci.getImaginary());
  402.             if (ni > normR) {
  403.                 normR = ni;
  404.                 normC = ci;
  405.             }
  406.         }
  407.         return normC;
  408.     }

  410.     /** Normalize a vector with respect to L∞ norm.
  411.      * @param v vector to normalized
  412.      */
  413.     private void normalize(final FieldVector<Complex> v) {
  414.         final Complex invNorm = getNorm(v).reciprocal();
  415.         for (int j = 0; j < v.getDimension(); ++j) {
  416.             v.setEntry(j, v.getEntry(j).multiply(invNorm));
  417.         }
  418.     }

  420.     /** Compute the separation between two normalized vectors (which may be in opposite directions).
  421.      * @param v1 first normalized vector
  422.      * @param v2 second normalized vector
  423.      * @return min (|v1 - v2|, |v1+v2|)
  424.      */
  425.     private double separation(final FieldVector<Complex> v1, final FieldVector<Complex> v2) {
  426.         double deltaPlus  = 0;
  427.         double deltaMinus = 0;
  428.         for (int j = 0; j < v1.getDimension(); ++j) {
  429.             final Complex bCurrj = v1.getEntry(j);
  430.             final Complex bNextj = v2.getEntry(j);
  431.             deltaPlus  = FastMath.max(deltaPlus,
  432.                                       FastMath.hypot(bNextj.getReal()      + bCurrj.getReal(),
  433.                                                      bNextj.getImaginary() + bCurrj.getImaginary()));
  434.             deltaMinus = FastMath.max(deltaMinus,
  435.                                       FastMath.hypot(bNextj.getReal()      - bCurrj.getReal(),
  436.                                                      bNextj.getImaginary() - bCurrj.getImaginary()));
  437.         }
  438.         return FastMath.min(deltaPlus, deltaMinus);
  439.     }

  441.     /**
  442.      * Check definition of the decomposition in runtime.
  443.      *
  444.      * @param matrix
  445.      *            matrix to be decomposed.
  446.      */
  447.     protected void checkDefinition(final RealMatrix matrix) {
  448.         FieldMatrix<Complex> matrixC = convertToFieldComplex(matrix);

  450.         // checking definition of the decomposition
  451.         // testing A*V = V*D
  452.         FieldMatrix<Complex> AV = matrixC.multiply(getV());
  453.         FieldMatrix<Complex> VD = getV().multiply(getD());
  454.         if (!equalsWithPrecision(AV, VD, epsilonAVVDCheck)) {
  455.             throw new MathRuntimeException(LocalizedCoreFormats.FAILED_DECOMPOSITION,
  456.                                            matrix.getRowDimension(), matrix.getColumnDimension());

  458.         }

  460.     }

  462.     /**
  463.      * Helper method that checks with two matrix is equals taking into account a
  464.      * given precision.
  465.      *
  466.      * @param matrix1 first matrix to compare
  467.      * @param matrix2 second matrix to compare
  468.      * @param tolerance tolerance on matrices entries
  469.      * @return true is matrices entries are equal within tolerance,
  470.      * false otherwise
  471.      */
  472.     private boolean equalsWithPrecision(final FieldMatrix<Complex> matrix1,
  473.                                         final FieldMatrix<Complex> matrix2, final double tolerance) {
  474.         boolean toRet = true;
  475.         for (int i = 0; i < matrix1.getRowDimension(); i++) {
  476.             for (int j = 0; j < matrix1.getColumnDimension(); j++) {
  477.                 Complex c1 = matrix1.getEntry(i, j);
  478.                 Complex c2 = matrix2.getEntry(i, j);
  479.                 if (c1.add(c2.negate()).norm() > tolerance) {
  480.                     toRet = false;
  481.                     break;
  482.                 }
  483.             }
  484.         }
  485.         return toRet;
  486.     }

  488.     /**
  489.      * It converts a real matrix into a complex field matrix.
  490.      *
  491.      * @param matrix
  492.      *            real matrix.
  493.      * @return complex matrix.
  494.      */
  495.     private FieldMatrix<Complex> convertToFieldComplex(RealMatrix matrix) {
  496.         final FieldMatrix<Complex> toRet =
  497.                         MatrixUtils.createFieldIdentityMatrix(ComplexField.getInstance(),
  498.                                                               matrix.getRowDimension());
  499.         for (int i = 0; i < toRet.getRowDimension(); i++) {
  500.             for (int j = 0; j < toRet.getColumnDimension(); j++) {
  501.                 toRet.setEntry(i, j, new Complex(matrix.getEntry(i, j)));
  502.             }
  503.         }
  504.         return toRet;
  505.     }
  506. }
Created with JaCoCo 0.8.12.202403310830