- * The eigen decomposition of matrix A is a set of two matrices: - * V and D such that A = V × D × VT. - * A, V and D are all m × m matrices. - *
- * This class is similar in spirit to the {@code EigenvalueDecomposition} - * class from the JAMA - * library, with the following changes: - *
- * As of 3.1, this class supports general real matrices (both symmetric and non-symmetric): - *
- * 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. - * {@code A = V.multiply(D.multiply(V.transpose()))} and - * {@code V.multiply(V.transpose())} equals the identity matrix. - *
- *- * 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: - *
- * [lambda, mu ] - * [ -mu, lambda] - *- * The columns of V represent the eigenvectors in the sense that {@code 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 {@code A = V*D*inverse(V)} depends upon the condition of V. - *
- * 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 MathWorld - * @see Wikipedia - * @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 static final byte MAX_ITER = 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. - *
- * 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 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); - } - - /** - * 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.createRealMatrixWithDiagonal(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. - *
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 ith eigenvalue of the original - * matrix. - * - * @param i index of the eigenvalue (counting from 0) - * @return real part of the ith 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 ith eigenvalue of the original - * matrix. - * - * @param i Index of the eigenvalue (counting from 0). - * @return the imaginary part of the ith 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 ith eigenvector of the original matrix. - * - * @param i Index of the eigenvector (counting from 0). - * @return a copy of the ith 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 × X = B solution in exact - * linear sense. - *
- * 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 final double[] realEigenvalues; - /** Imaginary part of the realEigenvalues. */ - private final 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 × X = B for symmetric matrices A. - *
- * This method only finds exact linear solutions, i.e. solutions for - * which ||A × X - B|| is exactly 0. - *
- * - * @param b Right-hand side of the equation A × X = B. - * @return a Vector X that minimizes the two norm of A × X - B. - * - * @throws DimensionMismatchException if the matrices dimensions do not match. - * @throws SingularMatrixException if the decomposed matrix is singular. - */ - @Override - 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} */ - @Override - 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. - */ - @Override - 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. - */ - @Override - 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 == MAX_ITER) { - throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED, - MAX_ITER); - } - 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 Complex.ofCartesian(xr, xi).divide(Complex.ofCartesian(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++) { - System.out.println("Eigenvector " + i + ": "); // XXX - for (int j = 0; j < n; j++) { - tmp[j] = matrixP[j][i]; - System.out.print(tmp[j] + "\t"); // XXX - } - System.out.println(); // XXX - eigenvectors[i] = new ArrayRealVector(tmp); - } - } -} diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizer.java.DEBUG b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizer.java.DEBUG deleted file mode 100644 index c6b442093..000000000 --- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizer.java.DEBUG +++ /dev/null @@ -1,339 +0,0 @@ -/* - * 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.math4.optim.nonlinear.scalar.noderiv; - -import java.util.Comparator; - -import org.apache.commons.rng.UniformRandomProvider; -import org.apache.commons.rng.simple.RandomSource; -import org.apache.commons.math4.analysis.MultivariateFunction; -import org.apache.commons.math4.exception.MathUnsupportedOperationException; -import org.apache.commons.math4.exception.NullArgumentException; -import org.apache.commons.math4.exception.util.LocalizedFormats; -import org.apache.commons.math4.optim.ConvergenceChecker; -import org.apache.commons.math4.optim.OptimizationData; -import org.apache.commons.math4.optim.PointValuePair; -import org.apache.commons.math4.optim.SimpleValueChecker; -import org.apache.commons.math4.optim.nonlinear.scalar.GoalType; -import org.apache.commons.math4.optim.nonlinear.scalar.SimulatedAnnealing; -import org.apache.commons.math4.optim.nonlinear.scalar.MultivariateOptimizer; - -/** - * This class implements simplex-based direct search optimization. - * - *- * Direct search methods only use objective function values, they do - * not need derivatives and don't either try to compute approximation - * of the derivatives. According to a 1996 paper by Margaret H. Wright - * (Direct - * Search Methods: Once Scorned, Now Respectable), they are used - * when either the computation of the derivative is impossible (noisy - * functions, unpredictable discontinuities) or difficult (complexity, - * computation cost). In the first cases, rather than an optimum, a - * not too bad point is desired. In the latter cases, an - * optimum is desired but cannot be reasonably found. In all cases - * direct search methods can be useful. - *
- *- * Simplex-based direct search methods are based on comparison of - * the objective function values at the vertices of a simplex (which is a - * set of n+1 points in dimension n) that is updated by the algorithms - * steps. - *
- *- * The simplex update procedure ({@link NelderMeadSimplex} or - * {@link MultiDirectionalSimplex}) must be passed to the - * {@code optimize} method. - *
- *- * Each call to {@code optimize} will re-use the start configuration of - * the current simplex and move it such that its first vertex is at the - * provided start point of the optimization. - * If the {@code optimize} method is called to solve a different problem - * and the number of parameters change, the simplex must be re-initialized - * to one with the appropriate dimensions. - *
- *- * Convergence is checked by providing the worst points of - * previous and current simplex to the convergence checker, not the best - * ones. - *
- *
- * This simplex optimizer implementation does not directly support constrained
- * optimization with simple bounds; so, for such optimizations, either a more
- * dedicated algorithm must be used like
- * {@link CMAESOptimizer} or {@link BOBYQAOptimizer}, or the objective
- * function must be wrapped in an adapter like
- * {@link org.apache.commons.math4.optim.nonlinear.scalar.MultivariateFunctionMappingAdapter
- * MultivariateFunctionMappingAdapter} or
- * {@link org.apache.commons.math4.optim.nonlinear.scalar.MultivariateFunctionPenaltyAdapter
- * MultivariateFunctionPenaltyAdapter}.
- *
- * The call to {@link #optimize(OptimizationData[]) optimize} will throw
- * {@link MathUnsupportedOperationException} if bounds are passed to it.
- *