Gilles Sadowski
parent
aea043cffd
commit
9b7a2c8edc
|
|
@ -1,945 +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.linear;
|
||||
|
||||
import org.apache.commons.numbers.complex.Complex;
|
||||
import org.apache.commons.numbers.core.Precision;
|
||||
import org.apache.commons.math4.exception.DimensionMismatchException;
|
||||
import org.apache.commons.math4.exception.MathArithmeticException;
|
||||
import org.apache.commons.math4.exception.MathUnsupportedOperationException;
|
||||
import org.apache.commons.math4.exception.MaxCountExceededException;
|
||||
import org.apache.commons.math4.exception.util.LocalizedFormats;
|
||||
import org.apache.commons.math4.util.FastMath;
|
||||
|
||||
/**
|
||||
* 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 × D × V<sup>T</sup>.
|
||||
* A, V and D are all m × m matrices.
|
||||
* <p>
|
||||
* This class is similar in spirit to the {@code EigenvalueDecomposition}
|
||||
* 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>
|
||||
* <li>two {@link #getRealEigenvalue(int) getRealEigenvalue} and
|
||||
* {@link #getImagEigenvalue(int) getImagEigenvalue} methods to pick up a
|
||||
* single eigenvalue have been added,</li>
|
||||
* <li>a {@link #getEigenvector(int) getEigenvector} method to pick up a
|
||||
* single eigenvector has been added,</li>
|
||||
* <li>a {@link #getDeterminant() getDeterminant} method has been added.</li>
|
||||
* <li>a {@link #getSolver() getSolver} method has been added.</li>
|
||||
* </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.
|
||||
* {@code A = V.multiply(D.multiply(V.transpose()))} and
|
||||
* {@code V.multiply(V.transpose())} equals the identity matrix.
|
||||
* </p>
|
||||
* <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 {@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.
|
||||
* <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 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.
|
||||
* <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 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.
|
||||
* <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 × 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 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.
|
||||
* <p>
|
||||
* This method only finds exact linear solutions, i.e. solutions for
|
||||
* which ||A × X - B|| is exactly 0.
|
||||
* </p>
|
||||
*
|
||||
* @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);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -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.
|
||||
*
|
||||
* <p>
|
||||
* 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
|
||||
* (<a href="http://cm.bell-labs.com/cm/cs/doc/96/4-02.ps.gz">Direct
|
||||
* Search Methods: Once Scorned, Now Respectable</a>), 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
|
||||
* <em>not too bad</em> 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.
|
||||
* </p>
|
||||
* <p>
|
||||
* 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.
|
||||
* </p>
|
||||
* <p>
|
||||
* The simplex update procedure ({@link NelderMeadSimplex} or
|
||||
* {@link MultiDirectionalSimplex}) must be passed to the
|
||||
* {@code optimize} method.
|
||||
* </p>
|
||||
* <p>
|
||||
* 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.
|
||||
* </p>
|
||||
* <p>
|
||||
* Convergence is checked by providing the <em>worst</em> points of
|
||||
* previous and current simplex to the convergence checker, not the best
|
||||
* ones.
|
||||
* </p>
|
||||
* <p>
|
||||
* 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}.
|
||||
* <br>
|
||||
* The call to {@link #optimize(OptimizationData[]) optimize} will throw
|
||||
* {@link MathUnsupportedOperationException} if bounds are passed to it.
|
||||
* </p>
|
||||
*
|
||||
* @since 3.0
|
||||
*/
|
||||
public class SimplexOptimizer extends MultivariateOptimizer {
|
||||
/** Simplex update rule. */
|
||||
private AbstractSimplex simplex;
|
||||
/** Simulated annealing setup. */
|
||||
private SimulatedAnnealing annealing;
|
||||
/** Overall best. */
|
||||
private PointValuePair best;
|
||||
|
||||
/**
|
||||
* @param checker Convergence checker.
|
||||
*/
|
||||
public SimplexOptimizer(ConvergenceChecker<PointValuePair> checker) {
|
||||
super(checker);
|
||||
}
|
||||
|
||||
/**
|
||||
* @param rel Relative threshold.
|
||||
* @param abs Absolute threshold.
|
||||
*/
|
||||
public SimplexOptimizer(double rel, double abs) {
|
||||
this(new SimpleValueChecker(rel, abs));
|
||||
}
|
||||
|
||||
/**
|
||||
* {@inheritDoc}
|
||||
*
|
||||
* @param optData Optimization data. In addition to those documented in
|
||||
* {@link MultivariateOptimizer#parseOptimizationData(OptimizationData[])
|
||||
* MultivariateOptimizer}, this method will register the following data:
|
||||
* <ul>
|
||||
* <li>{@link AbstractSimplex}</li>
|
||||
* <li>{@link SimulatedAnnealing}</li>
|
||||
* </ul>
|
||||
* @return {@inheritDoc}
|
||||
*/
|
||||
@Override
|
||||
public PointValuePair optimize(OptimizationData... optData) {
|
||||
// Set up base class and perform computation.
|
||||
return super.optimize(optData);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
@Override
|
||||
protected PointValuePair doOptimize() {
|
||||
checkParameters();
|
||||
|
||||
// Indirect call to "computeObjectiveValue" in order to update the
|
||||
// evaluations counter.
|
||||
final MultivariateFunction evalFunc
|
||||
= new MultivariateFunction() {
|
||||
/** {@inheritDoc} */
|
||||
@Override
|
||||
public double value(double[] point) {
|
||||
return computeObjectiveValue(point);
|
||||
}
|
||||
};
|
||||
|
||||
final boolean isMinim = getGoalType() == GoalType.MINIMIZE;
|
||||
final Comparator<PointValuePair> comparator
|
||||
= new Comparator<PointValuePair>() {
|
||||
/** {@inheritDoc} */
|
||||
@Override
|
||||
public int compare(final PointValuePair o1,
|
||||
final PointValuePair o2) {
|
||||
final double v1 = o1.getValue();
|
||||
final double v2 = o2.getValue();
|
||||
return isMinim ? Double.compare(v1, v2) : Double.compare(v2, v1);
|
||||
}
|
||||
};
|
||||
|
||||
// Initialize search.
|
||||
simplex.build(getStartPoint());
|
||||
simplex.evaluate(evalFunc, comparator);
|
||||
final UniformRandomProvider rng = annealing != null ?
|
||||
RandomSource.create(RandomSource.KISS) :
|
||||
null;
|
||||
|
||||
PointValuePair[] previous = null;
|
||||
int iteration = 0;
|
||||
final ConvergenceChecker<PointValuePair> checker = getConvergenceChecker();
|
||||
while (true) {
|
||||
iteration = getIterations();
|
||||
if (iteration > 0) {
|
||||
boolean converged = true;
|
||||
for (int i = 0; i < simplex.getSize(); i++) {
|
||||
PointValuePair prev = previous[i];
|
||||
converged = converged &&
|
||||
checker.converged(iteration, prev, simplex.getPoint(i));
|
||||
|
||||
if (!converged) {
|
||||
// Short circuit, since "converged" will stay "false".
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (converged) {
|
||||
System.out.println(" saAcceptCount=" + saAcceptCount); // XXX
|
||||
// We have found an optimum.
|
||||
return best;
|
||||
}
|
||||
}
|
||||
|
||||
// We still need to search.
|
||||
previous = simplex.getPoints();
|
||||
simplex.iterate(evalFunc, comparator);
|
||||
|
||||
// Track best point.
|
||||
final int bestIndex = 0; // Index of best point.
|
||||
if (best == null ||
|
||||
comparator.compare(best, simplex.getPoint(bestIndex)) > 0) {
|
||||
best = simplex.getPoint(bestIndex);
|
||||
}
|
||||
|
||||
if (annealing != null) {
|
||||
// Simulated annealing step.
|
||||
simulatedAnnealing(iteration,
|
||||
evalFunc,
|
||||
isMinim,
|
||||
rng);
|
||||
}
|
||||
|
||||
incrementIterationCount();
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Scans the list of (required and optional) optimization data that
|
||||
* characterize the problem.
|
||||
*
|
||||
* @param optData Optimization data.
|
||||
* The following data will be looked for:
|
||||
* <ul>
|
||||
* <li>{@link AbstractSimplex}</li>
|
||||
* <li>{@link SimulatedAnnealing}</li>
|
||||
* </ul>
|
||||
*/
|
||||
@Override
|
||||
protected void parseOptimizationData(OptimizationData... optData) {
|
||||
// Allow base class to register its own data.
|
||||
super.parseOptimizationData(optData);
|
||||
|
||||
// The existing values (as set by the previous call) are reused if
|
||||
// not provided in the argument list.
|
||||
for (OptimizationData data : optData) {
|
||||
if (data instanceof AbstractSimplex) {
|
||||
simplex = (AbstractSimplex) data;
|
||||
continue;
|
||||
}
|
||||
if (data instanceof SimulatedAnnealing) {
|
||||
annealing = (SimulatedAnnealing) data;
|
||||
continue;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* @throws MathUnsupportedOperationException if bounds were passed to the
|
||||
* {@link #optimize(OptimizationData[]) optimize} method.
|
||||
* @throws NullArgumentException if no initial simplex was passed to the
|
||||
* {@link #optimize(OptimizationData[]) optimize} method.
|
||||
*/
|
||||
private void checkParameters() {
|
||||
if (simplex == null) {
|
||||
throw new NullArgumentException();
|
||||
}
|
||||
if (getLowerBound() != null ||
|
||||
getUpperBound() != null) {
|
||||
throw new MathUnsupportedOperationException(LocalizedFormats.CONSTRAINT);
|
||||
}
|
||||
}
|
||||
|
||||
private int saAcceptCount = 0; // XXX
|
||||
/**
|
||||
* Perform the annealing step (possibly replacing simplex's best point).
|
||||
*
|
||||
* @param iteration Current iteration.
|
||||
* @param evalFunc Evaluation function.
|
||||
* @param isMinim Whether a minimization is performed.
|
||||
* @param rng RNG.
|
||||
*/
|
||||
private void simulatedAnnealing(int iteration,
|
||||
MultivariateFunction evalFunc,
|
||||
boolean isMinim,
|
||||
UniformRandomProvider rng) {
|
||||
if (iteration > annealing.getIterations()) {
|
||||
return; // Do nothing.
|
||||
}
|
||||
|
||||
// Construct alternative state.
|
||||
final int bestIndex = 0; // Index of best point.
|
||||
final PointValuePair alt = alternativeState(simplex,
|
||||
bestIndex,
|
||||
evalFunc,
|
||||
rng);
|
||||
|
||||
if (annealing.accept(simplex.getPoint(bestIndex).getValue(),
|
||||
alt.getValue(),
|
||||
isMinim,
|
||||
iteration)) {
|
||||
System.out.println("eO=" + simplex.getPoint(bestIndex).getValue()); // XXX
|
||||
System.out.println("eN=" + alt.getValue()); // XXX
|
||||
++saAcceptCount; // XXX
|
||||
|
||||
// Modify best point of the current simplex.
|
||||
simplex.setPoint(bestIndex, alt);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Creates a state that could replace the one stored at
|
||||
* {@code replaceIndex} in the given {@code simplex}.
|
||||
*
|
||||
* @param simplex Current simplex.
|
||||
* @param replaceIndex Index of the simplex point that will potentially be
|
||||
* replace by the new state.
|
||||
* @param evalFunc Evaluation function.
|
||||
* @param rng RNG.
|
||||
* @return a new state.
|
||||
*/
|
||||
private static PointValuePair alternativeState(AbstractSimplex simplex,
|
||||
int replaceIndex,
|
||||
MultivariateFunction evalFunc,
|
||||
UniformRandomProvider rng) {
|
||||
final PointValuePair[] points = simplex.getPoints();
|
||||
final int numPoints = points.length;
|
||||
final int spaceDim = numPoints - 1;
|
||||
|
||||
// Compute mean coordinate offsets from the point to replace
|
||||
// to all the other points.
|
||||
final double[] coord = new double[spaceDim];
|
||||
final double[] replaceCoord = points[replaceIndex].getPointRef();
|
||||
for (int j = 0; j < numPoints; j++) {
|
||||
if (j == replaceIndex) {
|
||||
continue;
|
||||
}
|
||||
final double[] c = points[j].getPointRef();
|
||||
for (int i = 0; i < spaceDim; i++) {
|
||||
coord[i] += c[i] - replaceCoord[i];
|
||||
}
|
||||
}
|
||||
|
||||
for (int i = 0; i < spaceDim; i++) {
|
||||
coord[i] /= numPoints; // Mean coordinate offset.
|
||||
coord[i] = replaceCoord[i] + (rng.nextDouble() - 0.5) * coord[i];
|
||||
}
|
||||
|
||||
return new PointValuePair(coord, evalFunc.value(coord), false);
|
||||
}
|
||||
}
|
||||
Loading…
Reference in New Issue