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 × 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 × X = B for symmetric matrices A.
272 * <p>
273 * This method only find exact linear solutions, i.e. solutions for
274 * which ||A × X - B|| is exactly 0.
275 * </p>
276 * @param b
277 * right-hand side of the equation A × X = B
278 * @return a vector X that minimizes the two norm of A × 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 × X = B for symmetric matrices A.
314 * <p>
315 * This method only find exact linear solutions, i.e. solutions for
316 * which ||A × X - B|| is exactly 0.
317 * </p>
318 * @param b
319 * right-hand side of the equation A × X = B
320 * @return a vector X that minimizes the two norm of A × 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 × X = B for symmetric matrices A.
356 * <p>
357 * This method only find exact linear solutions, i.e. solutions for
358 * which ||A × X - B|| is exactly 0.
359 * </p>
360 * @param b
361 * right-hand side of the equation A × X = B
362 * @return a matrix X that minimizes the two norm of A × 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 }