Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

added a first version of eigen decomposition implementation 

this version is not finished yet, but it does work when
eigenvalues are well separated and is faster than JAMA for
dimensions above 100.
It still needs work as the MRRR algorithm is not implemented
yet (only the basic parts with twisted factorization is there).
I continue working on this, but wanted to have a first version
committed to let people play with it and as a basis for comparison.

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/branches/MATH_2_0@701897 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2008-10-05 22:38:13 +00:00
parent dc4f70a777
commit 924041770d
4 changed files with 1002 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
src/java/org/apache/commons/math/linear/EigenDecomposition.java
View File

@ -35,7 +35,11 @@ package org.apache.commons.math.linear;
* been added (in the superinterface),</li> * been added (in the superinterface),</li>
* <li>a {@link DecompositionSolver#getInverse() getInverse} method has been * <li>a {@link DecompositionSolver#getInverse() getInverse} method has been
* added (in the superinterface),</li> * added (in the superinterface),</li>
* <li>a {@link #getVt() getVt} method has been added,</li> * <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>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 * <li>the <code>getRealEigenvalues</code> method has been renamed as {@link
* #getEigenValues() getEigenValues},</li> * #getEigenValues() getEigenValues},</li>
* <li>the <code>getImagEigenvalues</code> method has been removed</li> * <li>the <code>getImagEigenvalues</code> method has been removed</li>
@ -76,15 +80,43 @@ public interface EigenDecomposition extends DecompositionSolver {
* @exception IllegalStateException if {@link * @exception IllegalStateException if {@link
* DecompositionSolver#decompose(RealMatrix) decompose} has not been called * DecompositionSolver#decompose(RealMatrix) decompose} has not been called
*/ */
RealMatrix getVt() throws IllegalStateException; RealMatrix getVT() throws IllegalStateException;
/** /**
* Returns the eigenvalues of the original matrix. * Returns a copy of the eigenvalues of the original matrix.
* @return eigenvalues of the original matrix * @return a copy of the eigenvalues of the original matrix
* @exception IllegalStateException if {@link * @exception IllegalStateException if {@link
* DecompositionSolver#decompose(RealMatrix) decompose} has not been called * DecompositionSolver#decompose(RealMatrix) decompose} has not been called
* @see #getD() * @see #getD()
*/ */
double[] getEigenValues() throws IllegalStateException; double[] getEigenvalues() throws IllegalStateException;
/**
* Returns the i<sup>th</sup> eigenvalue of the original matrix.
* @return i<sup>th</sup> eigenvalue of the original matrix
* @exception IllegalStateException if {@link
* DecompositionSolver#decompose(RealMatrix) decompose} has not been called
* @exception ArrayIndexOutOfBoundsException if i is not
* @see #getD()
*/
double getEigenvalue(int i) throws IllegalStateException;
/**
* Returns a copy of the i<sup>th</sup> eigenvector of the original matrix.
* @return copy of the i<sup>th</sup> eigenvector of the original matrix
* @exception IllegalStateException if {@link
* DecompositionSolver#decompose(RealMatrix) decompose} has not been called
* @see #getD()
*/
RealVector getEigenvector(int i) throws IllegalStateException;
/**
* Return the determinant of the matrix
* @return determinant of the matrix
* @exception IllegalStateException if {@link
* DecompositionSolver#decompose(RealMatrix) decompose} has not been called
* @see #isNonSingular()
*/
double getDeterminant() throws IllegalStateException;
} }

658
src/java/org/apache/commons/math/linear/EigenDecompositionImpl.java Normal file
View File

@ -0,0 +1,658 @@
/*
* 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 java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.util.SortedSet;
import java.util.TreeSet;
import org.apache.commons.math.ConvergenceException;
/**
* Calculates the eigen decomposition of a matrix.
* <p>The eigen decomposition of matrix A is a set of two matrices:
* V and D such that A = V &times; D &times; V<sup>T</sup>.
* Let A be an m &times; n matrix, then V is an m &times; m orthogonal matrix
* and D is a m &times; n diagonal matrix.</p>
* <p>This implementation is based on Inderjit Singh Dhillon thesis
* <a href="http://www.cs.utexas.edu/users/inderjit/public_papers/thesis.pdf">A
* New O(n<sup>2</sup>) Algorithm for the Symmetric Tridiagonal Eigenvalue/Eigenvector
* Problem</a>.</p>
* @version $Revision$ $Date$
* @since 2.0
*/
public class EigenDecompositionImpl implements EigenDecomposition {
/** Serializable version identifier. */
private static final long serialVersionUID = -8550254713195393577L;
/** Eigenvalues. */
private double[] eigenvalues;
/** Eigenvectors. */
private RealVectorImpl[] eigenvectors;
/** Cached value of V. */
private RealMatrix cachedV;
/** Cached value of D. */
private RealMatrix cachedD;
/** Cached value of Vt. */
private RealMatrix cachedVt;
/**
* Build a new instance.
* <p>Note that {@link #decompose(RealMatrix)} <strong>must</strong> be called
* before any of the {@link #getV()}, {@link #getD()}, {@link #getVT()},
* {@link #getEignevalues()}, {@link #solve(double[])}, {@link #solve(RealMatrix)},
* {@link #solve(RealVector)} or {@link #solve(RealVectorImpl)} methods can be
* called.</p>
* @see #decompose(RealMatrix)
*/
public EigenDecompositionImpl() {
}
/**
* Calculates the eigen decomposition of the given matrix.
* <p>Calling this constructor is equivalent to first call the no-arguments
* constructor and then call {@link #decompose(RealMatrix)}.</p>
* @param matrix The matrix to decompose.
* @exception InvalidMatrixException (wrapping a {@link ConvergenceException}
* if algorithm fails to converge
*/
public EigenDecompositionImpl(RealMatrix matrix)
throws InvalidMatrixException {
decompose(matrix);
}
/** {@inheritDoc} */
public void decompose(RealMatrix matrix)
throws InvalidMatrixException {
cachedV = null;
cachedD = null;
cachedVt = null;
// transform the matrix to tridiagonal
TriDiagonalTransformer transformer = new TriDiagonalTransformer(matrix);
final double[] main = transformer.getMainDiagonalRef();
final double[] secondary = transformer.getSecondaryDiagonalRef();
final int m = main.length;
// pre-compute the square of the secondary diagonal
double[] squaredSecondary = new double[secondary.length];
for (int i = 0; i < squaredSecondary.length; ++i) {
final double s = secondary[i];
squaredSecondary[i] = s * s;
}
// compute the eigenvalues bounds
List<GershgorinCirclesUnion> bounds =
getEigenvaluesBounds(main, secondary);
// TODO this implementation is not finished yet
// the MRRR algorithm is NOT implemented, Gershgorin circles are
// merged together when they could be separated, we only perform blindly
// the basic steps, we search all eigenvalues with an arbitrary
// threshold, we use twisted factorization afterwards with no
// heuristic to speed up the selection of the twist index ...
// The decomposition does work in its current state and seems reasonably
// efficient when eigenvalues are separated. However, it is expected to
// fail in difficult cases and its performances can obviously be improved
// for now, it is slower than JAMA for dimensions below 100 and faster
// for dimensions above 100. The speed gain with respect to JAMA increase
// regularly with dimension
// find eigenvalues using bisection
eigenvalues = new double[m];
final double low = bounds.get(0).getLow();
final double high = bounds.get(bounds.size() - 1).getHigh();
final double threshold =
1.0e-15 * Math.max(Math.abs(low), Math.abs(high));
findEigenvalues(main, squaredSecondary, low, high, threshold, 0, m);
// find eigenvectors
eigenvectors = new RealVectorImpl[m];
final double[] eigenvector = new double[m];
final double[] lp = new double[m - 1];
final double[] dp = new double[m];
final double[] um = new double[m - 1];
final double[] dm = new double[m];
final double[] gamma = new double[m];
for (int i = 0; i < m; ++i) {
// find the eigenvector of the tridiagonal matrix
findEigenvector(eigenvalues[i], eigenvector,
main, secondary, lp, dp, um, dm, gamma);
// find the eigenvector of the original matrix
eigenvectors[i] =
new RealVectorImpl(transformer.getQ().operate(eigenvector), true);
}
}
/** {@inheritDoc} */
public RealMatrix getV()
throws InvalidMatrixException {
if (cachedV == null) {
cachedV = getVT().transpose();
}
// return the cached matrix
return cachedV;
}
/** {@inheritDoc} */
public RealMatrix getD()
throws InvalidMatrixException {
if (cachedD == null) {
checkDecomposed();
final int m = eigenvalues.length;
final double[][] sData = new double[m][m];
for (int i = 0; i < m; ++i) {
sData[i][i] = eigenvalues[i];
}
// cache the matrix for subsequent calls
cachedD = new RealMatrixImpl(sData, false);
}
return cachedD;
}
/** {@inheritDoc} */
public RealMatrix getVT()
throws InvalidMatrixException {
if (cachedVt == null) {
checkDecomposed();
final double[][] vtData = new double[eigenvectors.length][];
for (int k = 0; k < eigenvectors.length; ++k) {
vtData[k] = eigenvectors[k].getData();
}
// cache the matrix for subsequent calls
cachedVt = new RealMatrixImpl(vtData, false);
}
// return the cached matrix
return cachedVt;
}
/** {@inheritDoc} */
public double[] getEigenvalues()
throws InvalidMatrixException {
checkDecomposed();
return eigenvalues.clone();
}
/** {@inheritDoc} */
public double getEigenvalue(final int i)
throws InvalidMatrixException, ArrayIndexOutOfBoundsException {
checkDecomposed();
return eigenvalues[i];
}
/** {@inheritDoc} */
public RealVector getEigenvector(final int i)
throws InvalidMatrixException, ArrayIndexOutOfBoundsException {
checkDecomposed();
return eigenvectors[i].copy();
}
/** {@inheritDoc} */
public boolean isNonSingular()
throws IllegalStateException {
for (double lambda : eigenvalues) {
if (lambda == 0) {
return false;
}
}
return true;
}
/** {@inheritDoc} */
public double[] solve(final double[] b)
throws IllegalArgumentException, InvalidMatrixException {
checkNonSingular();
final int m = eigenvalues.length;
if (b.length != m) {
throw new IllegalArgumentException("constant vector has wrong length");
}
final double[] bp = new double[m];
for (int i = 0; i < m; ++i) {
final RealVectorImpl v = eigenvectors[i];
final double s = v.dotProduct(b) / eigenvalues[i];
final double[] vData = v.getDataRef();
for (int j = 0; j < m; ++j) {
bp[j] += s * vData[j];
}
}
return bp;
}
/** {@inheritDoc} */
public RealVector solve(final RealVector b)
throws IllegalArgumentException, InvalidMatrixException {
try {
return solve((RealVectorImpl) b);
} catch (ClassCastException cce) {
checkNonSingular();
final int m = eigenvalues.length;
if (b.getDimension() != m) {
throw new IllegalArgumentException("constant vector has wrong length");
}
final double[] bp = new double[m];
for (int i = 0; i < m; ++i) {
final RealVectorImpl v = eigenvectors[i];
final double s = v.dotProduct(b) / eigenvalues[i];
final double[] vData = v.getDataRef();
for (int j = 0; j < m; ++j) {
bp[j] += s * vData[j];
}
}
return new RealVectorImpl(bp, false);
}
}
/** Solve the linear equation A &times; X = B.
* <p>The A matrix is implicit here. It is </p>
* @param b right-hand side of the equation A &times; X = B
* @return a vector X such that A &times; X = B
* @throws IllegalArgumentException if matrices dimensions don't match
* @throws InvalidMatrixException if decomposed matrix is singular
*/
public RealVectorImpl solve(final RealVectorImpl b)
throws IllegalArgumentException, InvalidMatrixException {
return new RealVectorImpl(solve(b.getDataRef()), false);
}
/** {@inheritDoc} */
public RealMatrix solve(final RealMatrix b)
throws IllegalArgumentException, InvalidMatrixException {
checkNonSingular();
final int m = eigenvalues.length;
if (b.getRowDimension() != m) {
throw new IllegalArgumentException("Incorrect row dimension");
}
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 double[] vData = eigenvectors[i].getDataRef();
double s = 0;
for (int j = 0; j < m; ++j) {
s += vData[j] * b.getEntry(j, k);
}
s /= eigenvalues[i];
for (int j = 0; j < m; ++j) {
bp[j][k] += s * vData[j];
}
}
}
return new RealMatrixImpl(bp, false);
}
/** {@inheritDoc} */
public RealMatrix getInverse()
throws IllegalStateException, InvalidMatrixException {
checkNonSingular();
final int m = eigenvalues.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] / eigenvalues[k];
}
invI[j] = invIJ;
}
}
return new RealMatrixImpl(invData, false);
}
/** {@inheritDoc} */
public double getDeterminant()
throws IllegalStateException {
double determinant = 1;
for (double lambda : eigenvalues) {
determinant *= lambda;
}
return determinant;
}
/**
* Compute a set of possible bounding intervals for eigenvalues
* of a symmetric tridiagonal matrix.
* <p>The intervals are computed by applying the Gershgorin circle theorem.</p>
* @param main main diagonal
* @param secondary secondary diagonal of the tridiagonal matrix
* @return a collection of disjoint intervals where eigenvalues must lie,
* sorted in increasing order
*/
private List<GershgorinCirclesUnion> getEigenvaluesBounds(final double[] main,
final double[] secondary) {
final SortedSet<GershgorinCirclesUnion> rawCircles =
new TreeSet<GershgorinCirclesUnion>();
final int m = main.length;
// compute all the Gershgorin circles independently
rawCircles.add(new GershgorinCirclesUnion(main[0],
Math.abs(secondary[0])));
for (int i = 1; i < m - 1; ++i) {
rawCircles.add(new GershgorinCirclesUnion(main[i],
Math.abs(secondary[i - 1]) +
Math.abs(secondary[i])));
}
rawCircles.add(new GershgorinCirclesUnion(main[m - 1],
Math.abs(secondary[m - 2])));
// combine intersecting circles
final ArrayList<GershgorinCirclesUnion> combined =
new ArrayList<GershgorinCirclesUnion>();
GershgorinCirclesUnion current = null;
for (GershgorinCirclesUnion rawCircle : rawCircles) {
if (current == null) {
current = rawCircle;
} else if (current.intersects(rawCircle)) {
current.swallow(rawCircle);
} else {
combined.add(current);
current = rawCircle;
}
}
if (current != null) {
combined.add(current);
}
return combined;
}
/** Find eigenvalues in an interval.
* @param main main diagonal of the tridiagonal matrix
* @param squaredSecondary squared secondary diagonal of the tridiagonal matrix
* @param low lower bound of the search interval
* @param high higher bound of the search interval
* @param threshold convergence threshold
* @param iStart index of the first eigenvalue to find
* @param iEnd index one unit past the last eigenvalue to find
*/
private void findEigenvalues(final double[] main, final double[] squaredSecondary,
final double low, final double high, final double threshold,
final int iStart, final int iEnd) {
// use a simple loop to handle tail-recursion cases
double currentLow = low;
double currentHigh = high;
int currentStart = iStart;
while (true) {
final double middle = 0.5 * (currentLow + currentHigh);
if (currentHigh - currentLow < threshold) {
// we have found an elementary interval containing one or more eigenvalues
Arrays.fill(eigenvalues, currentStart, iEnd, middle);
return;
}
// compute the number of eigenvalues below the middle interval point
final int iMiddle = countEigenValues(main, squaredSecondary, middle);
if (iMiddle == currentStart) {
// all eigenvalues are in the upper half of the search interval
// update the interval and iterate
currentLow = middle;
} else if (iMiddle == iEnd) {
// all eigenvalues are in the lower half of the search interval
// update the interval and iterate
currentHigh = middle;
} else {
// split the interval and search eigenvalues in both sub-intervals
findEigenvalues(main, squaredSecondary, currentLow, middle, threshold,
currentStart, iMiddle);
currentLow = middle;
currentStart = iMiddle;
}
}
}
/**
* Count the number of eigenvalues below a point.
* @param main main diagonal of the tridiagonal matrix
* @param squaredSecondary squared secondary diagonal of the tridiagonal matrix
* @param mu value below which we must count the number of eigenvalues
* @return number of eigenvalues smaller than mu
*/
private int countEigenValues(final double[] main, final double[] squaredSecondary,
final double mu) {
double ratio = main[0] - mu;
int count = (ratio > 0) ? 0 : 1;
for (int i = 1; i < main.length; ++i) {
ratio = main[i] - squaredSecondary[i - 1] / ratio - mu;
if (ratio <= 0) {
++count;
}
}
return count;
}
/**
* Decompose the shifted tridiagonal matrix A - lambda I as L<sub>+</sub>
* &times; D<sub>+</sub> &times; U<sub>+</sub>.
* <p>A shifted symmetric tridiagonal matrix can be decomposed as
* L<sub>+</sub> &times; D<sub>+</sub> &times; U<sub>+</sub> where L<sub>+</sub>
* is a lower bi-diagonal matrix with unit diagonal, D<sub>+</sub> is a diagonal
* matrix and U<sub>+</sub> is the transpose of L<sub>+</sub>. The '+' indice
* comes from Dhillon's notation since decomposition is done in
* increasing rows order).</p>
* @param main main diagonal of the tridiagonal matrix
* @param secondary secondary diagonal of the tridiagonal matrix
* @param lambda shift to apply to the matrix before decomposing it
* @param r index at which factorization should stop (if r is
* <code>main.length</code>, complete factorization is performed)
* @param lp placeholder where to put the (r-1) first off-diagonal
* elements of the L<sub>+</sub> matrix
* @param dp placeholder where to put the r first diagonal elements
* of the D<sub>+</sub> matrix
*/
private void lduDecomposition(final double[] main, final double[] secondary,
final double lambda, final int r,
final double[] lp, final double[] dp) {
double di = main[0] - lambda;
dp[0] = di;
for (int i = 1; i < r; ++i) {
final double eiM1 = secondary[i - 1];
final double ratio = eiM1 / di;
di = main[i] - lambda - eiM1 * ratio;
lp[i - 1] = ratio;
dp[i] = di;
}
}
/**
* Decompose the shifted tridiagonal matrix A - lambda I as U<sub>-</sub>
* &times; D<sub>-</sub> &times; L<sub>-</sub>.
* <p>A shifted symmetric tridiagonal matrix can be decomposed as
* U<sub>-</sub> &times; D<sub>-</sub> &times; L<sub>-</sub> where U<sub>-</sub>
* is an upper bi-diagonal matrix with unit diagonal, D<sub>-</sub> is a diagonal
* matrix and L<sub>-</sub> is the transpose of U<sub>-</sub>. The '-' indice
* comes from Dhillon's notation since decomposition is done in
* decreasing rows order).</p>
* @param main main diagonal of the tridiagonal matrix
* @param secondary secondary diagonal of the tridiagonal matrix
* @param lambda shift to apply to the matrix before decomposing it
* @param r index at which factorization should stop (if r is 0, complete
* factorization is performed)
* @param um placeholder where to put the m-(r-1) last off-diagonal elements
* of the U<sub>-</sub> matrix, where m is the size of the original matrix
* @param dm placeholder where to put the m-r last diagonal elements
* of the D<sub>-</sub> matrix, where m is the size of the original matrix
*/
private void udlDecomposition(final double[] main, final double[] secondary,
final double lambda, final int r,
final double[] um, final double[] dm) {
final int mM1 = main.length - 1;
double di = main[mM1] - lambda;
dm[mM1] = di;
for (int i = mM1 - 1; i >= r; --i) {
final double ei = secondary[i];
final double ratio = ei / di;
di = main[i] - lambda - ei * ratio;
um[i] = ratio;
dm[i] = di;
}
}
/**
* Find an eigenvector corresponding to an eigenvalue.
* @param eigenvalue eigenvalue for which eigenvector is desired
* @param eigenvector placeholder where to put the eigenvector
* @param main main diagonal of the tridiagonal matrix
* @param secondary secondary diagonal of the tridiagonal matrix
* @param lp placeholder where to put the off-diagonal elements of the
* L<sub>+</sub> matrix
* @param dp placeholder where to put the diagonal elements of the
* D<sub>+</sub> matrix
* @param um placeholder where to put the off-diagonal elements of the
* U<sub>-</sub> matrix
* @param dm placeholder where to put the diagonal elements of the
* D<sub>-</sub> matrix
* @param gamma placeholder where to put the twist elements for all
* possible twist indices
*/
private void findEigenvector(final double eigenvalue, final double[] eigenvector,
final double[] main, final double[] secondary,
final double[] lp, final double[] dp,
final double[] um, final double[] dm,
final double[] gamma) {
// compute the LDU and UDL decomposition of the
// perfectly shifted tridiagonal matrix
final int m = main.length;
lduDecomposition(main, secondary, eigenvalue, m, lp, dp);
udlDecomposition(main, secondary, eigenvalue, 0, um, dm);
// select the twist index leading to
// the least diagonal element in the twisted factorization
int r = 0;
double g = dp[0] + dm[0] + eigenvalue - main[0];
gamma[0] = g;
double minG = Math.abs(g);
for (int i = 1; i < m; ++i) {
if (i < m - 1) {
g *= dm[i + 1] / dp[i];
} else {
g = dp[m - 1] + dm[m - 1] + eigenvalue - main[m - 1];
}
gamma[i] = g;
final double absG = Math.abs(g);
if (absG < minG) {
r = i;
minG = absG;
}
}
// solve the singular system by ignoring the equation
// at twist index and propagating upwards and downwards
double n2 = 1;
eigenvector[r] = 1;
double z = 1;
for (int i = r - 1; i >= 0; --i) {
z *= -lp[i];
eigenvector[i] = z;
n2 += z * z;
}
z = 1;
for (int i = r + 1; i < m; ++i) {
z *= -um[i-1];
eigenvector[i] = z;
n2 += z * z;
}
// normalize vector
final double inv = 1.0 / Math.sqrt(n2);
for (int i = 0; i < m; ++i) {
eigenvector[i] *= inv;
}
}
/**
* Check if decomposition has been performed.
* @exception IllegalStateException if {@link #decompose(RealMatrix) decompose}
* has not been called
*/
private void checkDecomposed()
throws IllegalStateException {
if (eigenvalues == null) {
throw new IllegalStateException("no matrix have been decomposed yet");
}
}
/**
* Check if decomposed matrix is non singular.
* @exception IllegalStateException if {@link #decompose(RealMatrix) decompose}
* has not been called
* @exception InvalidMatrixException if decomposed matrix is singular
*/
private void checkNonSingular()
throws IllegalStateException, InvalidMatrixException {
checkDecomposed();
if (!isNonSingular()) {
throw new IllegalStateException("matrix is singular");
}
}
}

