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 agents.org.apache.commons.math.linear;
|
---|
19 |
|
---|
20 | import agents.org.apache.commons.math.MathRuntimeException;
|
---|
21 | import agents.org.apache.commons.math.MaxIterationsExceededException;
|
---|
22 | import agents.org.apache.commons.math.exception.util.LocalizedFormats;
|
---|
23 | import agents.org.apache.commons.math.util.FastMath;
|
---|
24 | import agents.org.apache.commons.math.util.MathUtils;
|
---|
25 |
|
---|
26 | /**
|
---|
27 | * Calculates the eigen decomposition of a real <strong>symmetric</strong>
|
---|
28 | * matrix.
|
---|
29 | * <p>
|
---|
30 | * The eigen decomposition of matrix A is a set of two matrices: V and D such
|
---|
31 | * that A = V D V<sup>T</sup>. A, V and D are all m × m matrices.
|
---|
32 | * </p>
|
---|
33 | * <p>
|
---|
34 | * As of 2.0, this class supports only <strong>symmetric</strong> matrices, and
|
---|
35 | * hence computes only real realEigenvalues. This implies the D matrix returned
|
---|
36 | * by {@link #getD()} is always diagonal and the imaginary values returned
|
---|
37 | * {@link #getImagEigenvalue(int)} and {@link #getImagEigenvalues()} are always
|
---|
38 | * null.
|
---|
39 | * </p>
|
---|
40 | * <p>
|
---|
41 | * When called with a {@link RealMatrix} argument, this implementation only uses
|
---|
42 | * the upper part of the matrix, the part below the diagonal is not accessed at
|
---|
43 | * all.
|
---|
44 | * </p>
|
---|
45 | * <p>
|
---|
46 | * This implementation is based on the paper by A. Drubrulle, R.S. Martin and
|
---|
47 | * J.H. Wilkinson 'The Implicit QL Algorithm' in Wilksinson and Reinsch (1971)
|
---|
48 | * Handbook for automatic computation, vol. 2, Linear algebra, Springer-Verlag,
|
---|
49 | * New-York
|
---|
50 | * </p>
|
---|
51 | * @version $Revision: 1002040 $ $Date: 2010-09-28 09:18:31 +0200 (mar. 28 sept. 2010) $
|
---|
52 | * @since 2.0
|
---|
53 | */
|
---|
54 | public class EigenDecompositionImpl implements EigenDecomposition {
|
---|
55 |
|
---|
56 | /** Maximum number of iterations accepted in the implicit QL transformation */
|
---|
57 | private byte maxIter = 30;
|
---|
58 |
|
---|
59 | /** Main diagonal of the tridiagonal matrix. */
|
---|
60 | private double[] main;
|
---|
61 |
|
---|
62 | /** Secondary diagonal of the tridiagonal matrix. */
|
---|
63 | private double[] secondary;
|
---|
64 |
|
---|
65 | /**
|
---|
66 | * Transformer to tridiagonal (may be null if matrix is already
|
---|
67 | * tridiagonal).
|
---|
68 | */
|
---|
69 | private TriDiagonalTransformer transformer;
|
---|
70 |
|
---|
71 | /** Real part of the realEigenvalues. */
|
---|
72 | private double[] realEigenvalues;
|
---|
73 |
|
---|
74 | /** Imaginary part of the realEigenvalues. */
|
---|
75 | private double[] imagEigenvalues;
|
---|
76 |
|
---|
77 | /** Eigenvectors. */
|
---|
78 | private ArrayRealVector[] eigenvectors;
|
---|
79 |
|
---|
80 | /** Cached value of V. */
|
---|
81 | private RealMatrix cachedV;
|
---|
82 |
|
---|
83 | /** Cached value of D. */
|
---|
84 | private RealMatrix cachedD;
|
---|
85 |
|
---|
86 | /** Cached value of Vt. */
|
---|
87 | private RealMatrix cachedVt;
|
---|
88 |
|
---|
89 | /**
|
---|
90 | * Calculates the eigen decomposition of the given symmetric matrix.
|
---|
91 | * @param matrix The <strong>symmetric</strong> matrix to decompose.
|
---|
92 | * @param splitTolerance dummy parameter, present for backward compatibility only.
|
---|
93 | * @exception InvalidMatrixException (wrapping a
|
---|
94 | * {@link agents.org.apache.commons.math.ConvergenceException} if algorithm
|
---|
95 | * fails to converge
|
---|
96 | */
|
---|
97 | public EigenDecompositionImpl(final RealMatrix matrix,final double splitTolerance)
|
---|
98 | throws InvalidMatrixException {
|
---|
99 | 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 agents.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 × X = B for symmetric matrices A.
|
---|
274 | * <p>
|
---|
275 | * This method only find exact linear solutions, i.e. solutions for
|
---|
276 | * which ||A × X - B|| is exactly 0.
|
---|
277 | * </p>
|
---|
278 | * @param b
|
---|
279 | * right-hand side of the equation A × X = B
|
---|
280 | * @return a vector X that minimizes the two norm of A × 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 × X = B for symmetric matrices A.
|
---|
316 | * <p>
|
---|
317 | * This method only find exact linear solutions, i.e. solutions for
|
---|
318 | * which ||A × X - B|| is exactly 0.
|
---|
319 | * </p>
|
---|
320 | * @param b
|
---|
321 | * right-hand side of the equation A × X = B
|
---|
322 | * @return a vector X that minimizes the two norm of A × 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 × X = B for symmetric matrices A.
|
---|
358 | * <p>
|
---|
359 | * This method only find exact linear solutions, i.e. solutions for
|
---|
360 | * which ||A × X - B|| is exactly 0.
|
---|
361 | * </p>
|
---|
362 | * @param b
|
---|
363 | * right-hand side of the equation A × X = B
|
---|
364 | * @return a matrix X that minimizes the two norm of A × 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 | }
|
---|