• a {@link DecompositionSolver#getInverse() getInverse} method has been * added (in the superinterface),
• - *
• + *
• + *
• a {@link #getEigenvalue(int) getEigenvalue} method to pick up a single + * eigenvalue has been added,
• + *
• a {@link #getEigenvector(int) getEigenvector} method to pick up a single + * eigenvector has been added,
• *
• the `getRealEigenvalues` method has been renamed as {@link * #getEigenValues() getEigenValues},
• *
• the `getImagEigenvalues` method has been removed

The 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. + *

+ * @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 + List 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 × X = B. + *

The 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 List getEigenvaluesBounds(final double[] main, + final double[] secondary) { + + final SortedSet rawCircles = + new TreeSet(); + 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 combined = + new ArrayList(); + 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+ + * × D+ × U+. + *

A 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).