/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF 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
*
* http://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.apache.commons.math3.linear;
import org.apache.commons.math3.complex.Complex;
import org.apache.commons.math3.exception.MathArithmeticException;
import org.apache.commons.math3.exception.MathUnsupportedOperationException;
import org.apache.commons.math3.exception.MaxCountExceededException;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.util.Precision;
import org.apache.commons.math3.util.FastMath;
/**
* Calculates the eigen decomposition of a real matrix.
* <p>The eigen decomposition of matrix A is a set of two matrices:
* V and D such that A = V × D × V<sup>T</sup>.
* A, V and D are all m × m matrices.</p>
* <p>This class is similar in spirit to the <code>EigenvalueDecomposition</code>
* class from the <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a>
* library, with the following changes:</p>
* <ul>
* <li>a {@link #getVT() getVt} method has been added,</li>
* <li>two {@link #getRealEigenvalue(int) getRealEigenvalue} and {@link #getImagEigenvalue(int)
* getImagEigenvalue} methods to pick up a single eigenvalue have been added,</li>
* <li>a {@link #getEigenvector(int) getEigenvector} method to pick up a single
* eigenvector has been added,</li>
* <li>a {@link #getDeterminant() getDeterminant} method has been added.</li>
* <li>a {@link #getSolver() getSolver} method has been added.</li>
* </ul>
* <p>
* As of 3.1, this class supports general real matrices (both symmetric and non-symmetric):
* </p>
* <p>
* If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is diagonal and the eigenvector
* matrix V is orthogonal, i.e. A = V.multiply(D.multiply(V.transpose())) and
* V.multiply(V.transpose()) equals the identity matrix.
* </p>
* <p>
* If A is not symmetric, then the eigenvalue matrix D is block diagonal with the real eigenvalues
* in 1-by-1 blocks and any complex eigenvalues, lambda + i*mu, in 2-by-2 blocks:
* <pre>
* [lambda, mu ]
* [ -mu, lambda]
* </pre>
* The columns of V represent the eigenvectors in the sense that A*V = V*D,
* i.e. A.multiply(V) equals V.multiply(D).
* The matrix V may be badly conditioned, or even singular, so the validity of the equation
* A = V*D*inverse(V) depends upon the condition of V.
* </p>
* <p>
* This implementation is based on the paper by A. Drubrulle, R.S. Martin and
* J.H. Wilkinson "The Implicit QL Algorithm" in Wilksinson and Reinsch (1971)
* Handbook for automatic computation, vol. 2, Linear algebra, Springer-Verlag,
* New-York
* </p>
* @see <a href="http://mathworld.wolfram.com/EigenDecomposition.html">MathWorld</a>
* @see <a href="http://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix">Wikipedia</a>
* @since 2.0 (changed to concrete class in 3.0)
*/
public class EigenDecomposition {
/** Internally used epsilon criteria. */
private static final double EPSILON = 1e-12;
/** Maximum number of iterations accepted in the implicit QL transformation */
private byte maxIter = 30;
/** Main diagonal of the tridiagonal matrix. */
private double[] main;
/** Secondary diagonal of the tridiagonal matrix. */
private double[] secondary;
/**
* Transformer to tridiagonal (may be null if matrix is already
* tridiagonal).
*/
private TriDiagonalTransformer transformer;
/** Real part of the realEigenvalues. */
private double[] realEigenvalues;
/** Imaginary part of the realEigenvalues. */
private double[] imagEigenvalues;
/** Eigenvectors. */
private ArrayRealVector[] eigenvectors;
/** Cached value of V. */
private RealMatrix cachedV;
/** Cached value of D. */
private RealMatrix cachedD;
/** Cached value of Vt. */
private RealMatrix cachedVt;
/** Whether the matrix is symmetric. */
private final boolean isSymmetric;
/**
* Calculates the eigen decomposition of the given real matrix.
* <p>
* Supports decomposition of a general matrix since 3.1.
*
* @param matrix Matrix to decompose.
* @throws MaxCountExceededException if the algorithm fails to converge.
* @throws MathArithmeticException if the decomposition of a general matrix
* results in a matrix with zero norm
* @since 3.1
*/
public EigenDecomposition(final RealMatrix matrix)
throws MathArithmeticException {
final double symTol = 10 * matrix.getRowDimension() * matrix.getColumnDimension() * Precision.EPSILON;
isSymmetric = MatrixUtils.isSymmetric(matrix, symTol);
if (isSymmetric) {
transformToTridiagonal(matrix);
findEigenVectors(transformer.getQ().getData());
} else {
final SchurTransformer t = transformToSchur(matrix);
findEigenVectorsFromSchur(t);
}
}
/**
* Calculates the eigen decomposition of the given real matrix.
*
* @param matrix Matrix to decompose.
* @param splitTolerance Dummy parameter (present for backward
* compatibility only).
* @throws MathArithmeticException if the decomposition of a general matrix
* results in a matrix with zero norm
* @throws MaxCountExceededException if the algorithm fails to converge.
* @deprecated in 3.1 (to be removed in 4.0) due to unused parameter
*/
@Deprecated
public EigenDecomposition(final RealMatrix matrix,
final double splitTolerance)
throws MathArithmeticException {
this(matrix);
}
/**
* Calculates the eigen decomposition of the symmetric tridiagonal
* matrix. The Householder matrix is assumed to be the identity matrix.
*
* @param main Main diagonal of the symmetric tridiagonal form.
* @param secondary Secondary of the tridiagonal form.
* @throws MaxCountExceededException if the algorithm fails to converge.
* @since 3.1
*/
public EigenDecomposition(final double[] main, final double[] secondary) {
isSymmetric = true;
this.main = main.clone();
this.secondary = secondary.clone();
transformer = null;
final int size = main.length;
final double[][] z = new double[size][size];
for (int i = 0; i < size; i++) {
z[i][i] = 1.0;
}
findEigenVectors(z);
}
/**
* Calculates the eigen decomposition of the symmetric tridiagonal
* matrix. The Householder matrix is assumed to be the identity matrix.
*
* @param main Main diagonal of the symmetric tridiagonal form.
* @param secondary Secondary of the tridiagonal form.
* @param splitTolerance Dummy parameter (present for backward
* compatibility only).
* @throws MaxCountExceededException if the algorithm fails to converge.
* @deprecated in 3.1 (to be removed in 4.0) due to unused parameter
*/
@Deprecated
public EigenDecomposition(final double[] main, final double[] secondary,
final double splitTolerance) {
this(main, secondary);
}
/**
* Gets the matrix V of the decomposition.
* V is an orthogonal matrix, i.e. its transpose is also its inverse.
* The columns of V are the eigenvectors of the original matrix.
* No assumption is made about the orientation of the system axes formed
* by the columns of V (e.g. in a 3-dimension space, V can form a left-
* or right-handed system).
*
* @return the V matrix.
*/
public RealMatrix getV() {
if (cachedV == null) {
final int m = eigenvectors.length;
cachedV = MatrixUtils.createRealMatrix(m, m);
for (int k = 0; k < m; ++k) {
cachedV.setColumnVector(k, eigenvectors[k]);
}
}
// return the cached matrix
return cachedV;
}
/**
* Gets the block diagonal matrix D of the decomposition.
* D is a block diagonal matrix.
* Real eigenvalues are on the diagonal while complex values are on
* 2x2 blocks { {real +imaginary}, {-imaginary, real} }.
*
* @return the D matrix.
*
* @see #getRealEigenvalues()
* @see #getImagEigenvalues()
*/
public RealMatrix getD() {
if (cachedD == null) {
// cache the matrix for subsequent calls
cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues);
for (int i = 0; i < imagEigenvalues.length; i++) {
if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) > 0) {
cachedD.setEntry(i, i+1, imagEigenvalues[i]);
} else if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) {
cachedD.setEntry(i, i-1, imagEigenvalues[i]);
}
}
}
return cachedD;
}
/**
* Gets the transpose of the matrix V of the decomposition.
* V is an orthogonal matrix, i.e. its transpose is also its inverse.
* The columns of V are the eigenvectors of the original matrix.
* No assumption is made about the orientation of the system axes formed
* by the columns of V (e.g. in a 3-dimension space, V can form a left-
* or right-handed system).
*
* @return the transpose of the V matrix.
*/
public RealMatrix getVT() {
if (cachedVt == null) {
final int m = eigenvectors.length;
cachedVt = MatrixUtils.createRealMatrix(m, m);
for (int k = 0; k < m; ++k) {
cachedVt.setRowVector(k, eigenvectors[k]);
}
}
// return the cached matrix
return cachedVt;
}
/**
* Returns whether the calculated eigen values are complex or real.
* <p>The method performs a zero check for each element of the
* {@link #getImagEigenvalues()} array and returns {@code true} if any
* element is not equal to zero.
*
* @return {@code true} if the eigen values are complex, {@code false} otherwise
* @since 3.1
*/
public boolean hasComplexEigenvalues() {
for (int i = 0; i < imagEigenvalues.length; i++) {
if (!Precision.equals(imagEigenvalues[i], 0.0, EPSILON)) {
return true;
}
}
return false;
}
/**
* Gets a copy of the real parts of the eigenvalues of the original matrix.
*
* @return a copy of the real parts of the eigenvalues of the original matrix.
*
* @see #getD()
* @see #getRealEigenvalue(int)
* @see #getImagEigenvalues()
*/
public double[] getRealEigenvalues() {
return realEigenvalues.clone();
}
/**
* Returns the real part of the i<sup>th</sup> eigenvalue of the original
* matrix.
*
* @param i index of the eigenvalue (counting from 0)
* @return real part of the i<sup>th</sup> eigenvalue of the original
* matrix.
*
* @see #getD()
* @see #getRealEigenvalues()
* @see #getImagEigenvalue(int)
*/
public double getRealEigenvalue(final int i) {
return realEigenvalues[i];
}
/**
* Gets a copy of the imaginary parts of the eigenvalues of the original
* matrix.
*
* @return a copy of the imaginary parts of the eigenvalues of the original
* matrix.
*
* @see #getD()
* @see #getImagEigenvalue(int)
* @see #getRealEigenvalues()
*/
public double[] getImagEigenvalues() {
return imagEigenvalues.clone();
}
/**
* Gets the imaginary part of the i<sup>th</sup> eigenvalue of the original
* matrix.
*
* @param i Index of the eigenvalue (counting from 0).
* @return the imaginary part of the i<sup>th</sup> eigenvalue of the original
* matrix.
*
* @see #getD()
* @see #getImagEigenvalues()
* @see #getRealEigenvalue(int)
*/
public double getImagEigenvalue(final int i) {
return imagEigenvalues[i];
}
/**
* Gets a copy of the i<sup>th</sup> eigenvector of the original matrix.
*
* @param i Index of the eigenvector (counting from 0).
* @return a copy of the i<sup>th</sup> eigenvector of the original matrix.
* @see #getD()
*/
public RealVector getEigenvector(final int i) {
return eigenvectors[i].copy();
}
/**
* Computes the determinant of the matrix.
*
* @return the determinant of the matrix.
*/
public double getDeterminant() {
double determinant = 1;
for (double lambda : realEigenvalues) {
determinant *= lambda;
}
return determinant;
}
/**
* Computes the square-root of the matrix.
* This implementation assumes that the matrix is symmetric and positive
* definite.
*
* @return the square-root of the matrix.
* @throws MathUnsupportedOperationException if the matrix is not
* symmetric or not positive definite.
* @since 3.1
*/
public RealMatrix getSquareRoot() {
if (!isSymmetric) {
throw new MathUnsupportedOperationException();
}
final double[] sqrtEigenValues = new double[realEigenvalues.length];
for (int i = 0; i < realEigenvalues.length; i++) {
final double eigen = realEigenvalues[i];
if (eigen <= 0) {
throw new MathUnsupportedOperationException();
}
sqrtEigenValues[i] = FastMath.sqrt(eigen);
}
final RealMatrix sqrtEigen = MatrixUtils.createRealDiagonalMatrix(sqrtEigenValues);
final RealMatrix v = getV();
final RealMatrix vT = getVT();
return v.multiply(sqrtEigen).multiply(vT);
}
/**
* Gets a solver for finding the A × X = B solution in exact
* linear sense.
* <p>
* Since 3.1, eigen decomposition of a general matrix is supported,
* but the {@link DecompositionSolver} only supports real eigenvalues.
*
* @return a solver
* @throws MathUnsupportedOperationException if the decomposition resulted in
* complex eigenvalues
*/
public DecompositionSolver getSolver() {
if (hasComplexEigenvalues()) {
throw new MathUnsupportedOperationException();
}
return new Solver(realEigenvalues, imagEigenvalues, eigenvectors);
}
/** Specialized solver. */
private static class Solver implements DecompositionSolver {
/** Real part of the realEigenvalues. */
private double[] realEigenvalues;
/** Imaginary part of the realEigenvalues. */
private double[] imagEigenvalues;
/** Eigenvectors. */
private final ArrayRealVector[] eigenvectors;
/**
* Builds a solver from decomposed matrix.
*
* @param realEigenvalues Real parts of the eigenvalues.
* @param imagEigenvalues Imaginary parts of the eigenvalues.
* @param eigenvectors Eigenvectors.
*/
private Solver(final double[] realEigenvalues,
final double[] imagEigenvalues,
final ArrayRealVector[] eigenvectors) {
this.realEigenvalues = realEigenvalues;
this.imagEigenvalues = imagEigenvalues;
this.eigenvectors = eigenvectors;
}
/**
* Solves the linear equation A × X = B for symmetric matrices A.
* <p>
* This method only finds exact linear solutions, i.e. solutions for
* which ||A × X - B|| is exactly 0.
* </p>
*
* @param b Right-hand side of the equation A × X = B.
* @return a Vector X that minimizes the two norm of A × X - B.
*
* @throws DimensionMismatchException if the matrices dimensions do not match.
* @throws SingularMatrixException if the decomposed matrix is singular.
*/
public RealVector solve(final RealVector b) {
if (!isNonSingular()) {
throw new SingularMatrixException();
}
final int m = realEigenvalues.length;
if (b.getDimension() != m) {
throw new DimensionMismatchException(b.getDimension(), m);
}
final double[] bp = new double[m];
for (int i = 0; i < m; ++i) {
final ArrayRealVector v = eigenvectors[i];
final double[] vData = v.getDataRef();
final double s = v.dotProduct(b) / realEigenvalues[i];
for (int j = 0; j < m; ++j) {
bp[j] += s * vData[j];
}
}
return new ArrayRealVector(bp, false);
}
/** {@inheritDoc} */
public RealMatrix solve(RealMatrix b) {
if (!isNonSingular()) {
throw new SingularMatrixException();
}
final int m = realEigenvalues.length;
if (b.getRowDimension() != m) {
throw new DimensionMismatchException(b.getRowDimension(), m);
}
final int nColB = b.getColumnDimension();
final double[][] bp = new double[m][nColB];
final double[] tmpCol = new double[m];
for (int k = 0; k < nColB; ++k) {
for (int i = 0; i < m; ++i) {
tmpCol[i] = b.getEntry(i, k);
bp[i][k] = 0;
}
for (int i = 0; i < m; ++i) {
final ArrayRealVector v = eigenvectors[i];
final double[] vData = v.getDataRef();
double s = 0;
for (int j = 0; j < m; ++j) {
s += v.getEntry(j) * tmpCol[j];
}
s /= realEigenvalues[i];
for (int j = 0; j < m; ++j) {
bp[j][k] += s * vData[j];
}
}
}
return new Array2DRowRealMatrix(bp, false);
}
/**
* Checks whether the decomposed matrix is non-singular.
*
* @return true if the decomposed matrix is non-singular.
*/
public boolean isNonSingular() {
double largestEigenvalueNorm = 0.0;
// Looping over all values (in case they are not sorted in decreasing
// order of their norm).
for (int i = 0; i < realEigenvalues.length; ++i) {
largestEigenvalueNorm = FastMath.max(largestEigenvalueNorm, eigenvalueNorm(i));
}
// Corner case: zero matrix, all exactly 0 eigenvalues
if (largestEigenvalueNorm == 0.0) {
return false;
}
for (int i = 0; i < realEigenvalues.length; ++i) {
// Looking for eigenvalues that are 0, where we consider anything much much smaller
// than the largest eigenvalue to be effectively 0.
if (Precision.equals(eigenvalueNorm(i) / largestEigenvalueNorm, 0, EPSILON)) {
return false;
}
}
return true;
}
/**
* @param i which eigenvalue to find the norm of
* @return the norm of ith (complex) eigenvalue.
*/
private double eigenvalueNorm(int i) {
final double re = realEigenvalues[i];
final double im = imagEigenvalues[i];
return FastMath.sqrt(re * re + im * im);
}
/**
* Get the inverse of the decomposed matrix.
*
* @return the inverse matrix.
* @throws SingularMatrixException if the decomposed matrix is singular.
*/
public RealMatrix getInverse() {
if (!isNonSingular()) {
throw new SingularMatrixException();
}
final int m = realEigenvalues.length;
final double[][] invData = new double[m][m];
for (int i = 0; i < m; ++i) {
final double[] invI = invData[i];
for (int j = 0; j < m; ++j) {
double invIJ = 0;
for (int k = 0; k < m; ++k) {
final double[] vK = eigenvectors[k].getDataRef();
invIJ += vK[i] * vK[j] / realEigenvalues[k];
}
invI[j] = invIJ;
}
}
return MatrixUtils.createRealMatrix(invData);
}
}
/**
* Transforms the matrix to tridiagonal form.
*
* @param matrix Matrix to transform.
*/
private void transformToTridiagonal(final RealMatrix matrix) {
// transform the matrix to tridiagonal
transformer = new TriDiagonalTransformer(matrix);
main = transformer.getMainDiagonalRef();
secondary = transformer.getSecondaryDiagonalRef();
}
/**
* Find eigenvalues and eigenvectors (Dubrulle et al., 1971)
*
* @param householderMatrix Householder matrix of the transformation
* to tridiagonal form.
*/
private void findEigenVectors(final double[][] householderMatrix) {
final double[][]z = householderMatrix.clone();
final int n = main.length;
realEigenvalues = new double[n];
imagEigenvalues = new double[n];
final double[] e = new double[n];
for (int i = 0; i < n - 1; i++) {
realEigenvalues[i] = main[i];
e[i] = secondary[i];
}
realEigenvalues[n - 1] = main[n - 1];
e[n - 1] = 0;
// Determine the largest main and secondary value in absolute term.
double maxAbsoluteValue = 0;
for (int i = 0; i < n; i++) {
if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) {
maxAbsoluteValue = FastMath.abs(realEigenvalues[i]);
}
if (FastMath.abs(e[i]) > maxAbsoluteValue) {
maxAbsoluteValue = FastMath.abs(e[i]);
}
}
// Make null any main and secondary value too small to be significant
if (maxAbsoluteValue != 0) {
for (int i=0; i < n; i++) {
if (FastMath.abs(realEigenvalues[i]) <= Precision.EPSILON * maxAbsoluteValue) {
realEigenvalues[i] = 0;
}
if (FastMath.abs(e[i]) <= Precision.EPSILON * maxAbsoluteValue) {
e[i]=0;
}
}
}
for (int j = 0; j < n; j++) {
int its = 0;
int m;
do {
for (m = j; m < n - 1; m++) {
double delta = FastMath.abs(realEigenvalues[m]) +
FastMath.abs(realEigenvalues[m + 1]);
if (FastMath.abs(e[m]) + delta == delta) {
break;
}
}
if (m != j) {
if (its == maxIter) {
throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED,
maxIter);
}
its++;
double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]);
double t = FastMath.sqrt(1 + q * q);
if (q < 0.0) {
q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q - t);
} else {
q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q + t);
}
double u = 0.0;
double s = 1.0;
double c = 1.0;
int i;
for (i = m - 1; i >= j; i--) {
double p = s * e[i];
double h = c * e[i];
if (FastMath.abs(p) >= FastMath.abs(q)) {
c = q / p;
t = FastMath.sqrt(c * c + 1.0);
e[i + 1] = p * t;
s = 1.0 / t;
c *= s;
} else {
s = p / q;
t = FastMath.sqrt(s * s + 1.0);
e[i + 1] = q * t;
c = 1.0 / t;
s *= c;
}
if (e[i + 1] == 0.0) {
realEigenvalues[i + 1] -= u;
e[m] = 0.0;
break;
}
q = realEigenvalues[i + 1] - u;
t = (realEigenvalues[i] - q) * s + 2.0 * c * h;
u = s * t;
realEigenvalues[i + 1] = q + u;
q = c * t - h;
for (int ia = 0; ia < n; ia++) {
p = z[ia][i + 1];
z[ia][i + 1] = s * z[ia][i] + c * p;
z[ia][i] = c * z[ia][i] - s * p;
}
}
if (t == 0.0 && i >= j) {
continue;
}
realEigenvalues[j] -= u;
e[j] = q;
e[m] = 0.0;
}
} while (m != j);
}
//Sort the eigen values (and vectors) in increase order
for (int i = 0; i < n; i++) {
int k = i;
double p = realEigenvalues[i];
for (int j = i + 1; j < n; j++) {
if (realEigenvalues[j] > p) {
k = j;
p = realEigenvalues[j];
}
}
if (k != i) {
realEigenvalues[k] = realEigenvalues[i];
realEigenvalues[i] = p;
for (int j = 0; j < n; j++) {
p = z[j][i];
z[j][i] = z[j][k];
z[j][k] = p;
}
}
}
// Determine the largest eigen value in absolute term.
maxAbsoluteValue = 0;
for (int i = 0; i < n; i++) {
if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) {
maxAbsoluteValue=FastMath.abs(realEigenvalues[i]);
}
}
// Make null any eigen value too small to be significant
if (maxAbsoluteValue != 0.0) {
for (int i=0; i < n; i++) {
if (FastMath.abs(realEigenvalues[i]) < Precision.EPSILON * maxAbsoluteValue) {
realEigenvalues[i] = 0;
}
}
}
eigenvectors = new ArrayRealVector[n];
final double[] tmp = new double[n];
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
tmp[j] = z[j][i];
}
eigenvectors[i] = new ArrayRealVector(tmp);
}
}
/**
* Transforms the matrix to Schur form and calculates the eigenvalues.
*
* @param matrix Matrix to transform.
* @return the {@link SchurTransformer Shur transform} for this matrix
*/
private SchurTransformer transformToSchur(final RealMatrix matrix) {
final SchurTransformer schurTransform = new SchurTransformer(matrix);
final double[][] matT = schurTransform.getT().getData();
realEigenvalues = new double[matT.length];
imagEigenvalues = new double[matT.length];
for (int i = 0; i < realEigenvalues.length; i++) {
if (i == (realEigenvalues.length - 1) ||
Precision.equals(matT[i + 1][i], 0.0, EPSILON)) {
realEigenvalues[i] = 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]));
realEigenvalues[i] = x + p;
imagEigenvalues[i] = z;
realEigenvalues[i + 1] = x + p;
imagEigenvalues[i + 1] = -z;
i++;
}
}
return schurTransform;
}
/**
* Performs a division of two complex numbers.
*
* @param xr real part of the first number
* @param xi imaginary part of the first number
* @param yr real part of the second number
* @param yi imaginary part of the second number
* @return result of the complex division
*/
private Complex cdiv(final double xr, final double xi,
final double yr, final double yi) {
return new Complex(xr, xi).divide(new Complex(yr, yi));
}
/**
* Find eigenvectors from a matrix transformed to Schur form.
*
* @param schur the schur transformation of the matrix
* @throws MathArithmeticException if the Schur form has a norm of zero
*/
private void findEigenVectorsFromSchur(final SchurTransformer schur)
throws MathArithmeticException {
final double[][] matrixT = schur.getT().getData();
final double[][] matrixP = schur.getP().getData();
final int n = matrixT.length;
// compute matrix norm
double norm = 0.0;
for (int i = 0; i < n; i++) {
for (int j = FastMath.max(i - 1, 0); j < n; j++) {
norm += FastMath.abs(matrixT[i][j]);
}
}
// we can not handle a matrix with zero norm
if (Precision.equals(norm, 0.0, EPSILON)) {
throw new MathArithmeticException(LocalizedFormats.ZERO_NORM);
}
// Backsubstitute to find vectors of upper triangular form
double r = 0.0;
double s = 0.0;
double z = 0.0;
for (int idx = n - 1; idx >= 0; idx--) {
double p = realEigenvalues[idx];
double q = imagEigenvalues[idx];
if (Precision.equals(q, 0.0)) {
// Real vector
int l = idx;
matrixT[idx][idx] = 1.0;
for (int i = idx - 1; i >= 0; i--) {
double w = matrixT[i][i] - p;
r = 0.0;
for (int j = l; j <= idx; j++) {
r += matrixT[i][j] * matrixT[j][idx];
}
if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) {
z = w;
s = r;
} else {
l = i;
if (Precision.equals(imagEigenvalues[i], 0.0)) {
if (w != 0.0) {
matrixT[i][idx] = -r / w;
} else {
matrixT[i][idx] = -r / (Precision.EPSILON * norm);
}
} else {
// Solve real equations
double x = matrixT[i][i + 1];
double y = matrixT[i + 1][i];
q = (realEigenvalues[i] - p) * (realEigenvalues[i] - p) +
imagEigenvalues[i] * imagEigenvalues[i];
double t = (x * s - z * r) / q;
matrixT[i][idx] = t;
if (FastMath.abs(x) > FastMath.abs(z)) {
matrixT[i + 1][idx] = (-r - w * t) / x;
} else {
matrixT[i + 1][idx] = (-s - y * t) / z;
}
}
// Overflow control
double t = FastMath.abs(matrixT[i][idx]);
if ((Precision.EPSILON * t) * t > 1) {
for (int j = i; j <= idx; j++) {
matrixT[j][idx] /= t;
}
}
}
}
} else if (q < 0.0) {
// Complex vector
int l = idx - 1;
// Last vector component imaginary so matrix is triangular
if (FastMath.abs(matrixT[idx][idx - 1]) > FastMath.abs(matrixT[idx - 1][idx])) {
matrixT[idx - 1][idx - 1] = q / matrixT[idx][idx - 1];
matrixT[idx - 1][idx] = -(matrixT[idx][idx] - p) / matrixT[idx][idx - 1];
} else {
final Complex result = cdiv(0.0, -matrixT[idx - 1][idx],
matrixT[idx - 1][idx - 1] - p, q);
matrixT[idx - 1][idx - 1] = result.getReal();
matrixT[idx - 1][idx] = result.getImaginary();
}
matrixT[idx][idx - 1] = 0.0;
matrixT[idx][idx] = 1.0;
for (int i = idx - 2; i >= 0; i--) {
double ra = 0.0;
double sa = 0.0;
for (int j = l; j <= idx; j++) {
ra += matrixT[i][j] * matrixT[j][idx - 1];
sa += matrixT[i][j] * matrixT[j][idx];
}
double w = matrixT[i][i] - p;
if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) {
z = w;
r = ra;
s = sa;
} else {
l = i;
if (Precision.equals(imagEigenvalues[i], 0.0)) {
final Complex c = cdiv(-ra, -sa, w, q);
matrixT[i][idx - 1] = c.getReal();
matrixT[i][idx] = c.getImaginary();
} else {
// Solve complex equations
double x = matrixT[i][i + 1];
double y = matrixT[i + 1][i];
double vr = (realEigenvalues[i] - p) * (realEigenvalues[i] - p) +
imagEigenvalues[i] * imagEigenvalues[i] - q * q;
final double vi = (realEigenvalues[i] - p) * 2.0 * q;
if (Precision.equals(vr, 0.0) && Precision.equals(vi, 0.0)) {
vr = Precision.EPSILON * norm *
(FastMath.abs(w) + FastMath.abs(q) + FastMath.abs(x) +
FastMath.abs(y) + FastMath.abs(z));
}
final Complex c = cdiv(x * r - z * ra + q * sa,
x * s - z * sa - q * ra, vr, vi);
matrixT[i][idx - 1] = c.getReal();
matrixT[i][idx] = c.getImaginary();
if (FastMath.abs(x) > (FastMath.abs(z) + FastMath.abs(q))) {
matrixT[i + 1][idx - 1] = (-ra - w * matrixT[i][idx - 1] +
q * matrixT[i][idx]) / x;
matrixT[i + 1][idx] = (-sa - w * matrixT[i][idx] -
q * matrixT[i][idx - 1]) / x;
} else {
final Complex c2 = cdiv(-r - y * matrixT[i][idx - 1],
-s - y * matrixT[i][idx], z, q);
matrixT[i + 1][idx - 1] = c2.getReal();
matrixT[i + 1][idx] = c2.getImaginary();
}
}
// Overflow control
double t = FastMath.max(FastMath.abs(matrixT[i][idx - 1]),
FastMath.abs(matrixT[i][idx]));
if ((Precision.EPSILON * t) * t > 1) {
for (int j = i; j <= idx; j++) {
matrixT[j][idx - 1] /= t;
matrixT[j][idx] /= t;
}
}
}
}
}
}
// Back transformation to get eigenvectors of original matrix
for (int j = n - 1; j >= 0; j--) {
for (int i = 0; i <= n - 1; i++) {
z = 0.0;
for (int k = 0; k <= FastMath.min(j, n - 1); k++) {
z += matrixP[i][k] * matrixT[k][j];
}
matrixP[i][j] = z;
}
}
eigenvectors = new ArrayRealVector[n];
final double[] tmp = new double[n];
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
tmp[j] = matrixP[j][i];
}
eigenvectors[i] = new ArrayRealVector(tmp);
}
}
}