|
|
Java example source code file (EigenDecomposition.java)
This example Java source code file (EigenDecomposition.java) is included in the alvinalexander.com
"Java Source Code
Warehouse" project. The intent of this project is to help you "Learn
Java by Example" TM.
Learn more about this Java project at its project page.
The EigenDecomposition.java Java example source code
/*
* 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.
* A, V and D are all m × m matrices.</p>
* <p>This class is similar in spirit to the EigenvalueDecomposition
* class from the <a href="http://math.nist.gov/javanumerics/jama/">JAMA
* library, with the following changes:</p>
* <ul>
* <li>a {@link #getVT() getVt} method has been added,
* <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>a {@link #getSolver() getSolver} method has been added.
* </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
* @see <a href="http://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix">Wikipedia
* @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 eigenvalue of the original
* matrix.
*
* @param i index of the eigenvalue (counting from 0)
* @return real part of the i<sup>th 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 eigenvalue of the original
* matrix.
*
* @param i Index of the eigenvalue (counting from 0).
* @return the imaginary part of the i<sup>th 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 eigenvector of the original matrix.
*
* @param i Index of the eigenvector (counting from 0).
* @return a copy of the i<sup>th 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);
}
}
}
Other Java examples (source code examples)
Here is a short list of links related to this Java EigenDecomposition.java source code file:
|
The search page
Other Java source code examples at this package level
Click here to learn more about this project
