org.apache.commons.math.linear.EigenDecompositionImpl Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of commons-math Show documentation
Show all versions of commons-math Show documentation
The Math project is a library of lightweight, self-contained mathematics and statistics components addressing the most common practical problems not immediately available in the Java programming language or commons-lang.
The newest version!
/*
* 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.math.linear;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.util.MathUtils;
import org.apache.commons.math.util.FastMath;
/**
* Calculates the eigen decomposition of a real symmetric
* matrix.
*
* The eigen decomposition of matrix A is a set of two matrices: V and D such
* that A = V D VT. A, V and D are all m × m matrices.
*
*
* As of 2.0, this class supports only symmetric matrices, and
* hence computes only real realEigenvalues. This implies the D matrix returned
* by {@link #getD()} is always diagonal and the imaginary values returned
* {@link #getImagEigenvalue(int)} and {@link #getImagEigenvalues()} are always
* null.
*
*
* When called with a {@link RealMatrix} argument, this implementation only uses
* the upper part of the matrix, the part below the diagonal is not accessed at
* all.
*
*
* 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
*
* @version $Revision: 1002040 $ $Date: 2010-09-28 09:18:31 +0200 (mar. 28 sept. 2010) $
* @since 2.0
*/
public class EigenDecompositionImpl implements EigenDecomposition {
/** 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;
/**
* Calculates the eigen decomposition of the given symmetric matrix.
* @param matrix The symmetric matrix to decompose.
* @param splitTolerance dummy parameter, present for backward compatibility only.
* @exception InvalidMatrixException (wrapping a
* {@link org.apache.commons.math.ConvergenceException} if algorithm
* fails to converge
*/
public EigenDecompositionImpl(final RealMatrix matrix,final double splitTolerance)
throws InvalidMatrixException {
if (isSymmetric(matrix)) {
transformToTridiagonal(matrix);
findEigenVectors(transformer.getQ().getData());
} else {
// as of 2.0, non-symmetric matrices (i.e. complex eigenvalues) are
// NOT supported
// see issue https://issues.apache.org/jira/browse/MATH-235
throw new InvalidMatrixException(
LocalizedFormats.ASSYMETRIC_EIGEN_NOT_SUPPORTED);
}
}
/**
* 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 triadiagonal form
* @param secondary Secondary of the tridiagonal form
* @param splitTolerance dummy parameter, present for backward compatibility only.
* @exception InvalidMatrixException (wrapping a
* {@link org.apache.commons.math.ConvergenceException} if algorithm
* fails to converge
*/
public EigenDecompositionImpl(final double[] main,final double[] secondary,
final double splitTolerance)
throws InvalidMatrixException {
this.main = main.clone();
this.secondary = secondary.clone();
transformer = null;
final int size=main.length;
double[][] z = new double[size][size];
for (int i=0;i
(FastMath.max(FastMath.abs(mij), FastMath.abs(mji)) * eps)) {
return false;
}
}
}
return true;
}
/** {@inheritDoc} */
public RealMatrix getV() throws InvalidMatrixException {
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;
}
/** {@inheritDoc} */
public RealMatrix getD() throws InvalidMatrixException {
if (cachedD == null) {
// cache the matrix for subsequent calls
cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues);
}
return cachedD;
}
/** {@inheritDoc} */
public RealMatrix getVT() throws InvalidMatrixException {
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;
}
/** {@inheritDoc} */
public double[] getRealEigenvalues() throws InvalidMatrixException {
return realEigenvalues.clone();
}
/** {@inheritDoc} */
public double getRealEigenvalue(final int i) throws InvalidMatrixException,
ArrayIndexOutOfBoundsException {
return realEigenvalues[i];
}
/** {@inheritDoc} */
public double[] getImagEigenvalues() throws InvalidMatrixException {
return imagEigenvalues.clone();
}
/** {@inheritDoc} */
public double getImagEigenvalue(final int i) throws InvalidMatrixException,
ArrayIndexOutOfBoundsException {
return imagEigenvalues[i];
}
/** {@inheritDoc} */
public RealVector getEigenvector(final int i)
throws InvalidMatrixException, ArrayIndexOutOfBoundsException {
return eigenvectors[i].copy();
}
/**
* Return the determinant of the matrix
* @return determinant of the matrix
*/
public double getDeterminant() {
double determinant = 1;
for (double lambda : realEigenvalues) {
determinant *= lambda;
}
return determinant;
}
/** {@inheritDoc} */
public DecompositionSolver getSolver() {
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;
/**
* Build 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;
}
/**
* Solve the linear equation A × X = B for symmetric matrices A.
*
* This method only find exact linear solutions, i.e. solutions for
* which ||A × X - B|| is exactly 0.
*
* @param b
* right-hand side of the equation A × X = B
* @return a vector X that minimizes the two norm of A × X - B
* @exception IllegalArgumentException
* if matrices dimensions don't match
* @exception InvalidMatrixException
* if decomposed matrix is singular
*/
public double[] solve(final double[] b)
throws IllegalArgumentException, InvalidMatrixException {
if (!isNonSingular()) {
throw new SingularMatrixException();
}
final int m = realEigenvalues.length;
if (b.length != m) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.VECTOR_LENGTH_MISMATCH,
b.length, 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 bp;
}
/**
* Solve the linear equation A × X = B for symmetric matrices A.
*
* This method only find exact linear solutions, i.e. solutions for
* which ||A × X - B|| is exactly 0.
*
* @param b
* right-hand side of the equation A × X = B
* @return a vector X that minimizes the two norm of A × X - B
* @exception IllegalArgumentException
* if matrices dimensions don't match
* @exception InvalidMatrixException
* if decomposed matrix is singular
*/
public RealVector solve(final RealVector b)
throws IllegalArgumentException, InvalidMatrixException {
if (!isNonSingular()) {
throw new SingularMatrixException();
}
final int m = realEigenvalues.length;
if (b.getDimension() != m) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.VECTOR_LENGTH_MISMATCH, 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);
}
/**
* Solve the linear equation A × X = B for symmetric matrices A.
*
* This method only find exact linear solutions, i.e. solutions for
* which ||A × X - B|| is exactly 0.
*
* @param b
* right-hand side of the equation A × X = B
* @return a matrix X that minimizes the two norm of A × X - B
* @exception IllegalArgumentException
* if matrices dimensions don't match
* @exception InvalidMatrixException
* if decomposed matrix is singular
*/
public RealMatrix solve(final RealMatrix b)
throws IllegalArgumentException, InvalidMatrixException {
if (!isNonSingular()) {
throw new SingularMatrixException();
}
final int m = realEigenvalues.length;
if (b.getRowDimension() != m) {
throw MathRuntimeException
.createIllegalArgumentException(
LocalizedFormats.DIMENSIONS_MISMATCH_2x2,
b.getRowDimension(), b.getColumnDimension(), m,
"n");
}
final int nColB = b.getColumnDimension();
final double[][] bp = new double[m][nColB];
for (int k = 0; k < nColB; ++k) {
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) * b.getEntry(j, k);
}
s /= realEigenvalues[i];
for (int j = 0; j < m; ++j) {
bp[j][k] += s * vData[j];
}
}
}
return MatrixUtils.createRealMatrix(bp);
}
/**
* Check if the decomposed matrix is non-singular.
* @return true if the decomposed matrix is non-singular
*/
public boolean isNonSingular() {
for (int i = 0; i < realEigenvalues.length; ++i) {
if ((realEigenvalues[i] == 0) && (imagEigenvalues[i] == 0)) {
return false;
}
}
return true;
}
/**
* Get the inverse of the decomposed matrix.
* @return inverse matrix
* @throws InvalidMatrixException
* if decomposed matrix is singular
*/
public RealMatrix getInverse() throws InvalidMatrixException {
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);
}
}
/**
* Transform matrix to tridiagonal.
* @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 tri-diagonal form.
*/
private void findEigenVectors(double[][] householderMatrix) {
double[][]z = householderMatrix.clone();
final int n = main.length;
realEigenvalues = new double[n];
imagEigenvalues = new double[n];
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.0;
// Determine the largest main and secondary value in absolute term.
double maxAbsoluteValue=0.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.0) {
for (int i=0; i < n; i++) {
if (FastMath.abs(realEigenvalues[i])<=MathUtils.EPSILON*maxAbsoluteValue) {
realEigenvalues[i]=0.0;
}
if (FastMath.abs(e[i])<=MathUtils.EPSILON*maxAbsoluteValue) {
e[i]=0.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 InvalidMatrixException(
new MaxIterationsExceededException(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 = c * s;
} else {
s = p / q;
t = FastMath.sqrt(s * s + 1.0);
e[i + 1] = q * t;
c = 1.0 / t;
s = 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.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])