ComplexEigenDecomposition.java
/*
* Licensed to the Hipparchus project under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The Hipparchus project licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* https://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.hipparchus.linear;
import java.lang.reflect.Array;
import org.hipparchus.complex.Complex;
import org.hipparchus.complex.ComplexField;
import org.hipparchus.exception.LocalizedCoreFormats;
import org.hipparchus.exception.MathIllegalArgumentException;
import org.hipparchus.exception.MathRuntimeException;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.Precision;
/**
* Given a matrix A, it computes a complex eigen decomposition AV = VD.
*
* <p>
* Complex Eigen Decomposition differs from the {@link EigenDecompositionSymmetric} since it
* computes the eigen vectors as complex eigen vectors (if applicable).
* </p>
*
* <p>
* Beware that in the complex case, you do not always have \(V \times V^{T} = I\) or even a
* diagonal matrix, even if the eigenvectors that form the columns of the V
* matrix are independent. On example is the square matrix
* \[
* A = \left(\begin{matrix}
* 3 & -2\\
* 4 & -1
* \end{matrix}\right)
* \]
* which has two conjugate eigenvalues \(\lambda_1=1+2i\) and \(\lambda_2=1-2i\)
* with associated eigenvectors \(v_1^T = (1, 1-i)\) and \(v_2^T = (1, 1+i)\).
* \[
* V\timesV^T = \left(\begin{matrix}
* 2 & 2\\
* 2 & 0
* \end{matrix}\right)
* \]
* which is not the identity matrix. Therefore, despite \(A \times V = V \times D\),
* \(A \ne V \times D \time V^T\), which would hold for real eigendecomposition.
* </p>
* <p>
* Also note that for consistency with Wolfram langage
* <a href="https://reference.wolfram.com/language/ref/Eigenvectors.html">eigenvectors</a>,
* we add zero vectors when the geometric multiplicity of the eigenvalue is smaller than
* its algebraic multiplicity (hence the regular eigenvector matrix should be non-square).
* With these additional null vectors, the eigenvectors matrix becomes square. This happens
* for example with the square matrix
* \[
* A = \left(\begin{matrix}
* 1 & 0 & 0\\
* -2 & 1 & 0\\
* 0 & 0 & 1
* \end{matrix}\right)
* \]
* Its characteristic polynomial is \((1-\lambda)^3\), hence is has one eigen value \(\lambda=1\)
* with algebraic multiplicity 3. However, this eigenvalue leads to only two eigenvectors
* \(v_1=(0, 1, 0)\) and \(v_2=(0, 0, 1)\), hence its geometric multiplicity is only 2, not 3.
* So we add a third zero vector \(v_3=(0, 0, 0)\), in the same way Wolfram language does.
* </p>
*
* Compute complex eigen values from the Schur transform. Compute complex eigen
* vectors based on eigen values and the inverse iteration method.
*
* see: <a href="https://en.wikipedia.org/wiki/Inverse_iteration">Inverse iteration</a>
* <a href="https://en.wikiversity.org/wiki/Shifted_inverse_iteration">Shifted inverse iteration</a>
* <a href="https://www.robots.ox.ac.uk/~sjrob/Teaching/EngComp/ecl4.pdf">Computation of matrix eigenvalues and eigenvectors</a>
*
*/
public class ComplexEigenDecomposition {
/** Default threshold below which eigenvectors are considered equal. */
public static final double DEFAULT_EIGENVECTORS_EQUALITY = 1.0e-5;
/** Default value to use for internal epsilon. */
public static final double DEFAULT_EPSILON = 1e-12;
/** Internally used epsilon criteria for final AV=VD check. */
public static final double DEFAULT_EPSILON_AV_VD_CHECK = 1e-6;
/** Maximum number of inverse iterations. */
private static final int MAX_ITER = 10;
/** complex eigenvalues. */
private Complex[] eigenvalues;
/** Eigenvectors. */
private FieldVector<Complex>[] eigenvectors;
/** Cached value of V. */
private FieldMatrix<Complex> V;
/** Cached value of D. */
private FieldMatrix<Complex> D;
/** Internally used threshold below which eigenvectors are considered equal. */
private final double eigenVectorsEquality;
/** Internally used epsilon criteria. */
private final double epsilon;
/** Internally used epsilon criteria for final AV=VD check. */
private final double epsilonAVVDCheck;
/**
* Constructor for decomposition.
* <p>
* This constructor uses the default values {@link #DEFAULT_EIGENVECTORS_EQUALITY},
* {@link #DEFAULT_EPSILON} and {@link #DEFAULT_EPSILON_AV_VD_CHECK}
* </p>
* @param matrix
* real matrix.
*/
public ComplexEigenDecomposition(final RealMatrix matrix) {
this(matrix, DEFAULT_EIGENVECTORS_EQUALITY,
DEFAULT_EPSILON, DEFAULT_EPSILON_AV_VD_CHECK);
}
/**
* Constructor for decomposition.
* <p>
* The {@code eigenVectorsEquality} threshold is used to ensure the L∞-normalized
* eigenvectors found using inverse iteration are different from each other.
* if \(min(|e_i-e_j|,|e_i+e_j|)\) is smaller than this threshold, the algorithm
* considers it has found again an already known vector, so it drops it and attempts
* a new inverse iteration with a different start vector. This value should be
* much larger than {@code epsilon} which is used for convergence
* </p>
* @param matrix real matrix.
* @param eigenVectorsEquality threshold below which eigenvectors are considered equal
* @param epsilon Epsilon used for internal tests (e.g. is singular, eigenvalue ratio, etc.)
* @param epsilonAVVDCheck Epsilon criteria for final AV=VD check
* @since 1.8
*/
public ComplexEigenDecomposition(final RealMatrix matrix, final double eigenVectorsEquality,
final double epsilon, final double epsilonAVVDCheck) {
if (!matrix.isSquare()) {
throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SQUARE_MATRIX,
matrix.getRowDimension(), matrix.getColumnDimension());
}
this.eigenVectorsEquality = eigenVectorsEquality;
this.epsilon = epsilon;
this.epsilonAVVDCheck = epsilonAVVDCheck;
// computing the eigen values
findEigenValues(matrix);
// computing the eigen vectors
findEigenVectors(convertToFieldComplex(matrix));
// V
final int m = eigenvectors.length;
V = MatrixUtils.createFieldMatrix(ComplexField.getInstance(), m, m);
for (int k = 0; k < m; ++k) {
V.setColumnVector(k, eigenvectors[k]);
}
// D
D = MatrixUtils.createFieldDiagonalMatrix(eigenvalues);
checkDefinition(matrix);
}
/**
* Getter of the eigen values.
*
* @return eigen values.
*/
public Complex[] getEigenvalues() {
return eigenvalues.clone();
}
/**
* Getter of the eigen vectors.
*
* @param i
* which eigen vector.
* @return eigen vector.
*/
public FieldVector<Complex> getEigenvector(final int i) {
return eigenvectors[i].copy();
}
/** Reset eigenvalues and eigen vectors from matrices.
* <p>
* This method is intended to be called by sub-classes (mainly {@link OrderedComplexEigenDecomposition})
* that reorder the matrices elements. It rebuild the eigenvalues and eigen vectors arrays
* from the D and V matrices.
* </p>
* @since 2.1
*/
protected void matricesToEigenArrays() {
for (int i = 0; i < eigenvalues.length; ++i) {
eigenvalues[i] = D.getEntry(i, i);
}
for (int i = 0; i < eigenvectors.length; ++i) {
for (int j = 0; j < eigenvectors[i].getDimension(); ++j) {
eigenvectors[i].setEntry(j, V.getEntry(j, i));
}
}
}
/**
* Confirm if there are complex eigen values.
*
* @return true if there are complex eigen values.
*/
public boolean hasComplexEigenvalues() {
for (int i = 0; i < eigenvalues.length; i++) {
if (!Precision.equals(eigenvalues[i].getImaginary(), 0.0, epsilon)) {
return true;
}
}
return false;
}
/**
* Computes the determinant.
*
* @return the determinant.
*/
public double getDeterminant() {
Complex determinant = new Complex(1, 0);
for (Complex lambda : eigenvalues) {
determinant = determinant.multiply(lambda);
}
return determinant.getReal();
}
/**
* Getter V.
*
* @return V.
*/
public FieldMatrix<Complex> getV() {
return V;
}
/**
* Getter D.
*
* @return D.
*/
public FieldMatrix<Complex> getD() {
return D;
}
/**
* Getter VT.
*
* @return VT.
*/
public FieldMatrix<Complex> getVT() {
return V.transpose();
}
/**
* Compute eigen values using the Schur transform.
*
* @param matrix
* real matrix to compute eigen values.
*/
protected void findEigenValues(final RealMatrix matrix) {
final SchurTransformer schurTransform = new SchurTransformer(matrix);
final double[][] matT = schurTransform.getT().getData();
eigenvalues = new Complex[matT.length];
for (int i = 0; i < eigenvalues.length; i++) {
if (i == (eigenvalues.length - 1) || Precision.equals(matT[i + 1][i], 0.0, epsilon)) {
eigenvalues[i] = new Complex(matT[i][i]);
} else {
final double x = matT[i + 1][i + 1];
final double p = 0.5 * (matT[i][i] - x);
final double z = FastMath.sqrt(FastMath.abs(p * p + matT[i + 1][i] * matT[i][i + 1]));
eigenvalues[i] = new Complex(x + p, z);
eigenvalues[i + 1] = new Complex(x + p, -z);
i++;
}
}
}
/**
* Compute the eigen vectors using the inverse power method.
*
* @param matrix
* real matrix to compute eigen vectors.
*/
@SuppressWarnings("unchecked")
protected void findEigenVectors(final FieldMatrix<Complex> matrix) {
// number of eigen values/vectors
int n = eigenvalues.length;
// eigen vectors
eigenvectors = (FieldVector<Complex>[]) Array.newInstance(FieldVector.class, n);
// computing eigen vector based on eigen values and inverse iteration
for (int i = 0; i < eigenvalues.length; i++) {
// shifted non-singular matrix matrix A-(λ+ε)I that is close to the singular matrix A-λI
Complex mu = eigenvalues[i].add(epsilon);
final FieldMatrix<Complex> shifted = matrix.copy();
for (int k = 0; k < matrix.getColumnDimension(); ++k) {
shifted.setEntry(k, k, shifted.getEntry(k, k).subtract(mu));
}
// solver for linear system (A - (λ+ε)I) Bₖ₊₁ = Bₖ
FieldDecompositionSolver<Complex> solver = new FieldQRDecomposition<>(shifted).getSolver();
// loop over possible start vectors
for (int p = 0; eigenvectors[i] == null && p < matrix.getColumnDimension(); ++p) {
// find a vector to start iterations
FieldVector<Complex> b = findStart(p);
if (getNorm(b).norm() > Precision.SAFE_MIN) {
// start vector is a good candidate for inverse iteration
// perform inverse iteration
double delta = Double.POSITIVE_INFINITY;
for (int k = 0; delta > epsilon && k < MAX_ITER; k++) {
// solve (A - (λ+ε)) Bₖ₊₁ = Bₖ
final FieldVector<Complex> bNext = solver.solve(b);
// normalize according to L∞ norm
normalize(bNext);
// compute convergence criterion, comparing Bₖ and both ±Bₖ₊₁
// as iterations sometimes flip between two opposite vectors
delta = separation(b, bNext);
// prepare next iteration
b = bNext;
}
// check we have not found again an already known vector
for (int j = 0; b != null && j < i; ++j) {
if (separation(eigenvectors[j], b) <= eigenVectorsEquality) {
// the selected start vector leads us to found a known vector again,
// we must try another start
b = null;
}
}
eigenvectors[i] = b;
}
}
if (eigenvectors[i] == null) {
// for consistency with Wolfram langage
// https://reference.wolfram.com/language/ref/Eigenvectors.html
// we add zero vectors when the geometric multiplicity of the eigenvalue
// is smaller than its algebraic multiplicity (hence the regular eigenvector
// matrix should be non-square). With these additional null vectors, the
// eigenvectors matrix becomes square
eigenvectors[i] = MatrixUtils.createFieldVector(ComplexField.getInstance(), n);
}
}
}
/** Find a start vector orthogonal to all already found normalized eigenvectors.
* @param index index of the vector
* @return start vector
*/
private FieldVector<Complex> findStart(final int index) {
// create vector
final FieldVector<Complex> start =
MatrixUtils.createFieldVector(ComplexField.getInstance(),
eigenvalues.length);
// initialize with a canonical vector
start.setEntry(index, Complex.ONE);
return start;
}
/**
* Compute the L∞ norm of the a given vector.
*
* @param vector
* vector.
* @return L∞ norm.
*/
private Complex getNorm(FieldVector<Complex> vector) {
double normR = 0;
Complex normC = Complex.ZERO;
for (int i = 0; i < vector.getDimension(); i++) {
final Complex ci = vector.getEntry(i);
final double ni = FastMath.hypot(ci.getReal(), ci.getImaginary());
if (ni > normR) {
normR = ni;
normC = ci;
}
}
return normC;
}
/** Normalize a vector with respect to L∞ norm.
* @param v vector to normalized
*/
private void normalize(final FieldVector<Complex> v) {
final Complex invNorm = getNorm(v).reciprocal();
for (int j = 0; j < v.getDimension(); ++j) {
v.setEntry(j, v.getEntry(j).multiply(invNorm));
}
}
/** Compute the separation between two normalized vectors (which may be in opposite directions).
* @param v1 first normalized vector
* @param v2 second normalized vector
* @return min (|v1 - v2|, |v1+v2|)
*/
private double separation(final FieldVector<Complex> v1, final FieldVector<Complex> v2) {
double deltaPlus = 0;
double deltaMinus = 0;
for (int j = 0; j < v1.getDimension(); ++j) {
final Complex bCurrj = v1.getEntry(j);
final Complex bNextj = v2.getEntry(j);
deltaPlus = FastMath.max(deltaPlus,
FastMath.hypot(bNextj.getReal() + bCurrj.getReal(),
bNextj.getImaginary() + bCurrj.getImaginary()));
deltaMinus = FastMath.max(deltaMinus,
FastMath.hypot(bNextj.getReal() - bCurrj.getReal(),
bNextj.getImaginary() - bCurrj.getImaginary()));
}
return FastMath.min(deltaPlus, deltaMinus);
}
/**
* Check definition of the decomposition in runtime.
*
* @param matrix
* matrix to be decomposed.
*/
protected void checkDefinition(final RealMatrix matrix) {
FieldMatrix<Complex> matrixC = convertToFieldComplex(matrix);
// checking definition of the decomposition
// testing A*V = V*D
FieldMatrix<Complex> AV = matrixC.multiply(getV());
FieldMatrix<Complex> VD = getV().multiply(getD());
if (!equalsWithPrecision(AV, VD, epsilonAVVDCheck)) {
throw new MathRuntimeException(LocalizedCoreFormats.FAILED_DECOMPOSITION,
matrix.getRowDimension(), matrix.getColumnDimension());
}
}
/**
* Helper method that checks with two matrix is equals taking into account a
* given precision.
*
* @param matrix1 first matrix to compare
* @param matrix2 second matrix to compare
* @param tolerance tolerance on matrices entries
* @return true is matrices entries are equal within tolerance,
* false otherwise
*/
private boolean equalsWithPrecision(final FieldMatrix<Complex> matrix1,
final FieldMatrix<Complex> matrix2, final double tolerance) {
boolean toRet = true;
for (int i = 0; i < matrix1.getRowDimension(); i++) {
for (int j = 0; j < matrix1.getColumnDimension(); j++) {
Complex c1 = matrix1.getEntry(i, j);
Complex c2 = matrix2.getEntry(i, j);
if (c1.add(c2.negate()).norm() > tolerance) {
toRet = false;
break;
}
}
}
return toRet;
}
/**
* It converts a real matrix into a complex field matrix.
*
* @param matrix
* real matrix.
* @return complex matrix.
*/
private FieldMatrix<Complex> convertToFieldComplex(RealMatrix matrix) {
final FieldMatrix<Complex> toRet =
MatrixUtils.createFieldIdentityMatrix(ComplexField.getInstance(),
matrix.getRowDimension());
for (int i = 0; i < toRet.getRowDimension(); i++) {
for (int j = 0; j < toRet.getColumnDimension(); j++) {
toRet.setEntry(i, j, new Complex(matrix.getEntry(i, j)));
}
}
return toRet;
}
}