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