    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017
    018    package org.apache.commons.math.linear;
    019
    020    import org.apache.commons.math.MathRuntimeException;
    021    import org.apache.commons.math.MaxIterationsExceededException;
    022    import org.apache.commons.math.exception.util.LocalizedFormats;
    023    import org.apache.commons.math.util.MathUtils;
    024    import org.apache.commons.math.util.FastMath;
    025
    026    /**
    027     * Calculates the eigen decomposition of a real <strong>symmetric</strong>
    028     * matrix.
    029     * <p>
    030     * The eigen decomposition of matrix A is a set of two matrices: V and D such
    031     * that A = V D V<sup>T</sup>. A, V and D are all m &times; m matrices.
    032     * </p>
    033     * <p>
    034     * As of 2.0, this class supports only <strong>symmetric</strong> matrices, and
    035     * hence computes only real realEigenvalues. This implies the D matrix returned
    036     * by {@link #getD()} is always diagonal and the imaginary values returned
    037     * {@link #getImagEigenvalue(int)} and {@link #getImagEigenvalues()} are always
    038     * null.
    039     * </p>
    040     * <p>
    041     * When called with a {@link RealMatrix} argument, this implementation only uses
    042     * the upper part of the matrix, the part below the diagonal is not accessed at
    043     * all.
    044     * </p>
    045     * <p>
    046     * This implementation is based on the paper by A. Drubrulle, R.S. Martin and
    047     * J.H. Wilkinson 'The Implicit QL Algorithm' in Wilksinson and Reinsch (1971)
    048     * Handbook for automatic computation, vol. 2, Linear algebra, Springer-Verlag,
    049     * New-York
    050     * </p>
    051     * @version $Revision: 1002040 $ $Date: 2010-09-28 09:18:31 +0200 (mar. 28 sept. 2010) $
    052     * @since 2.0
    053     */
    054    public class EigenDecompositionImpl implements EigenDecomposition {
    055
    056        /** Maximum number of iterations accepted in the implicit QL transformation */
    057        private byte maxIter = 30;
    058
    059        /** Main diagonal of the tridiagonal matrix. */
    060        private double[] main;
    061
    062        /** Secondary diagonal of the tridiagonal matrix. */
    063        private double[] secondary;
    064
    065        /**
    066         * Transformer to tridiagonal (may be null if matrix is already
    067         * tridiagonal).
    068         */
    069        private TriDiagonalTransformer transformer;
    070
    071        /** Real part of the realEigenvalues. */
    072        private double[] realEigenvalues;
    073
    074        /** Imaginary part of the realEigenvalues. */
    075        private double[] imagEigenvalues;
    076
    077        /** Eigenvectors. */
    078        private ArrayRealVector[] eigenvectors;
    079
    080        /** Cached value of V. */
    081        private RealMatrix cachedV;
    082
    083        /** Cached value of D. */
    084        private RealMatrix cachedD;
    085
    086        /** Cached value of Vt. */
    087        private RealMatrix cachedVt;
    088
    089        /**
    090         * Calculates the eigen decomposition of the given symmetric matrix.
    091         * @param matrix The <strong>symmetric</strong> matrix to decompose.
    092         * @param splitTolerance dummy parameter, present for backward compatibility only.
    093         * @exception InvalidMatrixException (wrapping a
    094         * {@link org.apache.commons.math.ConvergenceException} if algorithm
    095         * fails to converge
    096         */
    097        public EigenDecompositionImpl(final RealMatrix matrix,final double splitTolerance)
    098                throws InvalidMatrixException {
    099            if (isSymmetric(matrix)) {
    100                transformToTridiagonal(matrix);
    101                findEigenVectors(transformer.getQ().getData());
    102            } else {
    103                // as of 2.0, non-symmetric matrices (i.e. complex eigenvalues) are
    104                // NOT supported
    105                // see issue https://issues.apache.org/jira/browse/MATH-235
    106                throw new InvalidMatrixException(
    107                        LocalizedFormats.ASSYMETRIC_EIGEN_NOT_SUPPORTED);
    108            }
    109        }
    110
    111        /**
    112         * Calculates the eigen decomposition of the symmetric tridiagonal
    113         * matrix.  The Householder matrix is assumed to be the identity matrix.
    114         * @param main Main diagonal of the symmetric triadiagonal form
    115         * @param secondary Secondary of the tridiagonal form
    116         * @param splitTolerance dummy parameter, present for backward compatibility only.
    117         * @exception InvalidMatrixException (wrapping a
    118         * {@link org.apache.commons.math.ConvergenceException} if algorithm
    119         * fails to converge
    120         */
    121        public EigenDecompositionImpl(final double[] main,final double[] secondary,
    122                final double splitTolerance)
    123                throws InvalidMatrixException {
    124            this.main      = main.clone();
    125            this.secondary = secondary.clone();
    126            transformer    = null;
    127            final int size=main.length;
    128            double[][] z = new double[size][size];
    129            for (int i=0;i<size;i++) {
    130                z[i][i]=1.0;
    131            }
    132            findEigenVectors(z);
    133        }
    134
    135        /**
    136         * Check if a matrix is symmetric.
    137         * @param matrix
    138         *            matrix to check
    139         * @return true if matrix is symmetric
    140         */
    141        private boolean isSymmetric(final RealMatrix matrix) {
    142            final int rows = matrix.getRowDimension();
    143            final int columns = matrix.getColumnDimension();
    144            final double eps = 10 * rows * columns * MathUtils.EPSILON;
    145            for (int i = 0; i < rows; ++i) {
    146                for (int j = i + 1; j < columns; ++j) {
    147                    final double mij = matrix.getEntry(i, j);
    148                    final double mji = matrix.getEntry(j, i);
    149                    if (FastMath.abs(mij - mji) >
    150                        (FastMath.max(FastMath.abs(mij), FastMath.abs(mji)) * eps)) {
    151                        return false;
    152                    }
    153                }
    154            }
    155            return true;
    156        }
    157
    158        /** {@inheritDoc} */
    159        public RealMatrix getV() throws InvalidMatrixException {
    160
    161            if (cachedV == null) {
    162                final int m = eigenvectors.length;
    163                cachedV = MatrixUtils.createRealMatrix(m, m);
    164                for (int k = 0; k < m; ++k) {
    165                    cachedV.setColumnVector(k, eigenvectors[k]);
    166                }
    167            }
    168            // return the cached matrix
    169            return cachedV;
    170
    171        }
    172
    173        /** {@inheritDoc} */
    174        public RealMatrix getD() throws InvalidMatrixException {
    175            if (cachedD == null) {
    176                // cache the matrix for subsequent calls
    177                cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues);
    178            }
    179            return cachedD;
    180        }
    181
    182        /** {@inheritDoc} */
    183        public RealMatrix getVT() throws InvalidMatrixException {
    184
    185            if (cachedVt == null) {
    186                final int m = eigenvectors.length;
    187                cachedVt = MatrixUtils.createRealMatrix(m, m);
    188                for (int k = 0; k < m; ++k) {
    189                    cachedVt.setRowVector(k, eigenvectors[k]);
    190                }
    191
    192            }
    193
    194            // return the cached matrix
    195            return cachedVt;
    196        }
    197
    198        /** {@inheritDoc} */
    199        public double[] getRealEigenvalues() throws InvalidMatrixException {
    200            return realEigenvalues.clone();
    201        }
    202
    203        /** {@inheritDoc} */
    204        public double getRealEigenvalue(final int i) throws InvalidMatrixException,
    205                ArrayIndexOutOfBoundsException {
    206            return realEigenvalues[i];
    207        }
    208
    209        /** {@inheritDoc} */
    210        public double[] getImagEigenvalues() throws InvalidMatrixException {
    211            return imagEigenvalues.clone();
    212        }
    213
    214        /** {@inheritDoc} */
    215        public double getImagEigenvalue(final int i) throws InvalidMatrixException,
    216                ArrayIndexOutOfBoundsException {
    217            return imagEigenvalues[i];
    218        }
    219
    220        /** {@inheritDoc} */
    221        public RealVector getEigenvector(final int i)
    222                throws InvalidMatrixException, ArrayIndexOutOfBoundsException {
    223            return eigenvectors[i].copy();
    224        }
    225
    226        /**
    227         * Return the determinant of the matrix
    228         * @return determinant of the matrix
    229         */
    230        public double getDeterminant() {
    231            double determinant = 1;
    232            for (double lambda : realEigenvalues) {
    233                determinant *= lambda;
    234            }
    235            return determinant;
    236        }
    237
    238        /** {@inheritDoc} */
    239        public DecompositionSolver getSolver() {
    240            return new Solver(realEigenvalues, imagEigenvalues, eigenvectors);
    241        }
    242
    243        /** Specialized solver. */
    244        private static class Solver implements DecompositionSolver {
    245
    246            /** Real part of the realEigenvalues. */
    247            private double[] realEigenvalues;
    248
    249            /** Imaginary part of the realEigenvalues. */
    250            private double[] imagEigenvalues;
    251
    252            /** Eigenvectors. */
    253            private final ArrayRealVector[] eigenvectors;
    254
    255            /**
    256             * Build a solver from decomposed matrix.
    257             * @param realEigenvalues
    258             *            real parts of the eigenvalues
    259             * @param imagEigenvalues
    260             *            imaginary parts of the eigenvalues
    261             * @param eigenvectors
    262             *            eigenvectors
    263             */
    264            private Solver(final double[] realEigenvalues,
    265                    final double[] imagEigenvalues,
    266                    final ArrayRealVector[] eigenvectors) {
    267                this.realEigenvalues = realEigenvalues;
    268                this.imagEigenvalues = imagEigenvalues;
    269                this.eigenvectors = eigenvectors;
    270            }
    271
    272            /**
    273             * Solve the linear equation A &times; X = B for symmetric matrices A.
    274             * <p>
    275             * This method only find exact linear solutions, i.e. solutions for
    276             * which ||A &times; X - B|| is exactly 0.
    277             * </p>
    278             * @param b
    279             *            right-hand side of the equation A &times; X = B
    280             * @return a vector X that minimizes the two norm of A &times; X - B
    281             * @exception IllegalArgumentException
    282             *                if matrices dimensions don't match
    283             * @exception InvalidMatrixException
    284             *                if decomposed matrix is singular
    285             */
    286            public double[] solve(final double[] b)
    287                    throws IllegalArgumentException, InvalidMatrixException {
    288
    289                if (!isNonSingular()) {
    290                    throw new SingularMatrixException();
    291                }
    292
    293                final int m = realEigenvalues.length;
    294                if (b.length != m) {
    295                    throw MathRuntimeException.createIllegalArgumentException(
    296                            LocalizedFormats.VECTOR_LENGTH_MISMATCH,
    297                            b.length, m);
    298                }
    299
    300                final double[] bp = new double[m];
    301                for (int i = 0; i < m; ++i) {
    302                    final ArrayRealVector v = eigenvectors[i];
    303                    final double[] vData = v.getDataRef();
    304                    final double s = v.dotProduct(b) / realEigenvalues[i];
    305                    for (int j = 0; j < m; ++j) {
    306                        bp[j] += s * vData[j];
    307                    }
    308                }
    309
    310                return bp;
    311
    312            }
    313
    314            /**
    315             * Solve the linear equation A &times; X = B for symmetric matrices A.
    316             * <p>
    317             * This method only find exact linear solutions, i.e. solutions for
    318             * which ||A &times; X - B|| is exactly 0.
    319             * </p>
    320             * @param b
    321             *            right-hand side of the equation A &times; X = B
    322             * @return a vector X that minimizes the two norm of A &times; X - B
    323             * @exception IllegalArgumentException
    324             *                if matrices dimensions don't match
    325             * @exception InvalidMatrixException
    326             *                if decomposed matrix is singular
    327             */
    328            public RealVector solve(final RealVector b)
    329                    throws IllegalArgumentException, InvalidMatrixException {
    330
    331                if (!isNonSingular()) {
    332                    throw new SingularMatrixException();
    333                }
    334
    335                final int m = realEigenvalues.length;
    336                if (b.getDimension() != m) {
    337                    throw MathRuntimeException.createIllegalArgumentException(
    338                            LocalizedFormats.VECTOR_LENGTH_MISMATCH, b
    339                                    .getDimension(), m);
    340                }
    341
    342                final double[] bp = new double[m];
    343                for (int i = 0; i < m; ++i) {
    344                    final ArrayRealVector v = eigenvectors[i];
    345                    final double[] vData = v.getDataRef();
    346                    final double s = v.dotProduct(b) / realEigenvalues[i];
    347                    for (int j = 0; j < m; ++j) {
    348                        bp[j] += s * vData[j];
    349                    }
    350                }
    351
    352                return new ArrayRealVector(bp, false);
    353
    354            }
    355
    356            /**
    357             * Solve the linear equation A &times; X = B for symmetric matrices A.
    358             * <p>
    359             * This method only find exact linear solutions, i.e. solutions for
    360             * which ||A &times; X - B|| is exactly 0.
    361             * </p>
    362             * @param b
    363             *            right-hand side of the equation A &times; X = B
    364             * @return a matrix X that minimizes the two norm of A &times; X - B
    365             * @exception IllegalArgumentException
    366             *                if matrices dimensions don't match
    367             * @exception InvalidMatrixException
    368             *                if decomposed matrix is singular
    369             */
    370            public RealMatrix solve(final RealMatrix b)
    371                    throws IllegalArgumentException, InvalidMatrixException {
    372
    373                if (!isNonSingular()) {
    374                    throw new SingularMatrixException();
    375                }
    376
    377                final int m = realEigenvalues.length;
    378                if (b.getRowDimension() != m) {
    379                    throw MathRuntimeException
    380                            .createIllegalArgumentException(
    381                                    LocalizedFormats.DIMENSIONS_MISMATCH_2x2,
    382                                    b.getRowDimension(), b.getColumnDimension(), m,
    383                                    "n");
    384                }
    385
    386                final int nColB = b.getColumnDimension();
    387                final double[][] bp = new double[m][nColB];
    388                for (int k = 0; k < nColB; ++k) {
    389                    for (int i = 0; i < m; ++i) {
    390                        final ArrayRealVector v = eigenvectors[i];
    391                        final double[] vData = v.getDataRef();
    392                        double s = 0;
    393                        for (int j = 0; j < m; ++j) {
    394                            s += v.getEntry(j) * b.getEntry(j, k);
    395                        }
    396                        s /= realEigenvalues[i];
    397                        for (int j = 0; j < m; ++j) {
    398                            bp[j][k] += s * vData[j];
    399                        }
    400                    }
    401                }
    402
    403                return MatrixUtils.createRealMatrix(bp);
    404
    405            }
    406
    407            /**
    408             * Check if the decomposed matrix is non-singular.
    409             * @return true if the decomposed matrix is non-singular
    410             */
    411            public boolean isNonSingular() {
    412                for (int i = 0; i < realEigenvalues.length; ++i) {
    413                    if ((realEigenvalues[i] == 0) && (imagEigenvalues[i] == 0)) {
    414                        return false;
    415                    }
    416                }
    417                return true;
    418            }
    419
    420            /**
    421             * Get the inverse of the decomposed matrix.
    422             * @return inverse matrix
    423             * @throws InvalidMatrixException
    424             *             if decomposed matrix is singular
    425             */
    426            public RealMatrix getInverse() throws InvalidMatrixException {
    427
    428                if (!isNonSingular()) {
    429                    throw new SingularMatrixException();
    430                }
    431
    432                final int m = realEigenvalues.length;
    433                final double[][] invData = new double[m][m];
    434
    435                for (int i = 0; i < m; ++i) {
    436                    final double[] invI = invData[i];
    437                    for (int j = 0; j < m; ++j) {
    438                        double invIJ = 0;
    439                        for (int k = 0; k < m; ++k) {
    440                            final double[] vK = eigenvectors[k].getDataRef();
    441                            invIJ += vK[i] * vK[j] / realEigenvalues[k];
    442                        }
    443                        invI[j] = invIJ;
    444                    }
    445                }
    446                return MatrixUtils.createRealMatrix(invData);
    447
    448            }
    449
    450        }
    451
    452        /**
    453         * Transform matrix to tridiagonal.
    454         * @param matrix
    455         *            matrix to transform
    456         */
    457        private void transformToTridiagonal(final RealMatrix matrix) {
    458
    459            // transform the matrix to tridiagonal
    460            transformer = new TriDiagonalTransformer(matrix);
    461            main = transformer.getMainDiagonalRef();
    462            secondary = transformer.getSecondaryDiagonalRef();
    463
    464        }
    465
    466        /**
    467         * Find eigenvalues and eigenvectors (Dubrulle et al., 1971)
    468         * @param householderMatrix Householder matrix of the transformation
    469         *  to tri-diagonal form.
    470         */
    471        private void findEigenVectors(double[][] householderMatrix) {
    472
    473            double[][]z = householderMatrix.clone();
    474            final int n = main.length;
    475            realEigenvalues = new double[n];
    476            imagEigenvalues = new double[n];
    477            double[] e = new double[n];
    478            for (int i = 0; i < n - 1; i++) {
    479                realEigenvalues[i] = main[i];
    480                e[i] = secondary[i];
    481            }
    482            realEigenvalues[n - 1] = main[n - 1];
    483            e[n - 1] = 0.0;
    484
    485            // Determine the largest main and secondary value in absolute term.
    486            double maxAbsoluteValue=0.0;
    487            for (int i = 0; i < n; i++) {
    488                if (FastMath.abs(realEigenvalues[i])>maxAbsoluteValue) {
    489                    maxAbsoluteValue=FastMath.abs(realEigenvalues[i]);
    490                }
    491                if (FastMath.abs(e[i])>maxAbsoluteValue) {
    492                    maxAbsoluteValue=FastMath.abs(e[i]);
    493                }
    494            }
    495            // Make null any main and secondary value too small to be significant
    496            if (maxAbsoluteValue!=0.0) {
    497                for (int i=0; i < n; i++) {
    498                    if (FastMath.abs(realEigenvalues[i])<=MathUtils.EPSILON*maxAbsoluteValue) {
    499                        realEigenvalues[i]=0.0;
    500                    }
    501                    if (FastMath.abs(e[i])<=MathUtils.EPSILON*maxAbsoluteValue) {
    502                        e[i]=0.0;
    503                    }
    504                }
    505            }
    506
    507            for (int j = 0; j < n; j++) {
    508                int its = 0;
    509                int m;
    510                do {
    511                    for (m = j; m < n - 1; m++) {
    512                        double delta = FastMath.abs(realEigenvalues[m]) + FastMath.abs(realEigenvalues[m + 1]);
    513                        if (FastMath.abs(e[m]) + delta == delta) {
    514                            break;
    515                        }
    516                    }
    517                    if (m != j) {
    518                        if (its == maxIter)
    519                            throw new InvalidMatrixException(
    520                                    new MaxIterationsExceededException(maxIter));
    521                        its++;
    522                        double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]);
    523                        double t = FastMath.sqrt(1 + q * q);
    524                        if (q < 0.0) {
    525                            q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q - t);
    526                        } else {
    527                            q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q + t);
    528                        }
    529                        double u = 0.0;
    530                        double s = 1.0;
    531                        double c = 1.0;
    532                        int i;
    533                        for (i = m - 1; i >= j; i--) {
    534                            double p = s * e[i];
    535                            double h = c * e[i];
    536                            if (FastMath.abs(p) >= FastMath.abs(q)) {
    537                                c = q / p;
    538                                t = FastMath.sqrt(c * c + 1.0);
    539                                e[i + 1] = p * t;
    540                                s = 1.0 / t;
    541                                c = c * s;
    542                            } else {
    543                                s = p / q;
    544                                t = FastMath.sqrt(s * s + 1.0);
    545                                e[i + 1] = q * t;
    546                                c = 1.0 / t;
    547                                s = s * c;
    548                            }
    549                            if (e[i + 1] == 0.0) {
    550                                realEigenvalues[i + 1] -= u;
    551                                e[m] = 0.0;
    552                                break;
    553                            }
    554                            q = realEigenvalues[i + 1] - u;
    555                            t = (realEigenvalues[i] - q) * s + 2.0 * c * h;
    556                            u = s * t;
    557                            realEigenvalues[i + 1] = q + u;
    558                            q = c * t - h;
    559                            for (int ia = 0; ia < n; ia++) {
    560                                p = z[ia][i + 1];
    561                                z[ia][i + 1] = s * z[ia][i] + c * p;
    562                                z[ia][i] = c * z[ia][i] - s * p;
    563                            }
    564                        }
    565                        if (t == 0.0 && i >= j)
    566                            continue;
    567                        realEigenvalues[j] -= u;
    568                        e[j] = q;
    569                        e[m] = 0.0;
    570                    }
    571                } while (m != j);
    572            }
    573
    574            //Sort the eigen values (and vectors) in increase order
    575            for (int i = 0; i < n; i++) {
    576                int k = i;
    577                double p = realEigenvalues[i];
    578                for (int j = i + 1; j < n; j++) {
    579                    if (realEigenvalues[j] > p) {
    580                        k = j;
    581                        p = realEigenvalues[j];
    582                    }
    583                }
    584                if (k != i) {
    585                    realEigenvalues[k] = realEigenvalues[i];
    586                    realEigenvalues[i] = p;
    587                    for (int j = 0; j < n; j++) {
    588                        p = z[j][i];
    589                        z[j][i] = z[j][k];
    590                        z[j][k] = p;
    591                    }
    592                }
    593            }
    594
    595            // Determine the largest eigen value in absolute term.
    596            maxAbsoluteValue=0.0;
    597            for (int i = 0; i < n; i++) {
    598                if (FastMath.abs(realEigenvalues[i])>maxAbsoluteValue) {
    599                    maxAbsoluteValue=FastMath.abs(realEigenvalues[i]);
    600                }
    601            }
    602            // Make null any eigen value too small to be significant
    603            if (maxAbsoluteValue!=0.0) {
    604                for (int i=0; i < n; i++) {
    605                    if (FastMath.abs(realEigenvalues[i])<MathUtils.EPSILON*maxAbsoluteValue) {
    606                        realEigenvalues[i]=0.0;
    607                    }
    608                }
    609            }
    610            eigenvectors = new ArrayRealVector[n];
    611            double[] tmp = new double[n];
    612            for (int i = 0; i < n; i++) {
    613                for (int j = 0; j < n; j++) {
    614                    tmp[j] = z[j][i];
    615                }
    616                eigenvectors[i] = new ArrayRealVector(tmp);
    617            }
    618        }
    619    }