getRealEigenvalues method has been renamed as {@link
* #getEigenValues() getEigenValues},getImagEigenvalues method has been removedThe eigen decomposition of matrix A is a set of two matrices: + * V and D such that A = V × D × VT. + * Let A be an m × n matrix, then V is an m × m orthogonal matrix + * and D is a m × n diagonal matrix.
+ *This implementation is based on Inderjit Singh Dhillon thesis + * A + * New O(n2) Algorithm for the Symmetric Tridiagonal Eigenvalue/Eigenvector + * Problem.
+ * @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. + *Note that {@link #decompose(RealMatrix)} must 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.
+ * @see #decompose(RealMatrix) + */ + public EigenDecompositionImpl() { + } + + /** + * Calculates the eigen decomposition of the given matrix. + *Calling this constructor is equivalent to first call the no-arguments + * constructor and then call {@link #decompose(RealMatrix)}.
+ * @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 + ListThe A matrix is implicit here. It is
+ * @param b right-hand side of the equation A × X = B + * @return a vector X such that A × 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. + *The intervals are computed by applying the Gershgorin circle theorem.
+ * @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 ListA shifted symmetric tridiagonal matrix can be decomposed as + * L+ × D+ × U+ where L+ + * is a lower bi-diagonal matrix with unit diagonal, D+ is a diagonal + * matrix and U+ is the transpose of L+. The '+' indice + * comes from Dhillon's notation since decomposition is done in + * increasing rows order).
+ * @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 + *main.length, complete factorization is performed)
+ * @param lp placeholder where to put the (r-1) first off-diagonal
+ * elements of the L+ matrix
+ * @param dp placeholder where to put the r first diagonal elements
+ * of the D+ 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-
+ * × D- × L-.
+ * A shifted symmetric tridiagonal matrix can be decomposed as + * U- × D- × L- where U- + * is an upper bi-diagonal matrix with unit diagonal, D- is a diagonal + * matrix and L- is the transpose of U-. The '-' indice + * comes from Dhillon's notation since decomposition is done in + * decreasing rows order).
+ * @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- 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- 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+ matrix + * @param dp placeholder where to put the diagonal elements of the + * D+ matrix + * @param um placeholder where to put the off-diagonal elements of the + * U- matrix + * @param dm placeholder where to put the diagonal elements of the + * D- 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"); + } + } + +} diff --git a/src/java/org/apache/commons/math/linear/GershgorinCirclesUnion.java b/src/java/org/apache/commons/math/linear/GershgorinCirclesUnion.java new file mode 100644 index 000000000..2f88aa89e --- /dev/null +++ b/src/java/org/apache/commons/math/linear/GershgorinCirclesUnion.java @@ -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. + *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.
+ * @see EigenDecompositionImpl + * @version $Revision$ $Date$ + * @since 2.0 + */ +class GershgorinCirclesUnion implements ComparableSwallowing 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.
+ */ + 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); + } + +} diff --git a/src/test/org/apache/commons/math/linear/EigenDecompositionImplTest.java b/src/test/org/apache/commons/math/linear/EigenDecompositionImplTest.java new file mode 100644 index 000000000..d13e132a8 --- /dev/null +++ b/src/test/org/apache/commons/math/linear/EigenDecompositionImplTest.java @@ -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); + } + +}