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    }