mirror of
https://github.com/apache/commons-math.git
synced 2025-02-06 10:09:26 +00:00
prepared interface for general eigendecomposition
we currently support only symmetric matrices and real eigenvalues Jira issue MATH-235 asks for asymmetric matrices support, which needs an API adaptation. Preparing the interface yet avoids introducing incompatibility later. The current status is explained in javadoc and matrix symmetry is checked at runtime with an InvalidMatrixException triggered if asymmetry is detected. git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@728915 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
6916b0314d
commit
25a6b32fe1
@ -21,7 +21,7 @@ import java.io.Serializable;
|
||||
|
||||
/**
|
||||
* An interface to classes that implement an algorithm to calculate the
|
||||
* eigen decomposition of a real symmetric matrix.
|
||||
* eigen decomposition of a real matrix.
|
||||
* <p>The eigen decomposition of matrix A is a set of two matrices:
|
||||
* V and D such that A = V × D × V<sup>T</sup>.
|
||||
* A, V and D are all m × m matrices.</p>
|
||||
@ -30,13 +30,10 @@ import java.io.Serializable;
|
||||
* library, with the following changes:</p>
|
||||
* <ul>
|
||||
* <li>a {@link #getVT() getVt} method has been added,</li>
|
||||
* <li>a {@link #getEigenvalue(int) getEigenvalue} method to pick up a single
|
||||
* eigenvalue has been added,</li>
|
||||
* <li>two {@link #getRealEigenvalue(int) getRealEigenvalue} and {@link #getImagEigenvalue(int)
|
||||
* getImagEigenvalue} methods to pick up a single eigenvalue have been added,</li>
|
||||
* <li>a {@link #getEigenvector(int) getEigenvector} method to pick up a single
|
||||
* eigenvector has been added,</li>
|
||||
* <li>the <code>getRealEigenvalues</code> method has been renamed as {@link
|
||||
* #getEigenValues() getEigenValues},</li>
|
||||
* <li>the <code>getImagEigenvalues</code> method has been removed</li>
|
||||
* <li>a {@link #getDeterminant() getDeterminant} method has been added.</li>
|
||||
* <li>a {@link #getSolver() getSolver} method has been added.</li>
|
||||
* </ul>
|
||||
@ -56,9 +53,10 @@ public interface EigenDecomposition extends Serializable {
|
||||
RealMatrix getV();
|
||||
|
||||
/**
|
||||
* Returns the diagonal matrix D of the decomposition.
|
||||
* <p>D is a diagonal matrix.</p>
|
||||
* <p>The values on the diagonal are the eigenvalues of the original matrix.</p>
|
||||
* Returns the block diagonal matrix D of the decomposition.
|
||||
* <p>D is a block diagonal matrix.</p>
|
||||
* <p>Real eigenvalues are on the diagonal while complex values are on
|
||||
* 2x2 blocks { {real +imaginary}, {-imaginary, real} }.</p>
|
||||
* @return the D matrix
|
||||
* @see #getEigenValues()
|
||||
*/
|
||||
@ -73,19 +71,42 @@ public interface EigenDecomposition extends Serializable {
|
||||
RealMatrix getVT();
|
||||
|
||||
/**
|
||||
* Returns a copy of the eigenvalues of the original matrix.
|
||||
* @return a copy of the eigenvalues of the original matrix
|
||||
* Returns 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()
|
||||
*/
|
||||
double[] getEigenvalues();
|
||||
double[] getRealEigenvalues();
|
||||
|
||||
/**
|
||||
* Returns the i<sup>th</sup> eigenvalue of the original matrix.
|
||||
* Returns the real part of the i<sup>th</sup> eigenvalue of the original matrix.
|
||||
* @param i index of the eigenvalue (counting from 0)
|
||||
* @return i<sup>th</sup> eigenvalue of the original matrix
|
||||
* @return real part of the i<sup>th</sup> eigenvalue of the original matrix
|
||||
* @see #getD()
|
||||
* @see #getRealEigenvalues()
|
||||
* @see #getImagEigenvalue(int)
|
||||
*/
|
||||
double getEigenvalue(int i);
|
||||
double getRealEigenvalue(int i);
|
||||
|
||||
/**
|
||||
* Returns 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()
|
||||
*/
|
||||
double[] getImagEigenvalues();
|
||||
|
||||
/**
|
||||
* Returns the imaginary part of the i<sup>th</sup> eigenvalue of the original matrix.
|
||||
* @param i index of the eigenvalue (counting from 0)
|
||||
* @return imaginary part of the i<sup>th</sup> eigenvalue of the original matrix
|
||||
* @see #getD()
|
||||
* @see #getImagEigenvalues()
|
||||
* @see #getRealEigenvalue(int)
|
||||
*/
|
||||
double getImagEigenvalue(int i);
|
||||
|
||||
/**
|
||||
* Returns a copy of the i<sup>th</sup> eigenvector of the original matrix.
|
||||
@ -102,7 +123,7 @@ public interface EigenDecomposition extends Serializable {
|
||||
double getDeterminant();
|
||||
|
||||
/**
|
||||
* Get a solver for A × X = B.
|
||||
* Get a solver for finding the A × X = B solution in exact linear sense.
|
||||
* @return a solver
|
||||
*/
|
||||
DecompositionSolver getSolver();
|
||||
|
||||
@ -27,9 +27,13 @@ import org.apache.commons.math.util.MathUtils;
|
||||
|
||||
/**
|
||||
* Calculates the eigen decomposition of a <strong>symmetric</strong> matrix.
|
||||
* <p>The eigen decomposition of symmetric matrix A is a set of two matrices:
|
||||
* <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>Eigenvalues are computed as soon as the matrix is decomposed, but eigenvectors
|
||||
@ -132,8 +136,11 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
/** Shift ratio with respect to dMin used when tType == 6. */
|
||||
private double g;
|
||||
|
||||
/** Eigenvalues. */
|
||||
private double[] eigenvalues;
|
||||
/** Real part of the realEigenvalues. */
|
||||
private double[] realEigenvalues;
|
||||
|
||||
/** Imaginary part of the realEigenvalues. */
|
||||
private double[] imagEigenvalues;
|
||||
|
||||
/** Eigenvectors. */
|
||||
private RealVectorImpl[] eigenvectors;
|
||||
@ -163,9 +170,16 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
public EigenDecompositionImpl(final RealMatrix matrix,
|
||||
final double splitTolerance)
|
||||
throws InvalidMatrixException {
|
||||
this.splitTolerance = splitTolerance;
|
||||
transformToTridiagonal(matrix);
|
||||
decompose();
|
||||
if (isSymmetric(matrix)) {
|
||||
this.splitTolerance = splitTolerance;
|
||||
transformToTridiagonal(matrix);
|
||||
decompose();
|
||||
} 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("eigen decomposition of assymetric matrices not supported yet",
|
||||
null);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
@ -200,6 +214,27 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
* Check if a matrix is symmetric.
|
||||
* @param matrix matrix to check
|
||||
* @return true if matrix is symmetric
|
||||
*/
|
||||
private boolean isSymmetric(final RealMatrix matrix) {
|
||||
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 (Math.abs(mij - mji) > (Math.max(Math.abs(mij), Math.abs(mji)) * eps)) {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
/**
|
||||
* Decompose a tridiagonal symmetric matrix.
|
||||
* @exception InvalidMatrixException (wrapping a {@link ConvergenceException}
|
||||
@ -215,7 +250,7 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
// compute the Gershgorin circles
|
||||
computeGershgorinCircles();
|
||||
|
||||
// find all the eigenvalues
|
||||
// find all the realEigenvalues
|
||||
findEigenvalues();
|
||||
|
||||
// we will search for eigenvectors only if required
|
||||
@ -251,7 +286,7 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
throws InvalidMatrixException {
|
||||
if (cachedD == null) {
|
||||
// cache the matrix for subsequent calls
|
||||
cachedD = MatrixUtils.createRealDiagonalMatrix(eigenvalues);
|
||||
cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues);
|
||||
}
|
||||
return cachedD;
|
||||
}
|
||||
@ -280,15 +315,27 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double[] getEigenvalues()
|
||||
public double[] getRealEigenvalues()
|
||||
throws InvalidMatrixException {
|
||||
return eigenvalues.clone();
|
||||
return realEigenvalues.clone();
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double getEigenvalue(final int i)
|
||||
public double getRealEigenvalue(final int i)
|
||||
throws InvalidMatrixException, ArrayIndexOutOfBoundsException {
|
||||
return eigenvalues[i];
|
||||
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} */
|
||||
@ -307,7 +354,7 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
*/
|
||||
public double getDeterminant() {
|
||||
double determinant = 1;
|
||||
for (double lambda : eigenvalues) {
|
||||
for (double lambda : realEigenvalues) {
|
||||
determinant *= lambda;
|
||||
}
|
||||
return determinant;
|
||||
@ -318,7 +365,7 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
if (eigenvectors == null) {
|
||||
findEigenVectors();
|
||||
}
|
||||
return new Solver(eigenvalues, eigenvectors);
|
||||
return new Solver(realEigenvalues, eigenvectors);
|
||||
}
|
||||
|
||||
/** Specialized solver. */
|
||||
@ -335,7 +382,7 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
|
||||
/**
|
||||
* Build a solver from decomposed matrix.
|
||||
* @param eigenvalues eigenvalues
|
||||
* @param realEigenvalues realEigenvalues
|
||||
* @param eigenvectors eigenvectors
|
||||
*/
|
||||
private Solver(final double[] eigenvalues, final RealVectorImpl[] eigenvectors) {
|
||||
@ -555,7 +602,7 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
}
|
||||
|
||||
/**
|
||||
* Find the eigenvalues.
|
||||
* Find the realEigenvalues.
|
||||
* @exception InvalidMatrixException if a block cannot be diagonalized
|
||||
*/
|
||||
private void findEigenvalues()
|
||||
@ -564,8 +611,9 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
// compute splitting points
|
||||
List<Integer> splitIndices = computeSplits();
|
||||
|
||||
// find eigenvalues in each block
|
||||
eigenvalues = new double[main.length];
|
||||
// find realEigenvalues in each block
|
||||
realEigenvalues = new double[main.length];
|
||||
imagEigenvalues = new double[main.length];
|
||||
int begin = 0;
|
||||
for (final int end : splitIndices) {
|
||||
final int n = end - begin;
|
||||
@ -605,14 +653,14 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
// apply general dqd/dqds method
|
||||
processGeneralBlock(n);
|
||||
|
||||
// extract eigenvalues
|
||||
// extract realEigenvalues
|
||||
if (chooseLeft) {
|
||||
for (int i = 0; i < n; ++i) {
|
||||
eigenvalues[begin + i] = lambda + work[4 * i];
|
||||
realEigenvalues[begin + i] = lambda + work[4 * i];
|
||||
}
|
||||
} else {
|
||||
for (int i = 0; i < n; ++i) {
|
||||
eigenvalues[begin + i] = lambda - work[4 * i];
|
||||
realEigenvalues[begin + i] = lambda - work[4 * i];
|
||||
}
|
||||
}
|
||||
|
||||
@ -620,12 +668,12 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
begin = end;
|
||||
}
|
||||
|
||||
// sort the eigenvalues in decreasing order
|
||||
Arrays.sort(eigenvalues);
|
||||
for (int i = 0, j = eigenvalues.length - 1; i < j; ++i, --j) {
|
||||
final double tmp = eigenvalues[i];
|
||||
eigenvalues[i] = eigenvalues[j];
|
||||
eigenvalues[j] = tmp;
|
||||
// sort the realEigenvalues in decreasing order
|
||||
Arrays.sort(realEigenvalues);
|
||||
for (int i = 0, j = realEigenvalues.length - 1; i < j; ++i, --j) {
|
||||
final double tmp = realEigenvalues[i];
|
||||
realEigenvalues[i] = realEigenvalues[j];
|
||||
realEigenvalues[j] = tmp;
|
||||
}
|
||||
|
||||
}
|
||||
@ -662,11 +710,11 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
* @param index index of the first row of the block
|
||||
*/
|
||||
private void process1RowBlock(final int index) {
|
||||
eigenvalues[index] = main[index];
|
||||
realEigenvalues[index] = main[index];
|
||||
}
|
||||
|
||||
/**
|
||||
* Find eigenvalues in a block with 2 rows.
|
||||
* Find realEigenvalues in a block with 2 rows.
|
||||
* <p>In low dimensions, we simply solve the characteristic polynomial.</p>
|
||||
* @param index index of the first row of the block
|
||||
* @exception InvalidMatrixException if characteristic polynomial cannot be solved
|
||||
@ -688,13 +736,13 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
}
|
||||
|
||||
final double largestRoot = 0.5 * (s + Math.sqrt(delta));
|
||||
eigenvalues[index] = largestRoot;
|
||||
eigenvalues[index + 1] = p / largestRoot;
|
||||
realEigenvalues[index] = largestRoot;
|
||||
realEigenvalues[index + 1] = p / largestRoot;
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
* Find eigenvalues in a block with 3 rows.
|
||||
* Find realEigenvalues in a block with 3 rows.
|
||||
* <p>In low dimensions, we simply solve the characteristic polynomial.</p>
|
||||
* @param index index of the first row of the block
|
||||
* @exception InvalidMatrixException if diagonal elements are not positive
|
||||
@ -722,7 +770,7 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
final double delta = q * q * q + r * r;
|
||||
if (delta >= 0) {
|
||||
// in fact, there are solutions to the equation, but in the context
|
||||
// of symmetric eigenvalues problem, there should be three distinct
|
||||
// of symmetric realEigenvalues problem, there should be three distinct
|
||||
// real roots, so we throw an error if this condition is not met
|
||||
throw new InvalidMatrixException("cannot solve degree {0} equation", new Object[] { 3 });
|
||||
}
|
||||
@ -749,14 +797,14 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
z0 = z1;
|
||||
z1 = t;
|
||||
}
|
||||
eigenvalues[index] = z0;
|
||||
eigenvalues[index + 1] = z1;
|
||||
eigenvalues[index + 2] = z2;
|
||||
realEigenvalues[index] = z0;
|
||||
realEigenvalues[index + 1] = z1;
|
||||
realEigenvalues[index + 2] = z2;
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
* Find eigenvalues using dqd/dqds algorithms.
|
||||
* Find realEigenvalues using dqd/dqds algorithms.
|
||||
* <p>This implementation is based on Beresford N. Parlett
|
||||
* and Osni A. Marques paper <a
|
||||
* href="http://www.netlib.org/lapack/lawnspdf/lawn155.pdf">An
|
||||
@ -944,7 +992,7 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
|
||||
g = 0.0;
|
||||
|
||||
// step 1: accepting eigenvalues
|
||||
// step 1: accepting realEigenvalues
|
||||
int deflatedEnd = end;
|
||||
for (boolean deflating = true; deflating;) {
|
||||
|
||||
@ -968,7 +1016,7 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
(work[k - 9] <= TOLERANCE_2 * sigma) ||
|
||||
(work[k - 2 * pingPong - 8] <= TOLERANCE_2 * work[k - 11])) {
|
||||
|
||||
// two eigenvalues found, deflate array
|
||||
// two realEigenvalues found, deflate array
|
||||
if (work[k - 3] > work[k - 7]) {
|
||||
final double tmp = work[k - 3];
|
||||
work[k - 3] = work[k - 7];
|
||||
@ -992,7 +1040,7 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
deflatedEnd -= 2;
|
||||
} else {
|
||||
|
||||
// no more eigenvalues found, we need to iterate
|
||||
// no more realEigenvalues found, we need to iterate
|
||||
deflating = false;
|
||||
|
||||
}
|
||||
@ -1097,10 +1145,10 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute an interval containing all eigenvalues of a block.
|
||||
* Compute an interval containing all realEigenvalues of a block.
|
||||
* @param index index of the first row of the block
|
||||
* @param n number of rows of the block
|
||||
* @return an interval containing the eigenvalues
|
||||
* @return an interval containing the realEigenvalues
|
||||
*/
|
||||
private double[] eigenvaluesRange(final int index, final int n) {
|
||||
|
||||
@ -1171,11 +1219,11 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
}
|
||||
|
||||
/**
|
||||
* Count the number of eigenvalues below a point.
|
||||
* @param t value below which we must count the number of eigenvalues
|
||||
* Count the number of realEigenvalues below a point.
|
||||
* @param t value below which we must count the number of realEigenvalues
|
||||
* @param index index of the first row of the block
|
||||
* @param n number of rows of the block
|
||||
* @return number of eigenvalues smaller than t
|
||||
* @return number of realEigenvalues smaller than t
|
||||
*/
|
||||
private int countEigenValues(final double t, final int index, final int n) {
|
||||
double ratio = main[index] - t;
|
||||
@ -1376,7 +1424,7 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
* <p>This implementation is a translation of the LAPACK routine DLAZQ4.</p>
|
||||
* @param start start index
|
||||
* @param end end index
|
||||
* @param deflated number of eigenvalues just deflated
|
||||
* @param deflated number of realEigenvalues just deflated
|
||||
*/
|
||||
private void computeShiftIncrement(final int start, final int end, final int deflated) {
|
||||
|
||||
@ -1395,7 +1443,7 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
int nn = 4 * end + pingPong - 1;
|
||||
switch (deflated) {
|
||||
|
||||
case 0 : // no eigenvalues deflated.
|
||||
case 0 : // no realEigenvalues deflated.
|
||||
if (dMin == dN || dMin == dN1) {
|
||||
|
||||
double b1 = Math.sqrt(work[nn - 3]) * Math.sqrt(work[nn - 5]);
|
||||
@ -1577,7 +1625,7 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
}
|
||||
break;
|
||||
|
||||
case 2 : // two eigenvalues deflated. use dMin2, dN2 for dMin and dN.
|
||||
case 2 : // two realEigenvalues deflated. use dMin2, dN2 for dMin and dN.
|
||||
|
||||
// cases 10 and 11.
|
||||
if (dMin2 == dN2 && 2 * work[nn - 5] < work[nn - 7]) {
|
||||
@ -1615,7 +1663,7 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
}
|
||||
break;
|
||||
|
||||
default : // case 12, more than two eigenvalues deflated. no information.
|
||||
default : // case 12, more than two realEigenvalues deflated. no information.
|
||||
tau = 0.0;
|
||||
tType = -12;
|
||||
}
|
||||
@ -1665,7 +1713,7 @@ public class EigenDecompositionImpl implements EigenDecomposition {
|
||||
|
||||
// compute eigenvectors
|
||||
for (int i = 0; i < m; ++i) {
|
||||
eigenvectors[i] = findEigenvector(eigenvalues[i], d, l);
|
||||
eigenvectors[i] = findEigenvector(realEigenvalues[i], d, l);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
@ -45,46 +45,46 @@ public class EigenDecompositionImplTest extends TestCase {
|
||||
RealMatrix matrix =
|
||||
MatrixUtils.createRealMatrix(new double[][] { { 1.5 } });
|
||||
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
|
||||
assertEquals(1.5, ed.getEigenvalue(0), 1.0e-15);
|
||||
assertEquals(1.5, ed.getRealEigenvalue(0), 1.0e-15);
|
||||
}
|
||||
|
||||
public void testDimension2() {
|
||||
RealMatrix matrix =
|
||||
MatrixUtils.createRealMatrix(new double[][] {
|
||||
{ 59.0, 12.0 },
|
||||
{ Double.NaN, 66.0 }
|
||||
{ 59.0, 12.0 },
|
||||
{ 12.0, 66.0 }
|
||||
});
|
||||
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
|
||||
assertEquals(75.0, ed.getEigenvalue(0), 1.0e-15);
|
||||
assertEquals(50.0, ed.getEigenvalue(1), 1.0e-15);
|
||||
assertEquals(75.0, ed.getRealEigenvalue(0), 1.0e-15);
|
||||
assertEquals(50.0, ed.getRealEigenvalue(1), 1.0e-15);
|
||||
}
|
||||
|
||||
public void testDimension3() {
|
||||
RealMatrix matrix =
|
||||
MatrixUtils.createRealMatrix(new double[][] {
|
||||
{ 39632.0, -4824.0, -16560.0 },
|
||||
{ Double.NaN, 8693.0, 7920.0 },
|
||||
{ Double.NaN, Double.NaN, 17300.0 }
|
||||
{ 39632.0, -4824.0, -16560.0 },
|
||||
{ -4824.0, 8693.0, 7920.0 },
|
||||
{ -16560.0, 7920.0, 17300.0 }
|
||||
});
|
||||
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
|
||||
assertEquals(50000.0, ed.getEigenvalue(0), 3.0e-11);
|
||||
assertEquals(12500.0, ed.getEigenvalue(1), 3.0e-11);
|
||||
assertEquals( 3125.0, ed.getEigenvalue(2), 3.0e-11);
|
||||
assertEquals(50000.0, ed.getRealEigenvalue(0), 3.0e-11);
|
||||
assertEquals(12500.0, ed.getRealEigenvalue(1), 3.0e-11);
|
||||
assertEquals( 3125.0, ed.getRealEigenvalue(2), 3.0e-11);
|
||||
}
|
||||
|
||||
public void testDimension4WithSplit() {
|
||||
RealMatrix matrix =
|
||||
MatrixUtils.createRealMatrix(new double[][] {
|
||||
{ 0.784, -0.288, 0.000, 0.000 },
|
||||
{ Double.NaN, 0.616, 0.000, 0.000 },
|
||||
{ Double.NaN, Double.NaN, 0.164, -0.048 },
|
||||
{ Double.NaN, Double.NaN, Double.NaN, 0.136 }
|
||||
{ 0.784, -0.288, 0.000, 0.000 },
|
||||
{ -0.288, 0.616, 0.000, 0.000 },
|
||||
{ 0.000, 0.000, 0.164, -0.048 },
|
||||
{ 0.000, 0.000, -0.048, 0.136 }
|
||||
});
|
||||
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
|
||||
assertEquals(1.0, ed.getEigenvalue(0), 1.0e-15);
|
||||
assertEquals(0.4, ed.getEigenvalue(1), 1.0e-15);
|
||||
assertEquals(0.2, ed.getEigenvalue(2), 1.0e-15);
|
||||
assertEquals(0.1, ed.getEigenvalue(3), 1.0e-15);
|
||||
assertEquals(1.0, ed.getRealEigenvalue(0), 1.0e-15);
|
||||
assertEquals(0.4, ed.getRealEigenvalue(1), 1.0e-15);
|
||||
assertEquals(0.2, ed.getRealEigenvalue(2), 1.0e-15);
|
||||
assertEquals(0.1, ed.getRealEigenvalue(3), 1.0e-15);
|
||||
}
|
||||
|
||||
public void testDimension4WithoutSplit() {
|
||||
@ -96,10 +96,10 @@ public class EigenDecompositionImplTest extends TestCase {
|
||||
{ -0.2976, 0.1152, -0.1344, 0.3872 }
|
||||
});
|
||||
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
|
||||
assertEquals(1.0, ed.getEigenvalue(0), 1.0e-15);
|
||||
assertEquals(0.4, ed.getEigenvalue(1), 1.0e-15);
|
||||
assertEquals(0.2, ed.getEigenvalue(2), 1.0e-15);
|
||||
assertEquals(0.1, ed.getEigenvalue(3), 1.0e-15);
|
||||
assertEquals(1.0, ed.getRealEigenvalue(0), 1.0e-15);
|
||||
assertEquals(0.4, ed.getRealEigenvalue(1), 1.0e-15);
|
||||
assertEquals(0.2, ed.getRealEigenvalue(2), 1.0e-15);
|
||||
assertEquals(0.1, ed.getRealEigenvalue(3), 1.0e-15);
|
||||
}
|
||||
|
||||
/** test a matrix already in tridiagonal form. */
|
||||
@ -120,7 +120,7 @@ public class EigenDecompositionImplTest extends TestCase {
|
||||
new EigenDecompositionImpl(t.getMainDiagonalRef(),
|
||||
t.getSecondaryDiagonalRef(),
|
||||
MathUtils.SAFE_MIN);
|
||||
double[] eigenValues = ed.getEigenvalues();
|
||||
double[] eigenValues = ed.getRealEigenvalues();
|
||||
assertEquals(ref.length, eigenValues.length);
|
||||
for (int i = 0; i < ref.length; ++i) {
|
||||
assertEquals(ref[ref.length - i - 1], eigenValues[i], 2.0e-14);
|
||||
@ -143,7 +143,7 @@ public class EigenDecompositionImplTest extends TestCase {
|
||||
/** test eigenvalues */
|
||||
public void testEigenvalues() {
|
||||
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
|
||||
double[] eigenValues = ed.getEigenvalues();
|
||||
double[] eigenValues = ed.getRealEigenvalues();
|
||||
assertEquals(refValues.length, eigenValues.length);
|
||||
for (int i = 0; i < refValues.length; ++i) {
|
||||
assertEquals(refValues[i], eigenValues[i], 3.0e-15);
|
||||
@ -160,7 +160,7 @@ public class EigenDecompositionImplTest extends TestCase {
|
||||
Arrays.sort(bigValues);
|
||||
EigenDecomposition ed =
|
||||
new EigenDecompositionImpl(createTestMatrix(r, bigValues), MathUtils.SAFE_MIN);
|
||||
double[] eigenValues = ed.getEigenvalues();
|
||||
double[] eigenValues = ed.getRealEigenvalues();
|
||||
assertEquals(bigValues.length, eigenValues.length);
|
||||
for (int i = 0; i < bigValues.length; ++i) {
|
||||
assertEquals(bigValues[bigValues.length - i - 1], eigenValues[i], 2.0e-14);
|
||||
@ -171,7 +171,7 @@ public class EigenDecompositionImplTest extends TestCase {
|
||||
public void testEigenvectors() {
|
||||
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
|
||||
for (int i = 0; i < matrix.getRowDimension(); ++i) {
|
||||
double lambda = ed.getEigenvalue(i);
|
||||
double lambda = ed.getRealEigenvalue(i);
|
||||
RealVector v = ed.getEigenvector(i);
|
||||
RealVector mV = matrix.operate(v);
|
||||
assertEquals(0, mV.subtract(v.mapMultiplyToSelf(lambda)).getNorm(), 1.0e-13);
|
||||
@ -201,10 +201,10 @@ public class EigenDecompositionImplTest extends TestCase {
|
||||
double[] diagonal = new double[] { -3.0, -2.0, 2.0, 5.0 };
|
||||
RealMatrix m = createDiagonalMatrix(diagonal, diagonal.length, diagonal.length);
|
||||
EigenDecomposition ed = new EigenDecompositionImpl(m, MathUtils.SAFE_MIN);
|
||||
assertEquals(diagonal[0], ed.getEigenvalue(3), 2.0e-15);
|
||||
assertEquals(diagonal[1], ed.getEigenvalue(2), 2.0e-15);
|
||||
assertEquals(diagonal[2], ed.getEigenvalue(1), 2.0e-15);
|
||||
assertEquals(diagonal[3], ed.getEigenvalue(0), 2.0e-15);
|
||||
assertEquals(diagonal[0], ed.getRealEigenvalue(3), 2.0e-15);
|
||||
assertEquals(diagonal[1], ed.getRealEigenvalue(2), 2.0e-15);
|
||||
assertEquals(diagonal[2], ed.getRealEigenvalue(1), 2.0e-15);
|
||||
assertEquals(diagonal[3], ed.getRealEigenvalue(0), 2.0e-15);
|
||||
}
|
||||
|
||||
/**
|
||||
@ -244,7 +244,7 @@ public class EigenDecompositionImplTest extends TestCase {
|
||||
*/
|
||||
protected void checkEigenValues(double[] targetValues,
|
||||
EigenDecomposition ed, double tolerance) {
|
||||
double[] observed = ed.getEigenvalues();
|
||||
double[] observed = ed.getRealEigenvalues();
|
||||
for (int i = 0; i < observed.length; i++) {
|
||||
assertTrue(isIncludedValue(observed[i], targetValues, tolerance));
|
||||
assertTrue(isIncludedValue(targetValues[i], observed, tolerance));
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user