Merged EigenDecomposition and EigenDecompositionImpl (see MATH-662).
git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1173965 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Sebastien Brisard
parent
166c691938
commit
041df9cd9b
|
|
@ -17,14 +17,19 @@
|
|||
|
||||
package org.apache.commons.math.linear;
|
||||
|
||||
import org.apache.commons.math.exception.MaxCountExceededException;
|
||||
import org.apache.commons.math.exception.DimensionMismatchException;
|
||||
import org.apache.commons.math.exception.util.LocalizedFormats;
|
||||
import org.apache.commons.math.util.MathUtils;
|
||||
import org.apache.commons.math.util.FastMath;
|
||||
|
||||
/**
|
||||
* An interface to classes that implement an algorithm to calculate the
|
||||
* eigen decomposition of a real matrix.
|
||||
* Calculates the eigen decomposition of a real <strong>symmetric</strong>
|
||||
* 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 interface is similar in spirit to the <code>EigenvalueDecomposition</code>
|
||||
* <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>
|
||||
|
|
@ -36,12 +41,134 @@ package org.apache.commons.math.linear;
|
|||
* <li>a {@link #getDeterminant() getDeterminant} method has been added.</li>
|
||||
* <li>a {@link #getSolver() getSolver} method has been added.</li>
|
||||
* </ul>
|
||||
* <p>
|
||||
* As of 2.0, this class supports only <strong>symmetric</strong> 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.
|
||||
* </p>
|
||||
* <p>
|
||||
* 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.
|
||||
* </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>
|
||||
* @version $Id$
|
||||
* @since 2.0
|
||||
* @since 2.0 (changed to concrete class in 3.0)
|
||||
*/
|
||||
public interface EigenDecomposition {
|
||||
public class 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 Matrix to decompose. It <em>must</em> be symmetric.
|
||||
* @param splitTolerance Dummy parameter (present for backward
|
||||
* compatibility only).
|
||||
* @throws NonSymmetricMatrixException if the matrix is not symmetric.
|
||||
* @throws MaxCountExceededException if the algorithm fails to converge.
|
||||
*/
|
||||
public EigenDecomposition(final RealMatrix matrix,
|
||||
final double splitTolerance) {
|
||||
if (isSymmetric(matrix, true)) {
|
||||
transformToTridiagonal(matrix);
|
||||
findEigenVectors(transformer.getQ().getData());
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* 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).
|
||||
* @throws MaxCountExceededException if the algorithm fails to converge.
|
||||
*/
|
||||
public EigenDecomposition(final double[] main,final double[] secondary,
|
||||
final double splitTolerance) {
|
||||
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<size;i++) {
|
||||
z[i][i]=1.0;
|
||||
}
|
||||
findEigenVectors(z);
|
||||
}
|
||||
|
||||
/**
|
||||
* Check if a matrix is symmetric.
|
||||
*
|
||||
* @param matrix Matrix to check.
|
||||
* @param raiseException If {@code true}, the method will throw an
|
||||
* exception if {@code matrix} is not symmetric.
|
||||
* @return {@code true} if {@code matrix} is symmetric.
|
||||
* @throws NonSymmetricMatrixException if the matrix is not symmetric and
|
||||
* {@code raiseException} is {@code true}.
|
||||
*/
|
||||
private boolean isSymmetric(final RealMatrix matrix,
|
||||
boolean raiseException) {
|
||||
final int rows = matrix.getRowDimension();
|
||||
final int columns = matrix.getColumnDimension();
|
||||
final double eps = 10 * rows * columns * MathUtils.EPSILON;
|
||||
for (int i = 0; i < rows; ++i) {
|
||||
for (int j = i + 1; j < columns; ++j) {
|
||||
final double mij = matrix.getEntry(i, j);
|
||||
final double mji = matrix.getEntry(j, i);
|
||||
if (FastMath.abs(mij - mji) >
|
||||
(FastMath.max(FastMath.abs(mij), FastMath.abs(mji)) * eps)) {
|
||||
if (raiseException) {
|
||||
throw new NonSymmetricMatrixException(i, j, eps);
|
||||
}
|
||||
return false;
|
||||
}
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the matrix V of the decomposition.
|
||||
|
|
@ -52,7 +179,19 @@ public interface EigenDecomposition {
|
|||
* or right-handed system).</p>
|
||||
* @return the V matrix
|
||||
*/
|
||||
RealMatrix getV();
|
||||
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;
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the block diagonal matrix D of the decomposition.
|
||||
|
|
@ -63,7 +202,13 @@ public interface EigenDecomposition {
|
|||
* @see #getRealEigenvalues()
|
||||
* @see #getImagEigenvalues()
|
||||
*/
|
||||
RealMatrix getD();
|
||||
public RealMatrix getD() {
|
||||
if (cachedD == null) {
|
||||
// cache the matrix for subsequent calls
|
||||
cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues);
|
||||
}
|
||||
return cachedD;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the transpose of the matrix V of the decomposition.
|
||||
|
|
@ -74,7 +219,20 @@ public interface EigenDecomposition {
|
|||
* or right-handed system).</p>
|
||||
* @return the transpose of the V matrix
|
||||
*/
|
||||
RealMatrix getVT();
|
||||
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 a copy of the real parts of the eigenvalues of the original matrix.
|
||||
|
|
@ -83,7 +241,9 @@ public interface EigenDecomposition {
|
|||
* @see #getRealEigenvalue(int)
|
||||
* @see #getImagEigenvalues()
|
||||
*/
|
||||
double[] getRealEigenvalues();
|
||||
public double[] getRealEigenvalues() {
|
||||
return realEigenvalues.clone();
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the real part of the i<sup>th</sup> eigenvalue of the original matrix.
|
||||
|
|
@ -93,7 +253,9 @@ public interface EigenDecomposition {
|
|||
* @see #getRealEigenvalues()
|
||||
* @see #getImagEigenvalue(int)
|
||||
*/
|
||||
double getRealEigenvalue(int i);
|
||||
public double getRealEigenvalue(final int i) {
|
||||
return realEigenvalues[i];
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns a copy of the imaginary parts of the eigenvalues of the original matrix.
|
||||
|
|
@ -102,7 +264,9 @@ public interface EigenDecomposition {
|
|||
* @see #getImagEigenvalue(int)
|
||||
* @see #getRealEigenvalues()
|
||||
*/
|
||||
double[] getImagEigenvalues();
|
||||
public double[] getImagEigenvalues() {
|
||||
return imagEigenvalues.clone();
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the imaginary part of the i<sup>th</sup> eigenvalue of the original matrix.
|
||||
|
|
@ -112,7 +276,9 @@ public interface EigenDecomposition {
|
|||
* @see #getImagEigenvalues()
|
||||
* @see #getRealEigenvalue(int)
|
||||
*/
|
||||
double getImagEigenvalue(int i);
|
||||
public double getImagEigenvalue(final int i) {
|
||||
return imagEigenvalues[i];
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns a copy of the i<sup>th</sup> eigenvector of the original matrix.
|
||||
|
|
@ -120,18 +286,338 @@ public interface EigenDecomposition {
|
|||
* @return copy of the i<sup>th</sup> eigenvector of the original matrix
|
||||
* @see #getD()
|
||||
*/
|
||||
RealVector getEigenvector(int i);
|
||||
public RealVector getEigenvector(final int i) {
|
||||
return eigenvectors[i].copy();
|
||||
}
|
||||
|
||||
/**
|
||||
* Return the determinant of the matrix
|
||||
* @return determinant of the matrix
|
||||
*/
|
||||
double getDeterminant();
|
||||
public double getDeterminant() {
|
||||
double determinant = 1;
|
||||
for (double lambda : realEigenvalues) {
|
||||
determinant *= lambda;
|
||||
}
|
||||
return determinant;
|
||||
}
|
||||
|
||||
/**
|
||||
* Get a solver for finding the A × X = B solution in exact linear sense.
|
||||
* @return a solver
|
||||
*/
|
||||
DecompositionSolver getSolver();
|
||||
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.
|
||||
* <p>
|
||||
* This method only find 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);
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
* 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 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);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* 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 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 = 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])<MathUtils.EPSILON*maxAbsoluteValue) {
|
||||
realEigenvalues[i]=0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
eigenvectors = new ArrayRealVector[n];
|
||||
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);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,552 +0,0 @@
|
|||
/*
|
||||
* 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.exception.MaxCountExceededException;
|
||||
import org.apache.commons.math.exception.DimensionMismatchException;
|
||||
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 <strong>symmetric</strong>
|
||||
* 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>
|
||||
* As of 2.0, this class supports only <strong>symmetric</strong> 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.
|
||||
* </p>
|
||||
* <p>
|
||||
* 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.
|
||||
* </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>
|
||||
* @version $Id$
|
||||
* @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 Matrix to decompose. It <em>must</em> be symmetric.
|
||||
* @param splitTolerance Dummy parameter (present for backward
|
||||
* compatibility only).
|
||||
* @throws NonSymmetricMatrixException if the matrix is not symmetric.
|
||||
* @throws MaxCountExceededException if the algorithm fails to converge.
|
||||
*/
|
||||
public EigenDecompositionImpl(final RealMatrix matrix,
|
||||
final double splitTolerance) {
|
||||
if (isSymmetric(matrix, true)) {
|
||||
transformToTridiagonal(matrix);
|
||||
findEigenVectors(transformer.getQ().getData());
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* 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).
|
||||
* @throws MaxCountExceededException if the algorithm fails to converge.
|
||||
*/
|
||||
public EigenDecompositionImpl(final double[] main,final double[] secondary,
|
||||
final double splitTolerance) {
|
||||
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<size;i++) {
|
||||
z[i][i]=1.0;
|
||||
}
|
||||
findEigenVectors(z);
|
||||
}
|
||||
|
||||
/**
|
||||
* Check if a matrix is symmetric.
|
||||
*
|
||||
* @param matrix Matrix to check.
|
||||
* @param raiseException If {@code true}, the method will throw an
|
||||
* exception if {@code matrix} is not symmetric.
|
||||
* @return {@code true} if {@code matrix} is symmetric.
|
||||
* @throws NonSymmetricMatrixException if the matrix is not symmetric and
|
||||
* {@code raiseException} is {@code true}.
|
||||
*/
|
||||
private boolean isSymmetric(final RealMatrix matrix,
|
||||
boolean raiseException) {
|
||||
final int rows = matrix.getRowDimension();
|
||||
final int columns = matrix.getColumnDimension();
|
||||
final double eps = 10 * rows * columns * MathUtils.EPSILON;
|
||||
for (int i = 0; i < rows; ++i) {
|
||||
for (int j = i + 1; j < columns; ++j) {
|
||||
final double mij = matrix.getEntry(i, j);
|
||||
final double mji = matrix.getEntry(j, i);
|
||||
if (FastMath.abs(mij - mji) >
|
||||
(FastMath.max(FastMath.abs(mij), FastMath.abs(mji)) * eps)) {
|
||||
if (raiseException) {
|
||||
throw new NonSymmetricMatrixException(i, j, eps);
|
||||
}
|
||||
return false;
|
||||
}
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
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;
|
||||
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix getD() {
|
||||
if (cachedD == null) {
|
||||
// cache the matrix for subsequent calls
|
||||
cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues);
|
||||
}
|
||||
return cachedD;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
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;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double[] getRealEigenvalues() {
|
||||
return realEigenvalues.clone();
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double getRealEigenvalue(final int i) {
|
||||
return realEigenvalues[i];
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double[] getImagEigenvalues() {
|
||||
return imagEigenvalues.clone();
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double getImagEigenvalue(final int i) {
|
||||
return imagEigenvalues[i];
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealVector getEigenvector(final int i) {
|
||||
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.
|
||||
* <p>
|
||||
* This method only find 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);
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
* 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 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);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* 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 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 = 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])<MathUtils.EPSILON*maxAbsoluteValue) {
|
||||
realEigenvalues[i]=0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
eigenvectors = new ArrayRealVector[n];
|
||||
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);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -28,7 +28,7 @@ import org.apache.commons.math.exception.NotPositiveException;
|
|||
import org.apache.commons.math.exception.OutOfRangeException;
|
||||
import org.apache.commons.math.exception.TooManyEvaluationsException;
|
||||
import org.apache.commons.math.linear.Array2DRowRealMatrix;
|
||||
import org.apache.commons.math.linear.EigenDecompositionImpl;
|
||||
import org.apache.commons.math.linear.EigenDecomposition;
|
||||
import org.apache.commons.math.linear.MatrixUtils;
|
||||
import org.apache.commons.math.linear.RealMatrix;
|
||||
import org.apache.commons.math.optimization.GoalType;
|
||||
|
|
@ -768,7 +768,7 @@ public class CMAESOptimizer extends
|
|||
// to achieve O(N^2)
|
||||
C = triu(C, 0).add(triu(C, 1).transpose());
|
||||
// enforce symmetry to prevent complex numbers
|
||||
EigenDecompositionImpl eig = new EigenDecompositionImpl(C, 1.0);
|
||||
EigenDecomposition eig = new EigenDecomposition(C, 1.0);
|
||||
B = eig.getV(); // eigen decomposition, B==normalized eigenvectors
|
||||
D = eig.getD();
|
||||
diagD = diag(D);
|
||||
|
|
|
|||
|
|
@ -37,8 +37,8 @@ public class EigenDecompositionTest {
|
|||
public void testDimension1() {
|
||||
RealMatrix matrix =
|
||||
MatrixUtils.createRealMatrix(new double[][] { { 1.5 } });
|
||||
EigenDecompositionImpl ed;
|
||||
ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
|
||||
EigenDecomposition ed;
|
||||
ed = new EigenDecomposition(matrix, MathUtils.SAFE_MIN);
|
||||
Assert.assertEquals(1.5, ed.getRealEigenvalue(0), 1.0e-15);
|
||||
}
|
||||
|
||||
|
|
@ -49,8 +49,8 @@ public class EigenDecompositionTest {
|
|||
{ 59.0, 12.0 },
|
||||
{ 12.0, 66.0 }
|
||||
});
|
||||
EigenDecompositionImpl ed;
|
||||
ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
|
||||
EigenDecomposition ed;
|
||||
ed = new EigenDecomposition(matrix, MathUtils.SAFE_MIN);
|
||||
Assert.assertEquals(75.0, ed.getRealEigenvalue(0), 1.0e-15);
|
||||
Assert.assertEquals(50.0, ed.getRealEigenvalue(1), 1.0e-15);
|
||||
}
|
||||
|
|
@ -63,8 +63,8 @@ public class EigenDecompositionTest {
|
|||
{ -4824.0, 8693.0, 7920.0 },
|
||||
{ -16560.0, 7920.0, 17300.0 }
|
||||
});
|
||||
EigenDecompositionImpl ed;
|
||||
ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
|
||||
EigenDecomposition ed;
|
||||
ed = new EigenDecomposition(matrix, MathUtils.SAFE_MIN);
|
||||
Assert.assertEquals(50000.0, ed.getRealEigenvalue(0), 3.0e-11);
|
||||
Assert.assertEquals(12500.0, ed.getRealEigenvalue(1), 3.0e-11);
|
||||
Assert.assertEquals( 3125.0, ed.getRealEigenvalue(2), 3.0e-11);
|
||||
|
|
@ -78,8 +78,8 @@ public class EigenDecompositionTest {
|
|||
{ 10, 20, 30 },
|
||||
{ 15, 30, 45 }
|
||||
});
|
||||
EigenDecompositionImpl ed;
|
||||
ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
|
||||
EigenDecomposition ed;
|
||||
ed = new EigenDecomposition(matrix, MathUtils.SAFE_MIN);
|
||||
Assert.assertEquals(70.0, ed.getRealEigenvalue(0), 3.0e-11);
|
||||
Assert.assertEquals(0.0, ed.getRealEigenvalue(1), 3.0e-11);
|
||||
Assert.assertEquals(0.0, ed.getRealEigenvalue(2), 3.0e-11);
|
||||
|
|
@ -94,8 +94,8 @@ public class EigenDecompositionTest {
|
|||
{ 0.000, 0.000, 0.164, -0.048 },
|
||||
{ 0.000, 0.000, -0.048, 0.136 }
|
||||
});
|
||||
EigenDecompositionImpl ed;
|
||||
ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
|
||||
EigenDecomposition ed;
|
||||
ed = new EigenDecomposition(matrix, MathUtils.SAFE_MIN);
|
||||
Assert.assertEquals(1.0, ed.getRealEigenvalue(0), 1.0e-15);
|
||||
Assert.assertEquals(0.4, ed.getRealEigenvalue(1), 1.0e-15);
|
||||
Assert.assertEquals(0.2, ed.getRealEigenvalue(2), 1.0e-15);
|
||||
|
|
@ -111,8 +111,8 @@ public class EigenDecompositionTest {
|
|||
{ 0.1152, -0.2304, 0.3088, -0.1344 },
|
||||
{ -0.2976, 0.1152, -0.1344, 0.3872 }
|
||||
});
|
||||
EigenDecompositionImpl ed;
|
||||
ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
|
||||
EigenDecomposition ed;
|
||||
ed = new EigenDecomposition(matrix, MathUtils.SAFE_MIN);
|
||||
Assert.assertEquals(1.0, ed.getRealEigenvalue(0), 1.0e-15);
|
||||
Assert.assertEquals(0.4, ed.getRealEigenvalue(1), 1.0e-15);
|
||||
Assert.assertEquals(0.2, ed.getRealEigenvalue(2), 1.0e-15);
|
||||
|
|
@ -143,8 +143,8 @@ public class EigenDecompositionTest {
|
|||
new ArrayRealVector(new double[] { -0.584677060845929, 0.367177264979103, 0.721453187784497, -0.052971054621812, 0.005740715188257 })
|
||||
};
|
||||
|
||||
EigenDecompositionImpl decomposition;
|
||||
decomposition = new EigenDecompositionImpl(mainTridiagonal,
|
||||
EigenDecomposition decomposition;
|
||||
decomposition = new EigenDecomposition(mainTridiagonal,
|
||||
secondaryTridiagonal,
|
||||
MathUtils.SAFE_MIN);
|
||||
|
||||
|
|
@ -188,8 +188,8 @@ public class EigenDecompositionTest {
|
|||
};
|
||||
|
||||
// the following line triggers the exception
|
||||
EigenDecompositionImpl decomposition;
|
||||
decomposition = new EigenDecompositionImpl(mainTridiagonal,
|
||||
EigenDecomposition decomposition;
|
||||
decomposition = new EigenDecomposition(mainTridiagonal,
|
||||
secondaryTridiagonal,
|
||||
MathUtils.SAFE_MIN);
|
||||
|
||||
|
|
@ -235,8 +235,8 @@ public class EigenDecompositionTest {
|
|||
};
|
||||
|
||||
// the following line triggers the exception
|
||||
EigenDecompositionImpl decomposition;
|
||||
decomposition = new EigenDecompositionImpl(mainTridiagonal,
|
||||
EigenDecomposition decomposition;
|
||||
decomposition = new EigenDecomposition(mainTridiagonal,
|
||||
secondaryTridiagonal,
|
||||
MathUtils.SAFE_MIN);
|
||||
|
||||
|
|
@ -267,8 +267,8 @@ public class EigenDecompositionTest {
|
|||
Arrays.sort(ref);
|
||||
TriDiagonalTransformer t =
|
||||
new TriDiagonalTransformer(createTestMatrix(r, ref));
|
||||
EigenDecompositionImpl ed;
|
||||
ed = new EigenDecompositionImpl(t.getMainDiagonalRef(),
|
||||
EigenDecomposition ed;
|
||||
ed = new EigenDecomposition(t.getMainDiagonalRef(),
|
||||
t.getSecondaryDiagonalRef(),
|
||||
MathUtils.SAFE_MIN);
|
||||
double[] eigenValues = ed.getRealEigenvalues();
|
||||
|
|
@ -283,8 +283,8 @@ public class EigenDecompositionTest {
|
|||
@Test
|
||||
public void testDimensions() {
|
||||
final int m = matrix.getRowDimension();
|
||||
EigenDecompositionImpl ed;
|
||||
ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
|
||||
EigenDecomposition ed;
|
||||
ed = new EigenDecomposition(matrix, MathUtils.SAFE_MIN);
|
||||
Assert.assertEquals(m, ed.getV().getRowDimension());
|
||||
Assert.assertEquals(m, ed.getV().getColumnDimension());
|
||||
Assert.assertEquals(m, ed.getD().getColumnDimension());
|
||||
|
|
@ -296,8 +296,8 @@ public class EigenDecompositionTest {
|
|||
/** test eigenvalues */
|
||||
@Test
|
||||
public void testEigenvalues() {
|
||||
EigenDecompositionImpl ed;
|
||||
ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
|
||||
EigenDecomposition ed;
|
||||
ed = new EigenDecomposition(matrix, MathUtils.SAFE_MIN);
|
||||
double[] eigenValues = ed.getRealEigenvalues();
|
||||
Assert.assertEquals(refValues.length, eigenValues.length);
|
||||
for (int i = 0; i < refValues.length; ++i) {
|
||||
|
|
@ -314,8 +314,8 @@ public class EigenDecompositionTest {
|
|||
bigValues[i] = 2 * r.nextDouble() - 1;
|
||||
}
|
||||
Arrays.sort(bigValues);
|
||||
EigenDecompositionImpl ed;
|
||||
ed = new EigenDecompositionImpl(createTestMatrix(r, bigValues),
|
||||
EigenDecomposition ed;
|
||||
ed = new EigenDecomposition(createTestMatrix(r, bigValues),
|
||||
MathUtils.SAFE_MIN);
|
||||
double[] eigenValues = ed.getRealEigenvalues();
|
||||
Assert.assertEquals(bigValues.length, eigenValues.length);
|
||||
|
|
@ -327,8 +327,8 @@ public class EigenDecompositionTest {
|
|||
/** test eigenvectors */
|
||||
@Test
|
||||
public void testEigenvectors() {
|
||||
EigenDecompositionImpl ed;
|
||||
ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
|
||||
EigenDecomposition ed;
|
||||
ed = new EigenDecomposition(matrix, MathUtils.SAFE_MIN);
|
||||
for (int i = 0; i < matrix.getRowDimension(); ++i) {
|
||||
double lambda = ed.getRealEigenvalue(i);
|
||||
RealVector v = ed.getEigenvector(i);
|
||||
|
|
@ -340,8 +340,8 @@ public class EigenDecompositionTest {
|
|||
/** test A = VDVt */
|
||||
@Test
|
||||
public void testAEqualVDVt() {
|
||||
EigenDecompositionImpl ed;
|
||||
ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
|
||||
EigenDecomposition ed;
|
||||
ed = new EigenDecomposition(matrix, MathUtils.SAFE_MIN);
|
||||
RealMatrix v = ed.getV();
|
||||
RealMatrix d = ed.getD();
|
||||
RealMatrix vT = ed.getVT();
|
||||
|
|
@ -352,7 +352,7 @@ public class EigenDecompositionTest {
|
|||
/** test that V is orthogonal */
|
||||
@Test
|
||||
public void testVOrthogonal() {
|
||||
RealMatrix v = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN).getV();
|
||||
RealMatrix v = new EigenDecomposition(matrix, MathUtils.SAFE_MIN).getV();
|
||||
RealMatrix vTv = v.transpose().multiply(v);
|
||||
RealMatrix id = MatrixUtils.createRealIdentityMatrix(vTv.getRowDimension());
|
||||
Assert.assertEquals(0, vTv.subtract(id).getNorm(), 2.0e-13);
|
||||
|
|
@ -363,8 +363,8 @@ public class EigenDecompositionTest {
|
|||
public void testDiagonal() {
|
||||
double[] diagonal = new double[] { -3.0, -2.0, 2.0, 5.0 };
|
||||
RealMatrix m = createDiagonalMatrix(diagonal, diagonal.length, diagonal.length);
|
||||
EigenDecompositionImpl ed;
|
||||
ed = new EigenDecompositionImpl(m, MathUtils.SAFE_MIN);
|
||||
EigenDecomposition ed;
|
||||
ed = new EigenDecomposition(m, MathUtils.SAFE_MIN);
|
||||
Assert.assertEquals(diagonal[0], ed.getRealEigenvalue(3), 2.0e-15);
|
||||
Assert.assertEquals(diagonal[1], ed.getRealEigenvalue(2), 2.0e-15);
|
||||
Assert.assertEquals(diagonal[2], ed.getRealEigenvalue(1), 2.0e-15);
|
||||
|
|
@ -381,8 +381,8 @@ public class EigenDecompositionTest {
|
|||
{2, 0, 2},
|
||||
{4, 2, 3}
|
||||
});
|
||||
EigenDecompositionImpl ed;
|
||||
ed = new EigenDecompositionImpl(repeated, MathUtils.SAFE_MIN);
|
||||
EigenDecomposition ed;
|
||||
ed = new EigenDecomposition(repeated, MathUtils.SAFE_MIN);
|
||||
checkEigenValues((new double[] {8, -1, -1}), ed, 1E-12);
|
||||
checkEigenVector((new double[] {2, 1, 2}), ed, 1E-12);
|
||||
}
|
||||
|
|
@ -397,8 +397,8 @@ public class EigenDecompositionTest {
|
|||
{1, 3, -4},
|
||||
{-4, -4, 8}
|
||||
});
|
||||
EigenDecompositionImpl ed;
|
||||
ed = new EigenDecompositionImpl(distinct, MathUtils.SAFE_MIN);
|
||||
EigenDecomposition ed;
|
||||
ed = new EigenDecomposition(distinct, MathUtils.SAFE_MIN);
|
||||
checkEigenValues((new double[] {2, 0, 12}), ed, 1E-12);
|
||||
checkEigenVector((new double[] {1, -1, 0}), ed, 1E-12);
|
||||
checkEigenVector((new double[] {1, 1, 1}), ed, 1E-12);
|
||||
|
|
@ -415,8 +415,8 @@ public class EigenDecompositionTest {
|
|||
{ 1.0, 1.0, 0.0 },
|
||||
{ -1.0,0.0, 1.0 }
|
||||
});
|
||||
EigenDecompositionImpl ed;
|
||||
ed = new EigenDecompositionImpl(indefinite, MathUtils.SAFE_MIN);
|
||||
EigenDecomposition ed;
|
||||
ed = new EigenDecomposition(indefinite, MathUtils.SAFE_MIN);
|
||||
checkEigenValues((new double[] {2, 1, -1}), ed, 1E-12);
|
||||
double isqrt3 = 1/FastMath.sqrt(3.0);
|
||||
checkEigenVector((new double[] {isqrt3,isqrt3,-isqrt3}), ed, 1E-12);
|
||||
|
|
@ -431,7 +431,7 @@ public class EigenDecompositionTest {
|
|||
* values to differ by tolerance.
|
||||
*/
|
||||
protected void checkEigenValues(double[] targetValues,
|
||||
EigenDecompositionImpl ed, double tolerance) {
|
||||
EigenDecomposition ed, double tolerance) {
|
||||
double[] observed = ed.getRealEigenvalues();
|
||||
for (int i = 0; i < observed.length; i++) {
|
||||
Assert.assertTrue(isIncludedValue(observed[i], targetValues, tolerance));
|
||||
|
|
@ -463,7 +463,7 @@ public class EigenDecompositionTest {
|
|||
* used to find vectors in one-dimensional eigenspaces.
|
||||
*/
|
||||
protected void checkEigenVector(double[] eigenVector,
|
||||
EigenDecompositionImpl ed, double tolerance) {
|
||||
EigenDecomposition ed, double tolerance) {
|
||||
Assert.assertTrue(isIncludedColumn(eigenVector, ed.getV(), tolerance));
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -33,7 +33,7 @@ public class EigenSolverTest {
|
|||
Random r = new Random(9994100315209l);
|
||||
RealMatrix m =
|
||||
EigenDecompositionTest.createTestMatrix(r, new double[] { 1.0, 0.0, -1.0, -2.0, -3.0 });
|
||||
DecompositionSolver es = new EigenDecompositionImpl(m, MathUtils.SAFE_MIN).getSolver();
|
||||
DecompositionSolver es = new EigenDecomposition(m, MathUtils.SAFE_MIN).getSolver();
|
||||
Assert.assertFalse(es.isNonSingular());
|
||||
try {
|
||||
es.getInverse();
|
||||
|
|
@ -49,7 +49,7 @@ public class EigenSolverTest {
|
|||
Random r = new Random(9994100315209l);
|
||||
RealMatrix m =
|
||||
EigenDecompositionTest.createTestMatrix(r, new double[] { 1.0, 0.5, -1.0, -2.0, -3.0 });
|
||||
DecompositionSolver es = new EigenDecompositionImpl(m, MathUtils.SAFE_MIN).getSolver();
|
||||
DecompositionSolver es = new EigenDecomposition(m, MathUtils.SAFE_MIN).getSolver();
|
||||
Assert.assertTrue(es.isNonSingular());
|
||||
RealMatrix inverse = es.getInverse();
|
||||
RealMatrix error =
|
||||
|
|
@ -65,7 +65,7 @@ public class EigenSolverTest {
|
|||
};
|
||||
final RealMatrix matrix = EigenDecompositionTest.createTestMatrix(new Random(35992629946426l), refValues);
|
||||
|
||||
DecompositionSolver es = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN).getSolver();
|
||||
DecompositionSolver es = new EigenDecomposition(matrix, MathUtils.SAFE_MIN).getSolver();
|
||||
RealMatrix b = MatrixUtils.createRealMatrix(new double[2][2]);
|
||||
try {
|
||||
es.solve(b);
|
||||
|
|
@ -98,7 +98,7 @@ public class EigenSolverTest {
|
|||
{ 40, 2, 21, 9, 51, 19 },
|
||||
{ 14, -1, 8, 0, 19, 14 }
|
||||
});
|
||||
DecompositionSolver es = new EigenDecompositionImpl(m, MathUtils.SAFE_MIN).getSolver();
|
||||
DecompositionSolver es = new EigenDecomposition(m, MathUtils.SAFE_MIN).getSolver();
|
||||
RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
|
||||
{ 1561, 269, 188 },
|
||||
{ 69, -21, 70 },
|
||||
|
|
|
|||
Loading…
Reference in New Issue