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.exception.MaxCountExceededException;
021 import org.apache.commons.math.exception.DimensionMismatchException;
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 × 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 $Id: EigenDecompositionImpl.java 1139906 2011-06-26 18:42:32Z luc $
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 *
092 * @param matrix Matrix to decompose. It <em>must</em> be symmetric.
093 * @param splitTolerance Dummy parameter (present for backward
094 * compatibility only).
095 * @throws NonSymmetricMatrixException if the matrix is not symmetric.
096 * @throws MaxCountExceededException if the algorithm fails to converge.
097 */
098 public EigenDecompositionImpl(final RealMatrix matrix,
099 final double splitTolerance) {
100 if (isSymmetric(matrix, true)) {
101 transformToTridiagonal(matrix);
102 findEigenVectors(transformer.getQ().getData());
103 }
104 }
105
106 /**
107 * Calculates the eigen decomposition of the symmetric tridiagonal
108 * matrix. The Householder matrix is assumed to be the identity matrix.
109 *
110 * @param main Main diagonal of the symmetric triadiagonal form
111 * @param secondary Secondary of the tridiagonal form
112 * @param splitTolerance Dummy parameter (present for backward
113 * compatibility only).
114 * @throws MaxCountExceededException if the algorithm fails to converge.
115 */
116 public EigenDecompositionImpl(final double[] main,final double[] secondary,
117 final double splitTolerance) {
118 this.main = main.clone();
119 this.secondary = secondary.clone();
120 transformer = null;
121 final int size=main.length;
122 double[][] z = new double[size][size];
123 for (int i=0;i<size;i++) {
124 z[i][i]=1.0;
125 }
126 findEigenVectors(z);
127 }
128
129 /**
130 * Check if a matrix is symmetric.
131 *
132 * @param matrix Matrix to check.
133 * @param raiseException If {@code true}, the method will throw an
134 * exception if {@code matrix} is not symmetric.
135 * @return {@code true} if {@code matrix} is symmetric.
136 * @throws NonSymmetricMatrixException if the matrix is not symmetric and
137 * {@code raiseException} is {@code true}.
138 */
139 private boolean isSymmetric(final RealMatrix matrix,
140 boolean raiseException) {
141 final int rows = matrix.getRowDimension();
142 final int columns = matrix.getColumnDimension();
143 final double eps = 10 * rows * columns * MathUtils.EPSILON;
144 for (int i = 0; i < rows; ++i) {
145 for (int j = i + 1; j < columns; ++j) {
146 final double mij = matrix.getEntry(i, j);
147 final double mji = matrix.getEntry(j, i);
148 if (FastMath.abs(mij - mji) >
149 (FastMath.max(FastMath.abs(mij), FastMath.abs(mji)) * eps)) {
150 if (raiseException) {
151 throw new NonSymmetricMatrixException(i, j, eps);
152 }
153 return false;
154 }
155 }
156 }
157 return true;
158 }
159
160 /** {@inheritDoc} */
161 public RealMatrix getV() {
162
163 if (cachedV == null) {
164 final int m = eigenvectors.length;
165 cachedV = MatrixUtils.createRealMatrix(m, m);
166 for (int k = 0; k < m; ++k) {
167 cachedV.setColumnVector(k, eigenvectors[k]);
168 }
169 }
170 // return the cached matrix
171 return cachedV;
172
173 }
174
175 /** {@inheritDoc} */
176 public RealMatrix getD() {
177 if (cachedD == null) {
178 // cache the matrix for subsequent calls
179 cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues);
180 }
181 return cachedD;
182 }
183
184 /** {@inheritDoc} */
185 public RealMatrix getVT() {
186
187 if (cachedVt == null) {
188 final int m = eigenvectors.length;
189 cachedVt = MatrixUtils.createRealMatrix(m, m);
190 for (int k = 0; k < m; ++k) {
191 cachedVt.setRowVector(k, eigenvectors[k]);
192 }
193
194 }
195
196 // return the cached matrix
197 return cachedVt;
198 }
199
200 /** {@inheritDoc} */
201 public double[] getRealEigenvalues() {
202 return realEigenvalues.clone();
203 }
204
205 /** {@inheritDoc} */
206 public double getRealEigenvalue(final int i) {
207 return realEigenvalues[i];
208 }
209
210 /** {@inheritDoc} */
211 public double[] getImagEigenvalues() {
212 return imagEigenvalues.clone();
213 }
214
215 /** {@inheritDoc} */
216 public double getImagEigenvalue(final int i) {
217 return imagEigenvalues[i];
218 }
219
220 /** {@inheritDoc} */
221 public RealVector getEigenvector(final int i) {
222 return eigenvectors[i].copy();
223 }
224
225 /**
226 * Return the determinant of the matrix
227 * @return determinant of the matrix
228 */
229 public double getDeterminant() {
230 double determinant = 1;
231 for (double lambda : realEigenvalues) {
232 determinant *= lambda;
233 }
234 return determinant;
235 }
236
237 /** {@inheritDoc} */
238 public DecompositionSolver getSolver() {
239 return new Solver(realEigenvalues, imagEigenvalues, eigenvectors);
240 }
241
242 /** Specialized solver. */
243 private static class Solver implements DecompositionSolver {
244
245 /** Real part of the realEigenvalues. */
246 private double[] realEigenvalues;
247
248 /** Imaginary part of the realEigenvalues. */
249 private double[] imagEigenvalues;
250
251 /** Eigenvectors. */
252 private final ArrayRealVector[] eigenvectors;
253
254 /**
255 * Build a solver from decomposed matrix.
256 * @param realEigenvalues
257 * real parts of the eigenvalues
258 * @param imagEigenvalues
259 * imaginary parts of the eigenvalues
260 * @param eigenvectors
261 * eigenvectors
262 */
263 private Solver(final double[] realEigenvalues,
264 final double[] imagEigenvalues,
265 final ArrayRealVector[] eigenvectors) {
266 this.realEigenvalues = realEigenvalues;
267 this.imagEigenvalues = imagEigenvalues;
268 this.eigenvectors = eigenvectors;
269 }
270
271 /**
272 * Solve the linear equation A × X = B for symmetric matrices A.
273 * <p>
274 * This method only find exact linear solutions, i.e. solutions for
275 * which ||A × X - B|| is exactly 0.
276 * </p>
277 * @param b Right-hand side of the equation A × X = B
278 * @return a Vector X that minimizes the two norm of A × X - B
279 * @throws DimensionMismatchException if the matrices dimensions do not match.
280 * @throws SingularMatrixException if the decomposed matrix is singular.
281 */
282 public double[] solve(final double[] b) {
283
284 if (!isNonSingular()) {
285 throw new SingularMatrixException();
286 }
287
288 final int m = realEigenvalues.length;
289 if (b.length != m) {
290 throw new DimensionMismatchException(b.length, m);
291 }
292
293 final double[] bp = new double[m];
294 for (int i = 0; i < m; ++i) {
295 final ArrayRealVector v = eigenvectors[i];
296 final double[] vData = v.getDataRef();
297 final double s = v.dotProduct(b) / realEigenvalues[i];
298 for (int j = 0; j < m; ++j) {
299 bp[j] += s * vData[j];
300 }
301 }
302
303 return bp;
304
305 }
306
307 /**
308 * Solve the linear equation A × X = B for symmetric matrices A.
309 * <p>
310 * This method only find exact linear solutions, i.e. solutions for
311 * which ||A × X - B|| is exactly 0.
312 * </p>
313 * @param b Right-hand side of the equation A × X = B
314 * @return a Vector X that minimizes the two norm of A × X - B
315 * @throws DimensionMismatchException if the matrices dimensions do not match.
316 * @throws SingularMatrixException if the decomposed matrix is singular.
317 */
318 public RealVector solve(final RealVector b) {
319 if (!isNonSingular()) {
320 throw new SingularMatrixException();
321 }
322
323 final int m = realEigenvalues.length;
324 if (b.getDimension() != m) {
325 throw new DimensionMismatchException(b.getDimension(), m);
326 }
327
328 final double[] bp = new double[m];
329 for (int i = 0; i < m; ++i) {
330 final ArrayRealVector v = eigenvectors[i];
331 final double[] vData = v.getDataRef();
332 final double s = v.dotProduct(b) / realEigenvalues[i];
333 for (int j = 0; j < m; ++j) {
334 bp[j] += s * vData[j];
335 }
336 }
337
338 return new ArrayRealVector(bp, false);
339 }
340
341 /** Solve the linear equation A × X = B for matrices A.
342 * <p>The A matrix is implicit, it is provided by the underlying
343 * decomposition algorithm.</p>
344 * @param b right-hand side of the equation A × X = B
345 * @param reuseB if true, the b array will be reused and returned,
346 * instead of being copied
347 * @return a matrix X that minimizes the two norm of A × X - B
348 * @throws org.apache.commons.math.exception.DimensionMismatchException
349 * if the matrices dimensions do not match.
350 * @throws SingularMatrixException
351 * if the decomposed matrix is singular.
352 */
353 private double[][] solve(double[][] b, boolean reuseB) {
354
355 if (!isNonSingular()) {
356 throw new SingularMatrixException();
357 }
358
359 final int m = realEigenvalues.length;
360 if (b.length != m) {
361 throw new DimensionMismatchException(b.length, m);
362 }
363
364 final int nColB = b[0].length;
365 final double[][] bp;
366 if (reuseB) {
367 bp = b;
368 } else {
369 bp = new double[m][nColB];
370 }
371 final double[] tmpCol = new double[m];
372 for (int k = 0; k < nColB; ++k) {
373 for (int i = 0; i < m; ++i) {
374 tmpCol[i] = b[i][k];
375 bp[i][k] = 0;
376 }
377 for (int i = 0; i < m; ++i) {
378 final ArrayRealVector v = eigenvectors[i];
379 final double[] vData = v.getDataRef();
380 double s = 0;
381 for (int j = 0; j < m; ++j) {
382 s += v.getEntry(j) * tmpCol[j];
383 }
384 s /= realEigenvalues[i];
385 for (int j = 0; j < m; ++j) {
386 bp[j][k] += s * vData[j];
387 }
388 }
389 }
390
391 return bp;
392
393 }
394
395 /** {@inheritDoc} */
396 public double[][] solve(double[][] b) {
397 return solve(b, false);
398 }
399
400 /** {@inheritDoc} */
401 public RealMatrix solve(RealMatrix b) {
402 return new Array2DRowRealMatrix(solve(b.getData(), true), false);
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 *
421 * @return the inverse matrix.
422 * @throws SingularMatrixException if the decomposed matrix is singular.
423 */
424 public RealMatrix getInverse() {
425 if (!isNonSingular()) {
426 throw new SingularMatrixException();
427 }
428
429 final int m = realEigenvalues.length;
430 final double[][] invData = new double[m][m];
431
432 for (int i = 0; i < m; ++i) {
433 final double[] invI = invData[i];
434 for (int j = 0; j < m; ++j) {
435 double invIJ = 0;
436 for (int k = 0; k < m; ++k) {
437 final double[] vK = eigenvectors[k].getDataRef();
438 invIJ += vK[i] * vK[j] / realEigenvalues[k];
439 }
440 invI[j] = invIJ;
441 }
442 }
443 return MatrixUtils.createRealMatrix(invData);
444 }
445 }
446
447 /**
448 * Transform matrix to tridiagonal.
449 *
450 * @param matrix Matrix to transform.
451 */
452 private void transformToTridiagonal(final RealMatrix matrix) {
453 // transform the matrix to tridiagonal
454 transformer = new TriDiagonalTransformer(matrix);
455 main = transformer.getMainDiagonalRef();
456 secondary = transformer.getSecondaryDiagonalRef();
457 }
458
459 /**
460 * Find eigenvalues and eigenvectors (Dubrulle et al., 1971)
461 *
462 * @param householderMatrix Householder matrix of the transformation
463 * to tri-diagonal form.
464 */
465 private void findEigenVectors(double[][] householderMatrix) {
466 double[][]z = householderMatrix.clone();
467 final int n = main.length;
468 realEigenvalues = new double[n];
469 imagEigenvalues = new double[n];
470 double[] e = new double[n];
471 for (int i = 0; i < n - 1; i++) {
472 realEigenvalues[i] = main[i];
473 e[i] = secondary[i];
474 }
475 realEigenvalues[n - 1] = main[n - 1];
476 e[n - 1] = 0.0;
477
478 // Determine the largest main and secondary value in absolute term.
479 double maxAbsoluteValue=0.0;
480 for (int i = 0; i < n; i++) {
481 if (FastMath.abs(realEigenvalues[i])>maxAbsoluteValue) {
482 maxAbsoluteValue=FastMath.abs(realEigenvalues[i]);
483 }
484 if (FastMath.abs(e[i])>maxAbsoluteValue) {
485 maxAbsoluteValue=FastMath.abs(e[i]);
486 }
487 }
488 // Make null any main and secondary value too small to be significant
489 if (maxAbsoluteValue!=0.0) {
490 for (int i=0; i < n; i++) {
491 if (FastMath.abs(realEigenvalues[i])<=MathUtils.EPSILON*maxAbsoluteValue) {
492 realEigenvalues[i]=0.0;
493 }
494 if (FastMath.abs(e[i])<=MathUtils.EPSILON*maxAbsoluteValue) {
495 e[i]=0.0;
496 }
497 }
498 }
499
500 for (int j = 0; j < n; j++) {
501 int its = 0;
502 int m;
503 do {
504 for (m = j; m < n - 1; m++) {
505 double delta = FastMath.abs(realEigenvalues[m]) + FastMath.abs(realEigenvalues[m + 1]);
506 if (FastMath.abs(e[m]) + delta == delta) {
507 break;
508 }
509 }
510 if (m != j) {
511 if (its == maxIter) {
512 throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED,
513 maxIter);
514 }
515 its++;
516 double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]);
517 double t = FastMath.sqrt(1 + q * q);
518 if (q < 0.0) {
519 q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q - t);
520 } else {
521 q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q + t);
522 }
523 double u = 0.0;
524 double s = 1.0;
525 double c = 1.0;
526 int i;
527 for (i = m - 1; i >= j; i--) {
528 double p = s * e[i];
529 double h = c * e[i];
530 if (FastMath.abs(p) >= FastMath.abs(q)) {
531 c = q / p;
532 t = FastMath.sqrt(c * c + 1.0);
533 e[i + 1] = p * t;
534 s = 1.0 / t;
535 c = c * s;
536 } else {
537 s = p / q;
538 t = FastMath.sqrt(s * s + 1.0);
539 e[i + 1] = q * t;
540 c = 1.0 / t;
541 s = s * c;
542 }
543 if (e[i + 1] == 0.0) {
544 realEigenvalues[i + 1] -= u;
545 e[m] = 0.0;
546 break;
547 }
548 q = realEigenvalues[i + 1] - u;
549 t = (realEigenvalues[i] - q) * s + 2.0 * c * h;
550 u = s * t;
551 realEigenvalues[i + 1] = q + u;
552 q = c * t - h;
553 for (int ia = 0; ia < n; ia++) {
554 p = z[ia][i + 1];
555 z[ia][i + 1] = s * z[ia][i] + c * p;
556 z[ia][i] = c * z[ia][i] - s * p;
557 }
558 }
559 if (t == 0.0 && i >= j) {
560 continue;
561 }
562 realEigenvalues[j] -= u;
563 e[j] = q;
564 e[m] = 0.0;
565 }
566 } while (m != j);
567 }
568
569 //Sort the eigen values (and vectors) in increase order
570 for (int i = 0; i < n; i++) {
571 int k = i;
572 double p = realEigenvalues[i];
573 for (int j = i + 1; j < n; j++) {
574 if (realEigenvalues[j] > p) {
575 k = j;
576 p = realEigenvalues[j];
577 }
578 }
579 if (k != i) {
580 realEigenvalues[k] = realEigenvalues[i];
581 realEigenvalues[i] = p;
582 for (int j = 0; j < n; j++) {
583 p = z[j][i];
584 z[j][i] = z[j][k];
585 z[j][k] = p;
586 }
587 }
588 }
589
590 // Determine the largest eigen value in absolute term.
591 maxAbsoluteValue=0.0;
592 for (int i = 0; i < n; i++) {
593 if (FastMath.abs(realEigenvalues[i])>maxAbsoluteValue) {
594 maxAbsoluteValue=FastMath.abs(realEigenvalues[i]);
595 }
596 }
597 // Make null any eigen value too small to be significant
598 if (maxAbsoluteValue!=0.0) {
599 for (int i=0; i < n; i++) {
600 if (FastMath.abs(realEigenvalues[i])<MathUtils.EPSILON*maxAbsoluteValue) {
601 realEigenvalues[i]=0.0;
602 }
603 }
604 }
605 eigenvectors = new ArrayRealVector[n];
606 double[] tmp = new double[n];
607 for (int i = 0; i < n; i++) {
608 for (int j = 0; j < n; j++) {
609 tmp[j] = z[j][i];
610 }
611 eigenvectors[i] = new ArrayRealVector(tmp);
612 }
613 }
614 }