Java tutorial
/* * Licensed to the Apache Software Foundation (ASF) under one or more * contributor license agreements. See the NOTICE file distributed with * this work for additional information regarding copyright ownership. * The ASF licenses this file to You under the Apache License, Version 2.0 * (the "License"); you may not use this file except in compliance with * the License. You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. */ package org.apache.commons.math.linear; import org.apache.commons.math.MathRuntimeException; import org.apache.commons.math.MaxIterationsExceededException; import org.apache.commons.math.exception.util.LocalizedFormats; import org.apache.commons.math.util.MathUtils; import org.apache.commons.math.util.FastMath; /** * Calculates the eigen decomposition of a real <strong>symmetric</strong> * matrix. * <p> * The eigen decomposition of matrix A is a set of two matrices: V and D such * that A = V D V<sup>T</sup>. A, V and D are all m × m matrices. * </p> * <p> * As of 2.0, this class supports only <strong>symmetric</strong> matrices, and * hence computes only real realEigenvalues. This implies the D matrix returned * by {@link #getD()} is always diagonal and the imaginary values returned * {@link #getImagEigenvalue(int)} and {@link #getImagEigenvalues()} are always * null. * </p> * <p> * When called with a {@link RealMatrix} argument, this implementation only uses * the upper part of the matrix, the part below the diagonal is not accessed at * all. * </p> * <p> * This implementation is based on the paper by A. Drubrulle, R.S. Martin and * J.H. Wilkinson 'The Implicit QL Algorithm' in Wilksinson and Reinsch (1971) * Handbook for automatic computation, vol. 2, Linear algebra, Springer-Verlag, * New-York * </p> * @version $Revision: 1002040 $ $Date: 2010-09-28 09:18:31 +0200 (mar. 28 sept. 2010) $ * @since 2.0 */ public class EigenDecompositionImpl implements EigenDecomposition { /** Maximum number of iterations accepted in the implicit QL transformation */ private byte maxIter = 30; /** Main diagonal of the tridiagonal matrix. */ private double[] main; /** Secondary diagonal of the tridiagonal matrix. */ private double[] secondary; /** * Transformer to tridiagonal (may be null if matrix is already * tridiagonal). */ private TriDiagonalTransformer transformer; /** Real part of the realEigenvalues. */ private double[] realEigenvalues; /** Imaginary part of the realEigenvalues. */ private double[] imagEigenvalues; /** Eigenvectors. */ private ArrayRealVector[] eigenvectors; /** Cached value of V. */ private RealMatrix cachedV; /** Cached value of D. */ private RealMatrix cachedD; /** Cached value of Vt. */ private RealMatrix cachedVt; /** * Calculates the eigen decomposition of the given symmetric matrix. * @param matrix The <strong>symmetric</strong> matrix to decompose. * @param splitTolerance dummy parameter, present for backward compatibility only. * @exception InvalidMatrixException (wrapping a * {@link org.apache.commons.math.ConvergenceException} if algorithm * fails to converge */ public EigenDecompositionImpl(final RealMatrix matrix, final double splitTolerance) throws InvalidMatrixException { if (isSymmetric(matrix)) { transformToTridiagonal(matrix); findEigenVectors(transformer.getQ().getData()); } 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(LocalizedFormats.ASSYMETRIC_EIGEN_NOT_SUPPORTED); } } /** * Calculates the eigen decomposition of the symmetric tridiagonal * matrix. The Householder matrix is assumed to be the identity matrix. * @param main Main diagonal of the symmetric triadiagonal form * @param secondary Secondary of the tridiagonal form * @param splitTolerance dummy parameter, present for backward compatibility only. * @exception InvalidMatrixException (wrapping a * {@link org.apache.commons.math.ConvergenceException} if algorithm * fails to converge */ public EigenDecompositionImpl(final double[] main, final double[] secondary, final double splitTolerance) throws InvalidMatrixException { this.main = main.clone(); this.secondary = secondary.clone(); transformer = null; final int size = main.length; double[][] z = new double[size][size]; for (int i = 0; i < size; i++) { z[i][i] = 1.0; } findEigenVectors(z); } /** * Check if a matrix is symmetric. * @param matrix * matrix to check * @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 (FastMath.abs(mij - mji) > (FastMath.max(FastMath.abs(mij), FastMath.abs(mji)) * eps)) { return false; } } } return true; } /** {@inheritDoc} */ public RealMatrix getV() throws InvalidMatrixException { if (cachedV == null) { final int m = eigenvectors.length; cachedV = MatrixUtils.createRealMatrix(m, m); for (int k = 0; k < m; ++k) { cachedV.setColumnVector(k, eigenvectors[k]); } } // return the cached matrix return cachedV; } /** {@inheritDoc} */ public RealMatrix getD() throws InvalidMatrixException { if (cachedD == null) { // cache the matrix for subsequent calls cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues); } return cachedD; } /** {@inheritDoc} */ public RealMatrix getVT() throws InvalidMatrixException { if (cachedVt == null) { final int m = eigenvectors.length; cachedVt = MatrixUtils.createRealMatrix(m, m); for (int k = 0; k < m; ++k) { cachedVt.setRowVector(k, eigenvectors[k]); } } // return the cached matrix return cachedVt; } /** {@inheritDoc} */ public double[] getRealEigenvalues() throws InvalidMatrixException { return realEigenvalues.clone(); } /** {@inheritDoc} */ public double getRealEigenvalue(final int i) throws InvalidMatrixException, ArrayIndexOutOfBoundsException { 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} */ public RealVector getEigenvector(final int i) throws InvalidMatrixException, ArrayIndexOutOfBoundsException { return eigenvectors[i].copy(); } /** * Return the determinant of the matrix * @return determinant of the matrix */ public double getDeterminant() { double determinant = 1; for (double lambda : realEigenvalues) { determinant *= lambda; } return determinant; } /** {@inheritDoc} */ public DecompositionSolver getSolver() { return new Solver(realEigenvalues, imagEigenvalues, eigenvectors); } /** Specialized solver. */ private static class Solver implements DecompositionSolver { /** Real part of the realEigenvalues. */ private double[] realEigenvalues; /** Imaginary part of the realEigenvalues. */ private double[] imagEigenvalues; /** Eigenvectors. */ private final ArrayRealVector[] eigenvectors; /** * Build a solver from decomposed matrix. * @param realEigenvalues * real parts of the eigenvalues * @param imagEigenvalues * imaginary parts of the eigenvalues * @param eigenvectors * eigenvectors */ private Solver(final double[] realEigenvalues, final double[] imagEigenvalues, final ArrayRealVector[] eigenvectors) { this.realEigenvalues = realEigenvalues; this.imagEigenvalues = imagEigenvalues; this.eigenvectors = eigenvectors; } /** * Solve the linear equation A × X = B for symmetric matrices A. * <p> * This method only find exact linear solutions, i.e. solutions for * which ||A × X - B|| is exactly 0. * </p> * @param b * right-hand side of the equation A × X = B * @return a vector X that minimizes the two norm of A × X - B * @exception IllegalArgumentException * if matrices dimensions don't match * @exception InvalidMatrixException * if decomposed matrix is singular */ public double[] solve(final double[] b) throws IllegalArgumentException, InvalidMatrixException { if (!isNonSingular()) { throw new SingularMatrixException(); } final int m = realEigenvalues.length; if (b.length != m) { throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.VECTOR_LENGTH_MISMATCH, b.length, m); } final double[] bp = new double[m]; for (int i = 0; i < m; ++i) { final ArrayRealVector v = eigenvectors[i]; final double[] vData = v.getDataRef(); final double s = v.dotProduct(b) / realEigenvalues[i]; for (int j = 0; j < m; ++j) { bp[j] += s * vData[j]; } } return bp; } /** * Solve the linear equation A × X = B for symmetric matrices A. * <p> * This method only find exact linear solutions, i.e. solutions for * which ||A × X - B|| is exactly 0. * </p> * @param b * right-hand side of the equation A × X = B * @return a vector X that minimizes the two norm of A × X - B * @exception IllegalArgumentException * if matrices dimensions don't match * @exception InvalidMatrixException * if decomposed matrix is singular */ public RealVector solve(final RealVector b) throws IllegalArgumentException, InvalidMatrixException { if (!isNonSingular()) { throw new SingularMatrixException(); } final int m = realEigenvalues.length; if (b.getDimension() != m) { throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.VECTOR_LENGTH_MISMATCH, b.getDimension(), m); } final double[] bp = new double[m]; for (int i = 0; i < m; ++i) { final ArrayRealVector v = eigenvectors[i]; final double[] vData = v.getDataRef(); final double s = v.dotProduct(b) / realEigenvalues[i]; for (int j = 0; j < m; ++j) { bp[j] += s * vData[j]; } } return new ArrayRealVector(bp, false); } /** * Solve the linear equation A × X = B for symmetric matrices A. * <p> * This method only find exact linear solutions, i.e. solutions for * which ||A × X - B|| is exactly 0. * </p> * @param b * right-hand side of the equation A × X = B * @return a matrix X that minimizes the two norm of A × X - B * @exception IllegalArgumentException * if matrices dimensions don't match * @exception InvalidMatrixException * if decomposed matrix is singular */ public RealMatrix solve(final RealMatrix b) throws IllegalArgumentException, InvalidMatrixException { if (!isNonSingular()) { throw new SingularMatrixException(); } final int m = realEigenvalues.length; if (b.getRowDimension() != m) { throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.DIMENSIONS_MISMATCH_2x2, b.getRowDimension(), b.getColumnDimension(), m, "n"); } 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 ArrayRealVector v = eigenvectors[i]; final double[] vData = v.getDataRef(); double s = 0; for (int j = 0; j < m; ++j) { s += v.getEntry(j) * b.getEntry(j, k); } s /= realEigenvalues[i]; for (int j = 0; j < m; ++j) { bp[j][k] += s * vData[j]; } } } return MatrixUtils.createRealMatrix(bp); } /** * Check if the decomposed matrix is non-singular. * @return true if the decomposed matrix is non-singular */ public boolean isNonSingular() { for (int i = 0; i < realEigenvalues.length; ++i) { if ((realEigenvalues[i] == 0) && (imagEigenvalues[i] == 0)) { return false; } } return true; } /** * Get the inverse of the decomposed matrix. * @return inverse matrix * @throws InvalidMatrixException * if decomposed matrix is singular */ public RealMatrix getInverse() throws InvalidMatrixException { if (!isNonSingular()) { throw new SingularMatrixException(); } final int m = realEigenvalues.length; final double[][] invData = new double[m][m]; for (int i = 0; i < m; ++i) { final double[] invI = invData[i]; for (int j = 0; j < m; ++j) { double invIJ = 0; for (int k = 0; k < m; ++k) { final double[] vK = eigenvectors[k].getDataRef(); invIJ += vK[i] * vK[j] / realEigenvalues[k]; } invI[j] = invIJ; } } return MatrixUtils.createRealMatrix(invData); } } /** * Transform matrix to tridiagonal. * @param matrix * matrix to transform */ private void transformToTridiagonal(final RealMatrix matrix) { // transform the matrix to tridiagonal transformer = new TriDiagonalTransformer(matrix); main = transformer.getMainDiagonalRef(); secondary = transformer.getSecondaryDiagonalRef(); } /** * Find eigenvalues and eigenvectors (Dubrulle et al., 1971) * @param householderMatrix Householder matrix of the transformation * to tri-diagonal form. */ private void findEigenVectors(double[][] householderMatrix) { double[][] z = householderMatrix.clone(); final int n = main.length; realEigenvalues = new double[n]; imagEigenvalues = new double[n]; double[] e = new double[n]; for (int i = 0; i < n - 1; i++) { realEigenvalues[i] = main[i]; e[i] = secondary[i]; } realEigenvalues[n - 1] = main[n - 1]; e[n - 1] = 0.0; // Determine the largest main and secondary value in absolute term. double maxAbsoluteValue = 0.0; for (int i = 0; i < n; i++) { if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) { maxAbsoluteValue = FastMath.abs(realEigenvalues[i]); } if (FastMath.abs(e[i]) > maxAbsoluteValue) { maxAbsoluteValue = FastMath.abs(e[i]); } } // Make null any main and secondary value too small to be significant if (maxAbsoluteValue != 0.0) { for (int i = 0; i < n; i++) { if (FastMath.abs(realEigenvalues[i]) <= MathUtils.EPSILON * maxAbsoluteValue) { realEigenvalues[i] = 0.0; } if (FastMath.abs(e[i]) <= MathUtils.EPSILON * maxAbsoluteValue) { e[i] = 0.0; } } } for (int j = 0; j < n; j++) { int its = 0; int m; do { for (m = j; m < n - 1; m++) { double delta = FastMath.abs(realEigenvalues[m]) + FastMath.abs(realEigenvalues[m + 1]); if (FastMath.abs(e[m]) + delta == delta) { break; } } if (m != j) { if (its == maxIter) throw new InvalidMatrixException(new MaxIterationsExceededException(maxIter)); its++; double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]); double t = FastMath.sqrt(1 + q * q); if (q < 0.0) { q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q - t); } else { q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q + t); } double u = 0.0; double s = 1.0; double c = 1.0; int i; for (i = m - 1; i >= j; i--) { double p = s * e[i]; double h = c * e[i]; if (FastMath.abs(p) >= FastMath.abs(q)) { c = q / p; t = FastMath.sqrt(c * c + 1.0); e[i + 1] = p * t; s = 1.0 / t; c = c * s; } else { s = p / q; t = FastMath.sqrt(s * s + 1.0); e[i + 1] = q * t; c = 1.0 / t; s = s * c; } if (e[i + 1] == 0.0) { realEigenvalues[i + 1] -= u; e[m] = 0.0; break; } q = realEigenvalues[i + 1] - u; t = (realEigenvalues[i] - q) * s + 2.0 * c * h; u = s * t; realEigenvalues[i + 1] = q + u; q = c * t - h; for (int ia = 0; ia < n; ia++) { p = z[ia][i + 1]; z[ia][i + 1] = s * z[ia][i] + c * p; z[ia][i] = c * z[ia][i] - s * p; } } if (t == 0.0 && i >= j) continue; realEigenvalues[j] -= u; e[j] = q; e[m] = 0.0; } } while (m != j); } //Sort the eigen values (and vectors) in increase order for (int i = 0; i < n; i++) { int k = i; double p = realEigenvalues[i]; for (int j = i + 1; j < n; j++) { if (realEigenvalues[j] > p) { k = j; p = realEigenvalues[j]; } } if (k != i) { realEigenvalues[k] = realEigenvalues[i]; realEigenvalues[i] = p; for (int j = 0; j < n; j++) { p = z[j][i]; z[j][i] = z[j][k]; z[j][k] = p; } } } // Determine the largest eigen value in absolute term. maxAbsoluteValue = 0.0; for (int i = 0; i < n; i++) { if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) { maxAbsoluteValue = FastMath.abs(realEigenvalues[i]); } } // Make null any eigen value too small to be significant if (maxAbsoluteValue != 0.0) { for (int i = 0; i < n; i++) { if (FastMath.abs(realEigenvalues[i]) < MathUtils.EPSILON * maxAbsoluteValue) { realEigenvalues[i] = 0.0; } } } eigenvectors = new ArrayRealVector[n]; double[] tmp = new double[n]; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { tmp[j] = z[j][i]; } eigenvectors[i] = new ArrayRealVector(tmp); } } }