89
src/java/org/apache/commons/math/linear/GershgorinCirclesUnion.java Normal file
View File

@ -0,0 +1,89 @@
/*
* 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;
/** Class representing a union of Gershgorin circles.
* <p>Gershgorin circles are bounding areas where eigenvalues must lie.
* They are used as starting values for eigen decomposition algorithms.
* In the real case, Gershgorin circles are simple intervals.</p>
* @see EigenDecompositionImpl
* @version $Revision$ $Date$
* @since 2.0
*/
class GershgorinCirclesUnion implements Comparable<GershgorinCirclesUnion> {
/** Lower bound of the interval. */
private double low;
/** Higher bound of the interval. */
private double high;
/** Create a simple Gershgorin circle.
* @param d diagonal element of the current row
* @param sum sum of the absolute values of the off-diagonal elements
* of the current row
*/
public GershgorinCirclesUnion(final double d, final double sum) {
low = d - sum;
high = d + sum;
}
/**
* Get the lower bound of the interval.
* @return lower bound of the interval
*/
public double getLow() {
return low;
}
/**
* Get the higher bound of the interval.
* @return higher bound of the interval
*/
public double getHigh() {
return high;
}
/**
* Check if a Gershgorin circles union intersects instance.
* @param other Gershgorin circles union to test against instance
* @return true if the other Gershgorin circles union intersects instance
*/
public boolean intersects(final GershgorinCirclesUnion other) {
return (other.low <= this.high) && (other.high >= this.low);
}
/**
* Swallow another Gershgorin circles union.
* <p>Swallowing another Gershgorin circles union changes the
* instance such that it contains everything that was formerly in
* either circles union. It is mainly intended for circles unions
* that {@link #intersects(GershgorinCirclesUnion) intersect}
* each other beforehand.</p>
*/
public void swallow(final GershgorinCirclesUnion other) {
low = Math.min(low, other.low);
high = Math.max(high, other.high);
}
/** Comparator class for sorting intervals. */
public int compareTo(GershgorinCirclesUnion other) {
return Double.compare(low, other.low);
}
}

