summaryrefslogtreecommitdiff
path: root/src/main/java/org/apache/commons/math3/linear/EigenDecomposition.java
diff options
context:
space:
mode:
Diffstat (limited to 'src/main/java/org/apache/commons/math3/linear/EigenDecomposition.java')
-rw-r--r--src/main/java/org/apache/commons/math3/linear/EigenDecomposition.java968
1 files changed, 968 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/linear/EigenDecomposition.java b/src/main/java/org/apache/commons/math3/linear/EigenDecomposition.java
new file mode 100644
index 0000000..505897f
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/linear/EigenDecomposition.java
@@ -0,0 +1,968 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math3.linear;
+
+import org.apache.commons.math3.complex.Complex;
+import org.apache.commons.math3.exception.DimensionMismatchException;
+import org.apache.commons.math3.exception.MathArithmeticException;
+import org.apache.commons.math3.exception.MathUnsupportedOperationException;
+import org.apache.commons.math3.exception.MaxCountExceededException;
+import org.apache.commons.math3.exception.util.LocalizedFormats;
+import org.apache.commons.math3.util.FastMath;
+import org.apache.commons.math3.util.Precision;
+
+/**
+ * Calculates the eigen decomposition of a real matrix.
+ *
+ * <p>The eigen decomposition of matrix A is a set of two matrices: V and D such that A = V &times;
+ * D &times; V<sup>T</sup>. A, V and D are all m &times; m matrices.
+ *
+ * <p>This class is similar in spirit to the <code>EigenvalueDecomposition</code> class from the <a
+ * href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library, with the following changes:
+ *
+ * <ul>
+ * <li>a {@link #getVT() getVt} method has been added,
+ * <li>two {@link #getRealEigenvalue(int) getRealEigenvalue} and {@link #getImagEigenvalue(int)
+ * getImagEigenvalue} methods to pick up a single eigenvalue have been added,
+ * <li>a {@link #getEigenvector(int) getEigenvector} method to pick up a single eigenvector has
+ * been added,
+ * <li>a {@link #getDeterminant() getDeterminant} method has been added.
+ * <li>a {@link #getSolver() getSolver} method has been added.
+ * </ul>
+ *
+ * <p>As of 3.1, this class supports general real matrices (both symmetric and non-symmetric):
+ *
+ * <p>If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is diagonal and the
+ * eigenvector matrix V is orthogonal, i.e. A = V.multiply(D.multiply(V.transpose())) and
+ * V.multiply(V.transpose()) equals the identity matrix.
+ *
+ * <p>If A is not symmetric, then the eigenvalue matrix D is block diagonal with the real
+ * eigenvalues in 1-by-1 blocks and any complex eigenvalues, lambda + i*mu, in 2-by-2 blocks:
+ *
+ * <pre>
+ * [lambda, mu ]
+ * [ -mu, lambda]
+ * </pre>
+ *
+ * The columns of V represent the eigenvectors in the sense that A*V = V*D, i.e. A.multiply(V)
+ * equals V.multiply(D). The matrix V may be badly conditioned, or even singular, so the validity of
+ * the equation A = V*D*inverse(V) depends upon the condition of V.
+ *
+ * <p>This implementation is based on the paper by A. Drubrulle, R.S. Martin and J.H. Wilkinson "The
+ * Implicit QL Algorithm" in Wilksinson and Reinsch (1971) Handbook for automatic computation, vol.
+ * 2, Linear algebra, Springer-Verlag, New-York
+ *
+ * @see <a href="http://mathworld.wolfram.com/EigenDecomposition.html">MathWorld</a>
+ * @see <a href="http://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix">Wikipedia</a>
+ * @since 2.0 (changed to concrete class in 3.0)
+ */
+public class EigenDecomposition {
+ /** Internally used epsilon criteria. */
+ private static final double EPSILON = 1e-12;
+
+ /** Maximum number of iterations accepted in the implicit QL transformation */
+ private byte maxIter = 30;
+
+ /** Main diagonal of the tridiagonal matrix. */
+ private double[] main;
+
+ /** Secondary diagonal of the tridiagonal matrix. */
+ private double[] secondary;
+
+ /** Transformer to tridiagonal (may be null if matrix is already tridiagonal). */
+ private TriDiagonalTransformer transformer;
+
+ /** Real part of the realEigenvalues. */
+ private double[] realEigenvalues;
+
+ /** Imaginary part of the realEigenvalues. */
+ private double[] imagEigenvalues;
+
+ /** Eigenvectors. */
+ private ArrayRealVector[] eigenvectors;
+
+ /** Cached value of V. */
+ private RealMatrix cachedV;
+
+ /** Cached value of D. */
+ private RealMatrix cachedD;
+
+ /** Cached value of Vt. */
+ private RealMatrix cachedVt;
+
+ /** Whether the matrix is symmetric. */
+ private final boolean isSymmetric;
+
+ /**
+ * Calculates the eigen decomposition of the given real matrix.
+ *
+ * <p>Supports decomposition of a general matrix since 3.1.
+ *
+ * @param matrix Matrix to decompose.
+ * @throws MaxCountExceededException if the algorithm fails to converge.
+ * @throws MathArithmeticException if the decomposition of a general matrix results in a matrix
+ * with zero norm
+ * @since 3.1
+ */
+ public EigenDecomposition(final RealMatrix matrix) throws MathArithmeticException {
+ final double symTol =
+ 10 * matrix.getRowDimension() * matrix.getColumnDimension() * Precision.EPSILON;
+ isSymmetric = MatrixUtils.isSymmetric(matrix, symTol);
+ if (isSymmetric) {
+ transformToTridiagonal(matrix);
+ findEigenVectors(transformer.getQ().getData());
+ } else {
+ final SchurTransformer t = transformToSchur(matrix);
+ findEigenVectorsFromSchur(t);
+ }
+ }
+
+ /**
+ * Calculates the eigen decomposition of the given real matrix.
+ *
+ * @param matrix Matrix to decompose.
+ * @param splitTolerance Dummy parameter (present for backward compatibility only).
+ * @throws MathArithmeticException if the decomposition of a general matrix results in a matrix
+ * with zero norm
+ * @throws MaxCountExceededException if the algorithm fails to converge.
+ * @deprecated in 3.1 (to be removed in 4.0) due to unused parameter
+ */
+ @Deprecated
+ public EigenDecomposition(final RealMatrix matrix, final double splitTolerance)
+ throws MathArithmeticException {
+ this(matrix);
+ }
+
+ /**
+ * Calculates the eigen decomposition of the symmetric tridiagonal matrix. The Householder
+ * matrix is assumed to be the identity matrix.
+ *
+ * @param main Main diagonal of the symmetric tridiagonal form.
+ * @param secondary Secondary of the tridiagonal form.
+ * @throws MaxCountExceededException if the algorithm fails to converge.
+ * @since 3.1
+ */
+ public EigenDecomposition(final double[] main, final double[] secondary) {
+ isSymmetric = true;
+ this.main = main.clone();
+ this.secondary = secondary.clone();
+ transformer = null;
+ final int size = main.length;
+ final double[][] z = new double[size][size];
+ for (int i = 0; i < size; i++) {
+ z[i][i] = 1.0;
+ }
+ findEigenVectors(z);
+ }
+
+ /**
+ * Calculates the eigen decomposition of the symmetric tridiagonal matrix. The Householder
+ * matrix is assumed to be the identity matrix.
+ *
+ * @param main Main diagonal of the symmetric tridiagonal form.
+ * @param secondary Secondary of the tridiagonal form.
+ * @param splitTolerance Dummy parameter (present for backward compatibility only).
+ * @throws MaxCountExceededException if the algorithm fails to converge.
+ * @deprecated in 3.1 (to be removed in 4.0) due to unused parameter
+ */
+ @Deprecated
+ public EigenDecomposition(
+ final double[] main, final double[] secondary, final double splitTolerance) {
+ this(main, secondary);
+ }
+
+ /**
+ * Gets the matrix V of the decomposition. V is an orthogonal matrix, i.e. its transpose is also
+ * its inverse. The columns of V are the eigenvectors of the original matrix. No assumption is
+ * made about the orientation of the system axes formed by the columns of V (e.g. in a
+ * 3-dimension space, V can form a left- or right-handed system).
+ *
+ * @return the V matrix.
+ */
+ public RealMatrix getV() {
+
+ if (cachedV == null) {
+ final int m = eigenvectors.length;
+ cachedV = MatrixUtils.createRealMatrix(m, m);
+ for (int k = 0; k < m; ++k) {
+ cachedV.setColumnVector(k, eigenvectors[k]);
+ }
+ }
+ // return the cached matrix
+ return cachedV;
+ }
+
+ /**
+ * Gets the block diagonal matrix D of the decomposition. D is a block diagonal matrix. Real
+ * eigenvalues are on the diagonal while complex values are on 2x2 blocks { {real +imaginary},
+ * {-imaginary, real} }.
+ *
+ * @return the D matrix.
+ * @see #getRealEigenvalues()
+ * @see #getImagEigenvalues()
+ */
+ public RealMatrix getD() {
+
+ if (cachedD == null) {
+ // cache the matrix for subsequent calls
+ cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues);
+
+ for (int i = 0; i < imagEigenvalues.length; i++) {
+ if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) > 0) {
+ cachedD.setEntry(i, i + 1, imagEigenvalues[i]);
+ } else if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) {
+ cachedD.setEntry(i, i - 1, imagEigenvalues[i]);
+ }
+ }
+ }
+ return cachedD;
+ }
+
+ /**
+ * Gets the transpose of the matrix V of the decomposition. V is an orthogonal matrix, i.e. its
+ * transpose is also its inverse. The columns of V are the eigenvectors of the original matrix.
+ * No assumption is made about the orientation of the system axes formed by the columns of V
+ * (e.g. in a 3-dimension space, V can form a left- or right-handed system).
+ *
+ * @return the transpose of the V matrix.
+ */
+ public RealMatrix getVT() {
+
+ if (cachedVt == null) {
+ final int m = eigenvectors.length;
+ cachedVt = MatrixUtils.createRealMatrix(m, m);
+ for (int k = 0; k < m; ++k) {
+ cachedVt.setRowVector(k, eigenvectors[k]);
+ }
+ }
+
+ // return the cached matrix
+ return cachedVt;
+ }
+
+ /**
+ * Returns whether the calculated eigen values are complex or real.
+ *
+ * <p>The method performs a zero check for each element of the {@link #getImagEigenvalues()}
+ * array and returns {@code true} if any element is not equal to zero.
+ *
+ * @return {@code true} if the eigen values are complex, {@code false} otherwise
+ * @since 3.1
+ */
+ public boolean hasComplexEigenvalues() {
+ for (int i = 0; i < imagEigenvalues.length; i++) {
+ if (!Precision.equals(imagEigenvalues[i], 0.0, EPSILON)) {
+ return true;
+ }
+ }
+ return false;
+ }
+
+ /**
+ * Gets a copy of the real parts of the eigenvalues of the original matrix.
+ *
+ * @return a copy of the real parts of the eigenvalues of the original matrix.
+ * @see #getD()
+ * @see #getRealEigenvalue(int)
+ * @see #getImagEigenvalues()
+ */
+ public double[] getRealEigenvalues() {
+ return realEigenvalues.clone();
+ }
+
+ /**
+ * Returns the real part of the i<sup>th</sup> eigenvalue of the original matrix.
+ *
+ * @param i index of the eigenvalue (counting from 0)
+ * @return real part of the i<sup>th</sup> eigenvalue of the original matrix.
+ * @see #getD()
+ * @see #getRealEigenvalues()
+ * @see #getImagEigenvalue(int)
+ */
+ public double getRealEigenvalue(final int i) {
+ return realEigenvalues[i];
+ }
+
+ /**
+ * Gets a copy of the imaginary parts of the eigenvalues of the original matrix.
+ *
+ * @return a copy of the imaginary parts of the eigenvalues of the original matrix.
+ * @see #getD()
+ * @see #getImagEigenvalue(int)
+ * @see #getRealEigenvalues()
+ */
+ public double[] getImagEigenvalues() {
+ return imagEigenvalues.clone();
+ }
+
+ /**
+ * Gets the imaginary part of the i<sup>th</sup> eigenvalue of the original matrix.
+ *
+ * @param i Index of the eigenvalue (counting from 0).
+ * @return the imaginary part of the i<sup>th</sup> eigenvalue of the original matrix.
+ * @see #getD()
+ * @see #getImagEigenvalues()
+ * @see #getRealEigenvalue(int)
+ */
+ public double getImagEigenvalue(final int i) {
+ return imagEigenvalues[i];
+ }
+
+ /**
+ * Gets a copy of the i<sup>th</sup> eigenvector of the original matrix.
+ *
+ * @param i Index of the eigenvector (counting from 0).
+ * @return a copy of the i<sup>th</sup> eigenvector of the original matrix.
+ * @see #getD()
+ */
+ public RealVector getEigenvector(final int i) {
+ return eigenvectors[i].copy();
+ }
+
+ /**
+ * Computes the determinant of the matrix.
+ *
+ * @return the determinant of the matrix.
+ */
+ public double getDeterminant() {
+ double determinant = 1;
+ for (double lambda : realEigenvalues) {
+ determinant *= lambda;
+ }
+ return determinant;
+ }
+
+ /**
+ * Computes the square-root of the matrix. This implementation assumes that the matrix is
+ * symmetric and positive definite.
+ *
+ * @return the square-root of the matrix.
+ * @throws MathUnsupportedOperationException if the matrix is not symmetric or not positive
+ * definite.
+ * @since 3.1
+ */
+ public RealMatrix getSquareRoot() {
+ if (!isSymmetric) {
+ throw new MathUnsupportedOperationException();
+ }
+
+ final double[] sqrtEigenValues = new double[realEigenvalues.length];
+ for (int i = 0; i < realEigenvalues.length; i++) {
+ final double eigen = realEigenvalues[i];
+ if (eigen <= 0) {
+ throw new MathUnsupportedOperationException();
+ }
+ sqrtEigenValues[i] = FastMath.sqrt(eigen);
+ }
+ final RealMatrix sqrtEigen = MatrixUtils.createRealDiagonalMatrix(sqrtEigenValues);
+ final RealMatrix v = getV();
+ final RealMatrix vT = getVT();
+
+ return v.multiply(sqrtEigen).multiply(vT);
+ }
+
+ /**
+ * Gets a solver for finding the A &times; X = B solution in exact linear sense.
+ *
+ * <p>Since 3.1, eigen decomposition of a general matrix is supported, but the {@link
+ * DecompositionSolver} only supports real eigenvalues.
+ *
+ * @return a solver
+ * @throws MathUnsupportedOperationException if the decomposition resulted in complex
+ * eigenvalues
+ */
+ public DecompositionSolver getSolver() {
+ if (hasComplexEigenvalues()) {
+ throw new MathUnsupportedOperationException();
+ }
+ return new Solver(realEigenvalues, imagEigenvalues, eigenvectors);
+ }
+
+ /** Specialized solver. */
+ private static class Solver implements DecompositionSolver {
+ /** Real part of the realEigenvalues. */
+ private double[] realEigenvalues;
+
+ /** Imaginary part of the realEigenvalues. */
+ private double[] imagEigenvalues;
+
+ /** Eigenvectors. */
+ private final ArrayRealVector[] eigenvectors;
+
+ /**
+ * Builds a solver from decomposed matrix.
+ *
+ * @param realEigenvalues Real parts of the eigenvalues.
+ * @param imagEigenvalues Imaginary parts of the eigenvalues.
+ * @param eigenvectors Eigenvectors.
+ */
+ private Solver(
+ final double[] realEigenvalues,
+ final double[] imagEigenvalues,
+ final ArrayRealVector[] eigenvectors) {
+ this.realEigenvalues = realEigenvalues;
+ this.imagEigenvalues = imagEigenvalues;
+ this.eigenvectors = eigenvectors;
+ }
+
+ /**
+ * Solves the linear equation A &times; X = B for symmetric matrices A.
+ *
+ * <p>This method only finds exact linear solutions, i.e. solutions for which ||A &times; X
+ * - B|| is exactly 0.
+ *
+ * @param b Right-hand side of the equation A &times; X = B.
+ * @return a Vector X that minimizes the two norm of A &times; X - B.
+ * @throws DimensionMismatchException if the matrices dimensions do not match.
+ * @throws SingularMatrixException if the decomposed matrix is singular.
+ */
+ public RealVector solve(final RealVector b) {
+ if (!isNonSingular()) {
+ throw new SingularMatrixException();
+ }
+
+ final int m = realEigenvalues.length;
+ if (b.getDimension() != m) {
+ throw new DimensionMismatchException(b.getDimension(), m);
+ }
+
+ final double[] bp = new double[m];
+ for (int i = 0; i < m; ++i) {
+ final ArrayRealVector v = eigenvectors[i];
+ final double[] vData = v.getDataRef();
+ final double s = v.dotProduct(b) / realEigenvalues[i];
+ for (int j = 0; j < m; ++j) {
+ bp[j] += s * vData[j];
+ }
+ }
+
+ return new ArrayRealVector(bp, false);
+ }
+
+ /** {@inheritDoc} */
+ public RealMatrix solve(RealMatrix b) {
+
+ if (!isNonSingular()) {
+ throw new SingularMatrixException();
+ }
+
+ final int m = realEigenvalues.length;
+ if (b.getRowDimension() != m) {
+ throw new DimensionMismatchException(b.getRowDimension(), m);
+ }
+
+ final int nColB = b.getColumnDimension();
+ final double[][] bp = new double[m][nColB];
+ final double[] tmpCol = new double[m];
+ for (int k = 0; k < nColB; ++k) {
+ for (int i = 0; i < m; ++i) {
+ tmpCol[i] = b.getEntry(i, k);
+ bp[i][k] = 0;
+ }
+ for (int i = 0; i < m; ++i) {
+ final ArrayRealVector v = eigenvectors[i];
+ final double[] vData = v.getDataRef();
+ double s = 0;
+ for (int j = 0; j < m; ++j) {
+ s += v.getEntry(j) * tmpCol[j];
+ }
+ s /= realEigenvalues[i];
+ for (int j = 0; j < m; ++j) {
+ bp[j][k] += s * vData[j];
+ }
+ }
+ }
+
+ return new Array2DRowRealMatrix(bp, false);
+ }
+
+ /**
+ * Checks whether the decomposed matrix is non-singular.
+ *
+ * @return true if the decomposed matrix is non-singular.
+ */
+ public boolean isNonSingular() {
+ double largestEigenvalueNorm = 0.0;
+ // Looping over all values (in case they are not sorted in decreasing
+ // order of their norm).
+ for (int i = 0; i < realEigenvalues.length; ++i) {
+ largestEigenvalueNorm = FastMath.max(largestEigenvalueNorm, eigenvalueNorm(i));
+ }
+ // Corner case: zero matrix, all exactly 0 eigenvalues
+ if (largestEigenvalueNorm == 0.0) {
+ return false;
+ }
+ for (int i = 0; i < realEigenvalues.length; ++i) {
+ // Looking for eigenvalues that are 0, where we consider anything much much smaller
+ // than the largest eigenvalue to be effectively 0.
+ if (Precision.equals(eigenvalueNorm(i) / largestEigenvalueNorm, 0, EPSILON)) {
+ return false;
+ }
+ }
+ return true;
+ }
+
+ /**
+ * @param i which eigenvalue to find the norm of
+ * @return the norm of ith (complex) eigenvalue.
+ */
+ private double eigenvalueNorm(int i) {
+ final double re = realEigenvalues[i];
+ final double im = imagEigenvalues[i];
+ return FastMath.sqrt(re * re + im * im);
+ }
+
+ /**
+ * Get the inverse of the decomposed matrix.
+ *
+ * @return the inverse matrix.
+ * @throws SingularMatrixException if the decomposed matrix is singular.
+ */
+ public RealMatrix getInverse() {
+ if (!isNonSingular()) {
+ throw new SingularMatrixException();
+ }
+
+ final int m = realEigenvalues.length;
+ final double[][] invData = new double[m][m];
+
+ for (int i = 0; i < m; ++i) {
+ final double[] invI = invData[i];
+ for (int j = 0; j < m; ++j) {
+ double invIJ = 0;
+ for (int k = 0; k < m; ++k) {
+ final double[] vK = eigenvectors[k].getDataRef();
+ invIJ += vK[i] * vK[j] / realEigenvalues[k];
+ }
+ invI[j] = invIJ;
+ }
+ }
+ return MatrixUtils.createRealMatrix(invData);
+ }
+ }
+
+ /**
+ * Transforms the matrix to tridiagonal form.
+ *
+ * @param matrix Matrix to transform.
+ */
+ private void transformToTridiagonal(final RealMatrix matrix) {
+ // transform the matrix to tridiagonal
+ transformer = new TriDiagonalTransformer(matrix);
+ main = transformer.getMainDiagonalRef();
+ secondary = transformer.getSecondaryDiagonalRef();
+ }
+
+ /**
+ * Find eigenvalues and eigenvectors (Dubrulle et al., 1971)
+ *
+ * @param householderMatrix Householder matrix of the transformation to tridiagonal form.
+ */
+ private void findEigenVectors(final double[][] householderMatrix) {
+ final double[][] z = householderMatrix.clone();
+ final int n = main.length;
+ realEigenvalues = new double[n];
+ imagEigenvalues = new double[n];
+ final double[] e = new double[n];
+ for (int i = 0; i < n - 1; i++) {
+ realEigenvalues[i] = main[i];
+ e[i] = secondary[i];
+ }
+ realEigenvalues[n - 1] = main[n - 1];
+ e[n - 1] = 0;
+
+ // Determine the largest main and secondary value in absolute term.
+ double maxAbsoluteValue = 0;
+ for (int i = 0; i < n; i++) {
+ if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) {
+ maxAbsoluteValue = FastMath.abs(realEigenvalues[i]);
+ }
+ if (FastMath.abs(e[i]) > maxAbsoluteValue) {
+ maxAbsoluteValue = FastMath.abs(e[i]);
+ }
+ }
+ // Make null any main and secondary value too small to be significant
+ if (maxAbsoluteValue != 0) {
+ for (int i = 0; i < n; i++) {
+ if (FastMath.abs(realEigenvalues[i]) <= Precision.EPSILON * maxAbsoluteValue) {
+ realEigenvalues[i] = 0;
+ }
+ if (FastMath.abs(e[i]) <= Precision.EPSILON * maxAbsoluteValue) {
+ e[i] = 0;
+ }
+ }
+ }
+
+ for (int j = 0; j < n; j++) {
+ int its = 0;
+ int m;
+ do {
+ for (m = j; m < n - 1; m++) {
+ double delta =
+ FastMath.abs(realEigenvalues[m]) + FastMath.abs(realEigenvalues[m + 1]);
+ if (FastMath.abs(e[m]) + delta == delta) {
+ break;
+ }
+ }
+ if (m != j) {
+ if (its == maxIter) {
+ throw new MaxCountExceededException(
+ LocalizedFormats.CONVERGENCE_FAILED, maxIter);
+ }
+ its++;
+ double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]);
+ double t = FastMath.sqrt(1 + q * q);
+ if (q < 0.0) {
+ q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q - t);
+ } else {
+ q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q + t);
+ }
+ double u = 0.0;
+ double s = 1.0;
+ double c = 1.0;
+ int i;
+ for (i = m - 1; i >= j; i--) {
+ double p = s * e[i];
+ double h = c * e[i];
+ if (FastMath.abs(p) >= FastMath.abs(q)) {
+ c = q / p;
+ t = FastMath.sqrt(c * c + 1.0);
+ e[i + 1] = p * t;
+ s = 1.0 / t;
+ c *= s;
+ } else {
+ s = p / q;
+ t = FastMath.sqrt(s * s + 1.0);
+ e[i + 1] = q * t;
+ c = 1.0 / t;
+ s *= c;
+ }
+ if (e[i + 1] == 0.0) {
+ realEigenvalues[i + 1] -= u;
+ e[m] = 0.0;
+ break;
+ }
+ q = realEigenvalues[i + 1] - u;
+ t = (realEigenvalues[i] - q) * s + 2.0 * c * h;
+ u = s * t;
+ realEigenvalues[i + 1] = q + u;
+ q = c * t - h;
+ for (int ia = 0; ia < n; ia++) {
+ p = z[ia][i + 1];
+ z[ia][i + 1] = s * z[ia][i] + c * p;
+ z[ia][i] = c * z[ia][i] - s * p;
+ }
+ }
+ if (t == 0.0 && i >= j) {
+ continue;
+ }
+ realEigenvalues[j] -= u;
+ e[j] = q;
+ e[m] = 0.0;
+ }
+ } while (m != j);
+ }
+
+ // Sort the eigen values (and vectors) in increase order
+ for (int i = 0; i < n; i++) {
+ int k = i;
+ double p = realEigenvalues[i];
+ for (int j = i + 1; j < n; j++) {
+ if (realEigenvalues[j] > p) {
+ k = j;
+ p = realEigenvalues[j];
+ }
+ }
+ if (k != i) {
+ realEigenvalues[k] = realEigenvalues[i];
+ realEigenvalues[i] = p;
+ for (int j = 0; j < n; j++) {
+ p = z[j][i];
+ z[j][i] = z[j][k];
+ z[j][k] = p;
+ }
+ }
+ }
+
+ // Determine the largest eigen value in absolute term.
+ maxAbsoluteValue = 0;
+ for (int i = 0; i < n; i++) {
+ if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) {
+ maxAbsoluteValue = FastMath.abs(realEigenvalues[i]);
+ }
+ }
+ // Make null any eigen value too small to be significant
+ if (maxAbsoluteValue != 0.0) {
+ for (int i = 0; i < n; i++) {
+ if (FastMath.abs(realEigenvalues[i]) < Precision.EPSILON * maxAbsoluteValue) {
+ realEigenvalues[i] = 0;
+ }
+ }
+ }
+ eigenvectors = new ArrayRealVector[n];
+ final double[] tmp = new double[n];
+ for (int i = 0; i < n; i++) {
+ for (int j = 0; j < n; j++) {
+ tmp[j] = z[j][i];
+ }
+ eigenvectors[i] = new ArrayRealVector(tmp);
+ }
+ }
+
+ /**
+ * Transforms the matrix to Schur form and calculates the eigenvalues.
+ *
+ * @param matrix Matrix to transform.
+ * @return the {@link SchurTransformer Shur transform} for this matrix
+ */
+ private SchurTransformer transformToSchur(final RealMatrix matrix) {
+ final SchurTransformer schurTransform = new SchurTransformer(matrix);
+ final double[][] matT = schurTransform.getT().getData();
+
+ realEigenvalues = new double[matT.length];
+ imagEigenvalues = new double[matT.length];
+
+ for (int i = 0; i < realEigenvalues.length; i++) {
+ if (i == (realEigenvalues.length - 1)
+ || Precision.equals(matT[i + 1][i], 0.0, EPSILON)) {
+ realEigenvalues[i] = matT[i][i];
+ } else {
+ final double x = matT[i + 1][i + 1];
+ final double p = 0.5 * (matT[i][i] - x);
+ final double z =
+ FastMath.sqrt(FastMath.abs(p * p + matT[i + 1][i] * matT[i][i + 1]));
+ realEigenvalues[i] = x + p;
+ imagEigenvalues[i] = z;
+ realEigenvalues[i + 1] = x + p;
+ imagEigenvalues[i + 1] = -z;
+ i++;
+ }
+ }
+ return schurTransform;
+ }
+
+ /**
+ * Performs a division of two complex numbers.
+ *
+ * @param xr real part of the first number
+ * @param xi imaginary part of the first number
+ * @param yr real part of the second number
+ * @param yi imaginary part of the second number
+ * @return result of the complex division
+ */
+ private Complex cdiv(final double xr, final double xi, final double yr, final double yi) {
+ return new Complex(xr, xi).divide(new Complex(yr, yi));
+ }
+
+ /**
+ * Find eigenvectors from a matrix transformed to Schur form.
+ *
+ * @param schur the schur transformation of the matrix
+ * @throws MathArithmeticException if the Schur form has a norm of zero
+ */
+ private void findEigenVectorsFromSchur(final SchurTransformer schur)
+ throws MathArithmeticException {
+ final double[][] matrixT = schur.getT().getData();
+ final double[][] matrixP = schur.getP().getData();
+
+ final int n = matrixT.length;
+
+ // compute matrix norm
+ double norm = 0.0;
+ for (int i = 0; i < n; i++) {
+ for (int j = FastMath.max(i - 1, 0); j < n; j++) {
+ norm += FastMath.abs(matrixT[i][j]);
+ }
+ }
+
+ // we can not handle a matrix with zero norm
+ if (Precision.equals(norm, 0.0, EPSILON)) {
+ throw new MathArithmeticException(LocalizedFormats.ZERO_NORM);
+ }
+
+ // Backsubstitute to find vectors of upper triangular form
+
+ double r = 0.0;
+ double s = 0.0;
+ double z = 0.0;
+
+ for (int idx = n - 1; idx >= 0; idx--) {
+ double p = realEigenvalues[idx];
+ double q = imagEigenvalues[idx];
+
+ if (Precision.equals(q, 0.0)) {
+ // Real vector
+ int l = idx;
+ matrixT[idx][idx] = 1.0;
+ for (int i = idx - 1; i >= 0; i--) {
+ double w = matrixT[i][i] - p;
+ r = 0.0;
+ for (int j = l; j <= idx; j++) {
+ r += matrixT[i][j] * matrixT[j][idx];
+ }
+ if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) {
+ z = w;
+ s = r;
+ } else {
+ l = i;
+ if (Precision.equals(imagEigenvalues[i], 0.0)) {
+ if (w != 0.0) {
+ matrixT[i][idx] = -r / w;
+ } else {
+ matrixT[i][idx] = -r / (Precision.EPSILON * norm);
+ }
+ } else {
+ // Solve real equations
+ double x = matrixT[i][i + 1];
+ double y = matrixT[i + 1][i];
+ q =
+ (realEigenvalues[i] - p) * (realEigenvalues[i] - p)
+ + imagEigenvalues[i] * imagEigenvalues[i];
+ double t = (x * s - z * r) / q;
+ matrixT[i][idx] = t;
+ if (FastMath.abs(x) > FastMath.abs(z)) {
+ matrixT[i + 1][idx] = (-r - w * t) / x;
+ } else {
+ matrixT[i + 1][idx] = (-s - y * t) / z;
+ }
+ }
+
+ // Overflow control
+ double t = FastMath.abs(matrixT[i][idx]);
+ if ((Precision.EPSILON * t) * t > 1) {
+ for (int j = i; j <= idx; j++) {
+ matrixT[j][idx] /= t;
+ }
+ }
+ }
+ }
+ } else if (q < 0.0) {
+ // Complex vector
+ int l = idx - 1;
+
+ // Last vector component imaginary so matrix is triangular
+ if (FastMath.abs(matrixT[idx][idx - 1]) > FastMath.abs(matrixT[idx - 1][idx])) {
+ matrixT[idx - 1][idx - 1] = q / matrixT[idx][idx - 1];
+ matrixT[idx - 1][idx] = -(matrixT[idx][idx] - p) / matrixT[idx][idx - 1];
+ } else {
+ final Complex result =
+ cdiv(0.0, -matrixT[idx - 1][idx], matrixT[idx - 1][idx - 1] - p, q);
+ matrixT[idx - 1][idx - 1] = result.getReal();
+ matrixT[idx - 1][idx] = result.getImaginary();
+ }
+
+ matrixT[idx][idx - 1] = 0.0;
+ matrixT[idx][idx] = 1.0;
+
+ for (int i = idx - 2; i >= 0; i--) {
+ double ra = 0.0;
+ double sa = 0.0;
+ for (int j = l; j <= idx; j++) {
+ ra += matrixT[i][j] * matrixT[j][idx - 1];
+ sa += matrixT[i][j] * matrixT[j][idx];
+ }
+ double w = matrixT[i][i] - p;
+
+ if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) {
+ z = w;
+ r = ra;
+ s = sa;
+ } else {
+ l = i;
+ if (Precision.equals(imagEigenvalues[i], 0.0)) {
+ final Complex c = cdiv(-ra, -sa, w, q);
+ matrixT[i][idx - 1] = c.getReal();
+ matrixT[i][idx] = c.getImaginary();
+ } else {
+ // Solve complex equations
+ double x = matrixT[i][i + 1];
+ double y = matrixT[i + 1][i];
+ double vr =
+ (realEigenvalues[i] - p) * (realEigenvalues[i] - p)
+ + imagEigenvalues[i] * imagEigenvalues[i]
+ - q * q;
+ final double vi = (realEigenvalues[i] - p) * 2.0 * q;
+ if (Precision.equals(vr, 0.0) && Precision.equals(vi, 0.0)) {
+ vr =
+ Precision.EPSILON
+ * norm
+ * (FastMath.abs(w)
+ + FastMath.abs(q)
+ + FastMath.abs(x)
+ + FastMath.abs(y)
+ + FastMath.abs(z));
+ }
+ final Complex c =
+ cdiv(x * r - z * ra + q * sa, x * s - z * sa - q * ra, vr, vi);
+ matrixT[i][idx - 1] = c.getReal();
+ matrixT[i][idx] = c.getImaginary();
+
+ if (FastMath.abs(x) > (FastMath.abs(z) + FastMath.abs(q))) {
+ matrixT[i + 1][idx - 1] =
+ (-ra - w * matrixT[i][idx - 1] + q * matrixT[i][idx]) / x;
+ matrixT[i + 1][idx] =
+ (-sa - w * matrixT[i][idx] - q * matrixT[i][idx - 1]) / x;
+ } else {
+ final Complex c2 =
+ cdiv(
+ -r - y * matrixT[i][idx - 1],
+ -s - y * matrixT[i][idx],
+ z,
+ q);
+ matrixT[i + 1][idx - 1] = c2.getReal();
+ matrixT[i + 1][idx] = c2.getImaginary();
+ }
+ }
+
+ // Overflow control
+ double t =
+ FastMath.max(
+ FastMath.abs(matrixT[i][idx - 1]),
+ FastMath.abs(matrixT[i][idx]));
+ if ((Precision.EPSILON * t) * t > 1) {
+ for (int j = i; j <= idx; j++) {
+ matrixT[j][idx - 1] /= t;
+ matrixT[j][idx] /= t;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // Back transformation to get eigenvectors of original matrix
+ for (int j = n - 1; j >= 0; j--) {
+ for (int i = 0; i <= n - 1; i++) {
+ z = 0.0;
+ for (int k = 0; k <= FastMath.min(j, n - 1); k++) {
+ z += matrixP[i][k] * matrixT[k][j];
+ }
+ matrixP[i][j] = z;
+ }
+ }
+
+ eigenvectors = new ArrayRealVector[n];
+ final double[] tmp = new double[n];
+ for (int i = 0; i < n; i++) {
+ for (int j = 0; j < n; j++) {
+ tmp[j] = matrixP[j][i];
+ }
+ eigenvectors[i] = new ArrayRealVector(tmp);
+ }
+ }
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-14 21:55:57 +0000