1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 package org.apache.commons.math.linear;
19
20 import org.apache.commons.math.exception.MaxCountExceededException;
21 import org.apache.commons.math.exception.DimensionMismatchException;
22 import org.apache.commons.math.exception.util.LocalizedFormats;
23 import org.apache.commons.math.util.MathUtils;
24 import org.apache.commons.math.util.FastMath;
25
26 /**
27 * Calculates the eigen decomposition of a real <strong>symmetric</strong>
28 * matrix.
29 * <p>The eigen decomposition of matrix A is a set of two matrices:
30 * V and D such that A = V × D × V<sup>T</sup>.
31 * A, V and D are all m × m matrices.</p>
32 * <p>This class is similar in spirit to the <code>EigenvalueDecomposition</code>
33 * class from the <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a>
34 * library, with the following changes:</p>
35 * <ul>
36 * <li>a {@link #getVT() getVt} method has been added,</li>
37 * <li>two {@link #getRealEigenvalue(int) getRealEigenvalue} and {@link #getImagEigenvalue(int)
38 * getImagEigenvalue} methods to pick up a single eigenvalue have been added,</li>
39 * <li>a {@link #getEigenvector(int) getEigenvector} method to pick up a single
40 * eigenvector has been added,</li>
41 * <li>a {@link #getDeterminant() getDeterminant} method has been added.</li>
42 * <li>a {@link #getSolver() getSolver} method has been added.</li>
43 * </ul>
44 * <p>
45 * As of 2.0, this class supports only <strong>symmetric</strong> matrices, and
46 * hence computes only real realEigenvalues. This implies the D matrix returned
47 * by {@link #getD()} is always diagonal and the imaginary values returned
48 * {@link #getImagEigenvalue(int)} and {@link #getImagEigenvalues()} are always
49 * null.
50 * </p>
51 * <p>
52 * When called with a {@link RealMatrix} argument, this implementation only uses
53 * the upper part of the matrix, the part below the diagonal is not accessed at
54 * all.
55 * </p>
56 * <p>
57 * This implementation is based on the paper by A. Drubrulle, R.S. Martin and
58 * J.H. Wilkinson 'The Implicit QL Algorithm' in Wilksinson and Reinsch (1971)
59 * Handbook for automatic computation, vol. 2, Linear algebra, Springer-Verlag,
60 * New-York
61 * </p>
62 * @see <a href="http://mathworld.wolfram.com/EigenDecomposition.html">MathWorld</a>
63 * @see <a href="http://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix">Wikipedia</a>
64 * @version $Id: EigenDecomposition.java 1177938 2011-10-01 07:57:19Z celestin $
65 * @since 2.0 (changed to concrete class in 3.0)
66 */
67 public class EigenDecomposition{
68
69 /** Maximum number of iterations accepted in the implicit QL transformation */
70 private byte maxIter = 30;
71
72 /** Main diagonal of the tridiagonal matrix. */
73 private double[] main;
74
75 /** Secondary diagonal of the tridiagonal matrix. */
76 private double[] secondary;
77
78 /**
79 * Transformer to tridiagonal (may be null if matrix is already
80 * tridiagonal).
81 */
82 private TriDiagonalTransformer transformer;
83
84 /** Real part of the realEigenvalues. */
85 private double[] realEigenvalues;
86
87 /** Imaginary part of the realEigenvalues. */
88 private double[] imagEigenvalues;
89
90 /** Eigenvectors. */
91 private ArrayRealVector[] eigenvectors;
92
93 /** Cached value of V. */
94 private RealMatrix cachedV;
95
96 /** Cached value of D. */
97 private RealMatrix cachedD;
98
99 /** Cached value of Vt. */
100 private RealMatrix cachedVt;
101
102 /**
103 * Calculates the eigen decomposition of the given symmetric matrix.
104 *
105 * @param matrix Matrix to decompose. It <em>must</em> be symmetric.
106 * @param splitTolerance Dummy parameter (present for backward
107 * compatibility only).
108 * @throws NonSymmetricMatrixException if the matrix is not symmetric.
109 * @throws MaxCountExceededException if the algorithm fails to converge.
110 */
111 public EigenDecomposition(final RealMatrix matrix,
112 final double splitTolerance) {
113 if (isSymmetric(matrix, true)) {
114 transformToTridiagonal(matrix);
115 findEigenVectors(transformer.getQ().getData());
116 }
117 }
118
119 /**
120 * Calculates the eigen decomposition of the symmetric tridiagonal
121 * matrix. The Householder matrix is assumed to be the identity matrix.
122 *
123 * @param main Main diagonal of the symmetric triadiagonal form
124 * @param secondary Secondary of the tridiagonal form
125 * @param splitTolerance Dummy parameter (present for backward
126 * compatibility only).
127 * @throws MaxCountExceededException if the algorithm fails to converge.
128 */
129 public EigenDecomposition(final double[] main,final double[] secondary,
130 final double splitTolerance) {
131 this.main = main.clone();
132 this.secondary = secondary.clone();
133 transformer = null;
134 final int size=main.length;
135 double[][] z = new double[size][size];
136 for (int i=0;i<size;i++) {
137 z[i][i]=1.0;
138 }
139 findEigenVectors(z);
140 }
141
142 /**
143 * Check if a matrix is symmetric.
144 *
145 * @param matrix Matrix to check.
146 * @param raiseException If {@code true}, the method will throw an
147 * exception if {@code matrix} is not symmetric.
148 * @return {@code true} if {@code matrix} is symmetric.
149 * @throws NonSymmetricMatrixException if the matrix is not symmetric and
150 * {@code raiseException} is {@code true}.
151 */
152 private boolean isSymmetric(final RealMatrix matrix,
153 boolean raiseException) {
154 final int rows = matrix.getRowDimension();
155 final int columns = matrix.getColumnDimension();
156 final double eps = 10 * rows * columns * MathUtils.EPSILON;
157 for (int i = 0; i < rows; ++i) {
158 for (int j = i + 1; j < columns; ++j) {
159 final double mij = matrix.getEntry(i, j);
160 final double mji = matrix.getEntry(j, i);
161 if (FastMath.abs(mij - mji) >
162 (FastMath.max(FastMath.abs(mij), FastMath.abs(mji)) * eps)) {
163 if (raiseException) {
164 throw new NonSymmetricMatrixException(i, j, eps);
165 }
166 return false;
167 }
168 }
169 }
170 return true;
171 }
172
173 /**
174 * Returns the matrix V of the decomposition.
175 * <p>V is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
176 * <p>The columns of V are the eigenvectors of the original matrix.</p>
177 * <p>No assumption is made about the orientation of the system axes formed
178 * by the columns of V (e.g. in a 3-dimension space, V can form a left-
179 * or right-handed system).</p>
180 * @return the V matrix
181 */
182 public RealMatrix getV() {
183
184 if (cachedV == null) {
185 final int m = eigenvectors.length;
186 cachedV = MatrixUtils.createRealMatrix(m, m);
187 for (int k = 0; k < m; ++k) {
188 cachedV.setColumnVector(k, eigenvectors[k]);
189 }
190 }
191 // return the cached matrix
192 return cachedV;
193
194 }
195
196 /**
197 * Returns the block diagonal matrix D of the decomposition.
198 * <p>D is a block diagonal matrix.</p>
199 * <p>Real eigenvalues are on the diagonal while complex values are on
200 * 2x2 blocks { {real +imaginary}, {-imaginary, real} }.</p>
201 * @return the D matrix
202 * @see #getRealEigenvalues()
203 * @see #getImagEigenvalues()
204 */
205 public RealMatrix getD() {
206 if (cachedD == null) {
207 // cache the matrix for subsequent calls
208 cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues);
209 }
210 return cachedD;
211 }
212
213 /**
214 * Returns the transpose of the matrix V of the decomposition.
215 * <p>V is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
216 * <p>The columns of V are the eigenvectors of the original matrix.</p>
217 * <p>No assumption is made about the orientation of the system axes formed
218 * by the columns of V (e.g. in a 3-dimension space, V can form a left-
219 * or right-handed system).</p>
220 * @return the transpose of the V matrix
221 */
222 public RealMatrix getVT() {
223
224 if (cachedVt == null) {
225 final int m = eigenvectors.length;
226 cachedVt = MatrixUtils.createRealMatrix(m, m);
227 for (int k = 0; k < m; ++k) {
228 cachedVt.setRowVector(k, eigenvectors[k]);
229 }
230
231 }
232
233 // return the cached matrix
234 return cachedVt;
235 }
236
237 /**
238 * Returns a copy of the real parts of the eigenvalues of the original matrix.
239 * @return a copy of the real parts of the eigenvalues of the original matrix
240 * @see #getD()
241 * @see #getRealEigenvalue(int)
242 * @see #getImagEigenvalues()
243 */
244 public double[] getRealEigenvalues() {
245 return realEigenvalues.clone();
246 }
247
248 /**
249 * Returns the real part of the i<sup>th</sup> eigenvalue of the original matrix.
250 * @param i index of the eigenvalue (counting from 0)
251 * @return real part of the i<sup>th</sup> eigenvalue of the original matrix
252 * @see #getD()
253 * @see #getRealEigenvalues()
254 * @see #getImagEigenvalue(int)
255 */
256 public double getRealEigenvalue(final int i) {
257 return realEigenvalues[i];
258 }
259
260 /**
261 * Returns a copy of the imaginary parts of the eigenvalues of the original matrix.
262 * @return a copy of the imaginary parts of the eigenvalues of the original matrix
263 * @see #getD()
264 * @see #getImagEigenvalue(int)
265 * @see #getRealEigenvalues()
266 */
267 public double[] getImagEigenvalues() {
268 return imagEigenvalues.clone();
269 }
270
271 /**
272 * Returns the imaginary part of the i<sup>th</sup> eigenvalue of the original matrix.
273 * @param i index of the eigenvalue (counting from 0)
274 * @return imaginary part of the i<sup>th</sup> eigenvalue of the original matrix
275 * @see #getD()
276 * @see #getImagEigenvalues()
277 * @see #getRealEigenvalue(int)
278 */
279 public double getImagEigenvalue(final int i) {
280 return imagEigenvalues[i];
281 }
282
283 /**
284 * Returns a copy of the i<sup>th</sup> eigenvector of the original matrix.
285 * @param i index of the eigenvector (counting from 0)
286 * @return copy of the i<sup>th</sup> eigenvector of the original matrix
287 * @see #getD()
288 */
289 public RealVector getEigenvector(final int i) {
290 return eigenvectors[i].copy();
291 }
292
293 /**
294 * Return the determinant of the matrix
295 * @return determinant of the matrix
296 */
297 public double getDeterminant() {
298 double determinant = 1;
299 for (double lambda : realEigenvalues) {
300 determinant *= lambda;
301 }
302 return determinant;
303 }
304
305 /**
306 * Get a solver for finding the A × X = B solution in exact linear sense.
307 * @return a solver
308 */
309 public DecompositionSolver getSolver() {
310 return new Solver(realEigenvalues, imagEigenvalues, eigenvectors);
311 }
312
313 /** Specialized solver. */
314 private static class Solver implements DecompositionSolver {
315
316 /** Real part of the realEigenvalues. */
317 private double[] realEigenvalues;
318
319 /** Imaginary part of the realEigenvalues. */
320 private double[] imagEigenvalues;
321
322 /** Eigenvectors. */
323 private final ArrayRealVector[] eigenvectors;
324
325 /**
326 * Build a solver from decomposed matrix.
327 * @param realEigenvalues
328 * real parts of the eigenvalues
329 * @param imagEigenvalues
330 * imaginary parts of the eigenvalues
331 * @param eigenvectors
332 * eigenvectors
333 */
334 private Solver(final double[] realEigenvalues,
335 final double[] imagEigenvalues,
336 final ArrayRealVector[] eigenvectors) {
337 this.realEigenvalues = realEigenvalues;
338 this.imagEigenvalues = imagEigenvalues;
339 this.eigenvectors = eigenvectors;
340 }
341
342 /**
343 * Solve the linear equation A × X = B for symmetric matrices A.
344 * <p>
345 * This method only find exact linear solutions, i.e. solutions for
346 * which ||A × X - B|| is exactly 0.
347 * </p>
348 * @param b Right-hand side of the equation A × X = B
349 * @return a Vector X that minimizes the two norm of A × X - B
350 * @throws DimensionMismatchException if the matrices dimensions do not match.
351 * @throws SingularMatrixException if the decomposed matrix is singular.
352 */
353 public RealVector solve(final RealVector b) {
354 if (!isNonSingular()) {
355 throw new SingularMatrixException();
356 }
357
358 final int m = realEigenvalues.length;
359 if (b.getDimension() != m) {
360 throw new DimensionMismatchException(b.getDimension(), m);
361 }
362
363 final double[] bp = new double[m];
364 for (int i = 0; i < m; ++i) {
365 final ArrayRealVector v = eigenvectors[i];
366 final double[] vData = v.getDataRef();
367 final double s = v.dotProduct(b) / realEigenvalues[i];
368 for (int j = 0; j < m; ++j) {
369 bp[j] += s * vData[j];
370 }
371 }
372
373 return new ArrayRealVector(bp, false);
374 }
375
376 /** {@inheritDoc} */
377 public RealMatrix solve(RealMatrix b) {
378
379 if (!isNonSingular()) {
380 throw new SingularMatrixException();
381 }
382
383 final int m = realEigenvalues.length;
384 if (b.getRowDimension() != m) {
385 throw new DimensionMismatchException(b.getRowDimension(), m);
386 }
387
388 final int nColB = b.getColumnDimension();
389 final double[][] bp = new double[m][nColB];
390 final double[] tmpCol = new double[m];
391 for (int k = 0; k < nColB; ++k) {
392 for (int i = 0; i < m; ++i) {
393 tmpCol[i] = b.getEntry(i, k);
394 bp[i][k] = 0;
395 }
396 for (int i = 0; i < m; ++i) {
397 final ArrayRealVector v = eigenvectors[i];
398 final double[] vData = v.getDataRef();
399 double s = 0;
400 for (int j = 0; j < m; ++j) {
401 s += v.getEntry(j) * tmpCol[j];
402 }
403 s /= realEigenvalues[i];
404 for (int j = 0; j < m; ++j) {
405 bp[j][k] += s * vData[j];
406 }
407 }
408 }
409
410 return new Array2DRowRealMatrix(bp, false);
411
412 }
413
414 /**
415 * Check if the decomposed matrix is non-singular.
416 * @return true if the decomposed matrix is non-singular
417 */
418 public boolean isNonSingular() {
419 for (int i = 0; i < realEigenvalues.length; ++i) {
420 if ((realEigenvalues[i] == 0) && (imagEigenvalues[i] == 0)) {
421 return false;
422 }
423 }
424 return true;
425 }
426
427 /**
428 * Get the inverse of the decomposed matrix.
429 *
430 * @return the inverse matrix.
431 * @throws SingularMatrixException if the decomposed matrix is singular.
432 */
433 public RealMatrix getInverse() {
434 if (!isNonSingular()) {
435 throw new SingularMatrixException();
436 }
437
438 final int m = realEigenvalues.length;
439 final double[][] invData = new double[m][m];
440
441 for (int i = 0; i < m; ++i) {
442 final double[] invI = invData[i];
443 for (int j = 0; j < m; ++j) {
444 double invIJ = 0;
445 for (int k = 0; k < m; ++k) {
446 final double[] vK = eigenvectors[k].getDataRef();
447 invIJ += vK[i] * vK[j] / realEigenvalues[k];
448 }
449 invI[j] = invIJ;
450 }
451 }
452 return MatrixUtils.createRealMatrix(invData);
453 }
454 }
455
456 /**
457 * Transform matrix to tridiagonal.
458 *
459 * @param matrix Matrix to transform.
460 */
461 private void transformToTridiagonal(final RealMatrix matrix) {
462 // transform the matrix to tridiagonal
463 transformer = new TriDiagonalTransformer(matrix);
464 main = transformer.getMainDiagonalRef();
465 secondary = transformer.getSecondaryDiagonalRef();
466 }
467
468 /**
469 * Find eigenvalues and eigenvectors (Dubrulle et al., 1971)
470 *
471 * @param householderMatrix Householder matrix of the transformation
472 * to tri-diagonal form.
473 */
474 private void findEigenVectors(double[][] householderMatrix) {
475 double[][]z = householderMatrix.clone();
476 final int n = main.length;
477 realEigenvalues = new double[n];
478 imagEigenvalues = new double[n];
479 double[] e = new double[n];
480 for (int i = 0; i < n - 1; i++) {
481 realEigenvalues[i] = main[i];
482 e[i] = secondary[i];
483 }
484 realEigenvalues[n - 1] = main[n - 1];
485 e[n - 1] = 0.0;
486
487 // Determine the largest main and secondary value in absolute term.
488 double maxAbsoluteValue=0.0;
489 for (int i = 0; i < n; i++) {
490 if (FastMath.abs(realEigenvalues[i])>maxAbsoluteValue) {
491 maxAbsoluteValue=FastMath.abs(realEigenvalues[i]);
492 }
493 if (FastMath.abs(e[i])>maxAbsoluteValue) {
494 maxAbsoluteValue=FastMath.abs(e[i]);
495 }
496 }
497 // Make null any main and secondary value too small to be significant
498 if (maxAbsoluteValue!=0.0) {
499 for (int i=0; i < n; i++) {
500 if (FastMath.abs(realEigenvalues[i])<=MathUtils.EPSILON*maxAbsoluteValue) {
501 realEigenvalues[i]=0.0;
502 }
503 if (FastMath.abs(e[i])<=MathUtils.EPSILON*maxAbsoluteValue) {
504 e[i]=0.0;
505 }
506 }
507 }
508
509 for (int j = 0; j < n; j++) {
510 int its = 0;
511 int m;
512 do {
513 for (m = j; m < n - 1; m++) {
514 double delta = FastMath.abs(realEigenvalues[m]) + FastMath.abs(realEigenvalues[m + 1]);
515 if (FastMath.abs(e[m]) + delta == delta) {
516 break;
517 }
518 }
519 if (m != j) {
520 if (its == maxIter) {
521 throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED,
522 maxIter);
523 }
524 its++;
525 double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]);
526 double t = FastMath.sqrt(1 + q * q);
527 if (q < 0.0) {
528 q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q - t);
529 } else {
530 q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q + t);
531 }
532 double u = 0.0;
533 double s = 1.0;
534 double c = 1.0;
535 int i;
536 for (i = m - 1; i >= j; i--) {
537 double p = s * e[i];
538 double h = c * e[i];
539 if (FastMath.abs(p) >= FastMath.abs(q)) {
540 c = q / p;
541 t = FastMath.sqrt(c * c + 1.0);
542 e[i + 1] = p * t;
543 s = 1.0 / t;
544 c = c * s;
545 } else {
546 s = p / q;
547 t = FastMath.sqrt(s * s + 1.0);
548 e[i + 1] = q * t;
549 c = 1.0 / t;
550 s = s * c;
551 }
552 if (e[i + 1] == 0.0) {
553 realEigenvalues[i + 1] -= u;
554 e[m] = 0.0;
555 break;
556 }
557 q = realEigenvalues[i + 1] - u;
558 t = (realEigenvalues[i] - q) * s + 2.0 * c * h;
559 u = s * t;
560 realEigenvalues[i + 1] = q + u;
561 q = c * t - h;
562 for (int ia = 0; ia < n; ia++) {
563 p = z[ia][i + 1];
564 z[ia][i + 1] = s * z[ia][i] + c * p;
565 z[ia][i] = c * z[ia][i] - s * p;
566 }
567 }
568 if (t == 0.0 && i >= j) {
569 continue;
570 }
571 realEigenvalues[j] -= u;
572 e[j] = q;
573 e[m] = 0.0;
574 }
575 } while (m != j);
576 }
577
578 //Sort the eigen values (and vectors) in increase order
579 for (int i = 0; i < n; i++) {
580 int k = i;
581 double p = realEigenvalues[i];
582 for (int j = i + 1; j < n; j++) {
583 if (realEigenvalues[j] > p) {
584 k = j;
585 p = realEigenvalues[j];
586 }
587 }
588 if (k != i) {
589 realEigenvalues[k] = realEigenvalues[i];
590 realEigenvalues[i] = p;
591 for (int j = 0; j < n; j++) {
592 p = z[j][i];
593 z[j][i] = z[j][k];
594 z[j][k] = p;
595 }
596 }
597 }
598
599 // Determine the largest eigen value in absolute term.
600 maxAbsoluteValue=0.0;
601 for (int i = 0; i < n; i++) {
602 if (FastMath.abs(realEigenvalues[i])>maxAbsoluteValue) {
603 maxAbsoluteValue=FastMath.abs(realEigenvalues[i]);
604 }
605 }
606 // Make null any eigen value too small to be significant
607 if (maxAbsoluteValue!=0.0) {
608 for (int i=0; i < n; i++) {
609 if (FastMath.abs(realEigenvalues[i])<MathUtils.EPSILON*maxAbsoluteValue) {
610 realEigenvalues[i]=0.0;
611 }
612 }
613 }
614 eigenvectors = new ArrayRealVector[n];
615 double[] tmp = new double[n];
616 for (int i = 0; i < n; i++) {
617 for (int j = 0; j < n; j++) {
618 tmp[j] = z[j][i];
619 }
620 eigenvectors[i] = new ArrayRealVector(tmp);
621 }
622 }
623 }