218
src/test/org/apache/commons/math/linear/EigenDecompositionImplTest.java Normal file
View File

@ -0,0 +1,218 @@
/*
* 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 java.util.Random;
import junit.framework.Test;
import junit.framework.TestCase;
import junit.framework.TestSuite;
public class EigenDecompositionImplTest extends TestCase {
private double[] refValues;
private RealMatrix matrix;
private static final double normTolerance = 1.e-10;
public EigenDecompositionImplTest(String name) {
super(name);
}
public static Test suite() {
TestSuite suite = new TestSuite(EigenDecompositionImplTest.class);
suite.setName("EigenDecompositionImpl Tests");
return suite;
}
/** test dimensions */
public void testDimensions() {
final int m = matrix.getRowDimension();
EigenDecomposition ed = new EigenDecompositionImpl(matrix);
assertEquals(m, ed.getV().getRowDimension());
assertEquals(m, ed.getV().getColumnDimension());
assertEquals(m, ed.getD().getColumnDimension());
assertEquals(m, ed.getD().getColumnDimension());
assertEquals(m, ed.getVT().getRowDimension());
assertEquals(m, ed.getVT().getColumnDimension());
}
/** test eigenvalues */
public void testEigenvalues() {
EigenDecomposition ed = new EigenDecompositionImpl(matrix);
double[] eigenValues = ed.getEigenvalues();
assertEquals(refValues.length, eigenValues.length);
for (int i = 0; i < refValues.length; ++i) {
assertEquals(refValues[i], eigenValues[eigenValues.length - 1 - i], 3.0e-15);
}
}
/** test eigenvectors */
public void testEigenvectors() {
EigenDecomposition ed = new EigenDecompositionImpl(matrix);
for (int i = 0; i < matrix.getRowDimension(); ++i) {
double lambda = ed.getEigenvalue(i);
RealVector v = ed.getEigenvector(i);
RealVector mV = matrix.operate(v);
assertEquals(0, mV.subtract(v.mapMultiplyToSelf(lambda)).getNorm(), 1.0e-13);
}
}
/** test A = VDVt */
public void testAEqualVDVt() {
EigenDecomposition ed = new EigenDecompositionImpl(matrix);
RealMatrix v = ed.getV();
RealMatrix d = ed.getD();
RealMatrix vT = ed.getVT();
double norm = v.multiply(d).multiply(vT).subtract(matrix).getNorm();
assertEquals(0, norm, normTolerance);
}
/** test that V is orthogonal */
public void testVOrthogonal() {
RealMatrix v = new EigenDecompositionImpl(matrix).getV();
RealMatrix vTv = v.transpose().multiply(v);
RealMatrix id = MatrixUtils.createRealIdentityMatrix(vTv.getRowDimension());
assertEquals(0, vTv.subtract(id).getNorm(), normTolerance);
}
/** test solve dimension errors */
public void testSolveDimensionErrors() {
EigenDecomposition ed = new EigenDecompositionImpl(matrix);
RealMatrix b = new RealMatrixImpl(new double[2][2]);
try {
ed.solve(b);
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
ed.solve(b.getColumn(0));
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
ed.solve(new RealVectorImplTest.RealVectorTestImpl(b.getColumn(0)));
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test solve */
public void testSolve() {
RealMatrix m = new RealMatrixImpl(new double[][] {
{ 91, 5, 29, 32, 40, 14 },
{ 5, 34, -1, 0, 2, -1 },
{ 29, -1, 12, 9, 21, 8 },
{ 32, 0, 9, 14, 9, 0 },
{ 40, 2, 21, 9, 51, 19 },
{ 14, -1, 8, 0, 19, 14 }
});
EigenDecomposition ed = new EigenDecompositionImpl(m);
assertEquals(184041, ed.getDeterminant(), 2.0e-8);
RealMatrix b = new RealMatrixImpl(new double[][] {
{ 1561, 269, 188 },
{ 69, -21, 70 },
{ 739, 108, 63 },
{ 324, 86, 59 },
{ 1624, 194, 107 },
{ 796, 69, 36 }
});
RealMatrix xRef = new RealMatrixImpl(new double[][] {
{ 1, 2, 1 },
{ 2, -1, 2 },
{ 4, 2, 3 },
{ 8, -1, 0 },
{ 16, 2, 0 },
{ 32, -1, 0 }
});
// using RealMatrix
assertEquals(0, ed.solve(b).subtract(xRef).getNorm(), normTolerance);
// using double[]
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
new RealVectorImpl(ed.solve(b.getColumn(i))).subtract(xRef.getColumnVector(i)).getNorm(),
2.0e-11);
}
// using RealMatrixImpl
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
ed.solve(b.getColumnVector(i)).subtract(xRef.getColumnVector(i)).getNorm(),
2.0e-11);
}
// using RealMatrix with an alternate implementation
for (int i = 0; i < b.getColumnDimension(); ++i) {
RealVectorImplTest.RealVectorTestImpl v =
new RealVectorImplTest.RealVectorTestImpl(b.getColumn(i));
assertEquals(0,
ed.solve(v).subtract(xRef.getColumnVector(i)).getNorm(),
2.0e-11);
}
}
public void setUp() {
refValues = new double[] {
2.003, 2.002, 2.001, 1.001, 1.000, 0.001
};
matrix = createTestMatrix(new Random(35992629946426l), refValues);
}
public void tearDown() {
refValues = null;
matrix = null;
}
private RealMatrix createTestMatrix(final Random r, final double[] eigenValues) {
final RealMatrix v = createOrthogonalMatrix(r, eigenValues.length);
final RealMatrix d = createDiagonalMatrix(eigenValues, eigenValues.length);
return v.multiply(d).multiply(v.transpose());
}
private RealMatrix createOrthogonalMatrix(final Random r, final int size) {
final double[][] data = new double[size][size];
for (int i = 0; i < size; ++i) {
for (int j = 0; j < size; ++j) {
data[i][j] = 2 * r.nextDouble() - 1;
}
}
final RealMatrix m = new RealMatrixImpl(data, false);
return new QRDecompositionImpl(m).getQ();
}
private RealMatrix createDiagonalMatrix(final double[] data, final int rows) {
final double[][] dData = new double[rows][rows];
for (int i = 0; i < data.length; ++i) {
dData[i][i] = data[i];
}
return new RealMatrixImpl(dData, false);
}
}