/*
* 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.math.optimization.general;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MaxEvaluationsExceededException;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.analysis.DifferentiableMultivariateVectorialFunction;
import org.apache.commons.math.analysis.MultivariateMatrixFunction;
import org.apache.commons.math.linear.InvalidMatrixException;
import org.apache.commons.math.linear.LUDecompositionImpl;
import org.apache.commons.math.linear.MatrixUtils;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.optimization.OptimizationException;
import org.apache.commons.math.optimization.SimpleVectorialValueChecker;
import org.apache.commons.math.optimization.VectorialConvergenceChecker;
import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer;
import org.apache.commons.math.optimization.VectorialPointValuePair;
/**
* Base class for implementing least squares optimizers.
* <p>This base class handles the boilerplate methods associated to thresholds
* settings, jacobian and error estimation.</p>
* @version $Revision: 786466 $ $Date: 2009-06-19 08:03:14 -0400 (Fri, 19 Jun 2009) $
* @since 1.2
*
*/
public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMultivariateVectorialOptimizer {
/** Default maximal number of iterations allowed. */
public static final int DEFAULT_MAX_ITERATIONS = 100;
/** Maximal number of iterations allowed. */
private int maxIterations;
/** Number of iterations already performed. */
private int iterations;
/** Maximal number of evaluations allowed. */
private int maxEvaluations;
/** Number of evaluations already performed. */
private int objectiveEvaluations;
/** Number of jacobian evaluations. */
private int jacobianEvaluations;
/** Convergence checker. */
protected VectorialConvergenceChecker checker;
/**
* Jacobian matrix.
* <p>This matrix is in canonical form just after the calls to
* {@link #updateJacobian()}, but may be modified by the solver
* in the derived class (the {@link LevenbergMarquardtOptimizer
* Levenberg-Marquardt optimizer} does this).</p>
*/
protected double[][] jacobian;
/** Number of columns of the jacobian matrix. */
protected int cols;
/** Number of rows of the jacobian matrix. */
protected int rows;
/** Objective function. */
private DifferentiableMultivariateVectorialFunction f;
/** Objective function derivatives. */
private MultivariateMatrixFunction jF;
/** Target value for the objective functions at optimum. */
protected double[] target;
/** Weight for the least squares cost computation. */
protected double[] weights;
/** Current point. */
protected double[] point;
/** Current objective function value. */
protected double[] objective;
/** Current residuals. */
protected double[] residuals;
/** Cost value (square root of the sum of the residuals). */
protected double cost;
/** Simple constructor with default settings.
* <p>The convergence check is set to a {@link SimpleVectorialValueChecker}
* and the maximal number of evaluation is set to its default value.</p>
*/
protected AbstractLeastSquaresOptimizer() {
setConvergenceChecker(new SimpleVectorialValueChecker());
setMaxIterations(DEFAULT_MAX_ITERATIONS);
setMaxEvaluations(Integer.MAX_VALUE);
}
/** {@inheritDoc} */
public void setMaxIterations(int maxIterations) {
this.maxIterations = maxIterations;
}
/** {@inheritDoc} */
public int getMaxIterations() {
return maxIterations;
}
/** {@inheritDoc} */
public int getIterations() {
return iterations;
}
/** {@inheritDoc} */
public void setMaxEvaluations(int maxEvaluations) {
this.maxEvaluations = maxEvaluations;
}
/** {@inheritDoc} */
public int getMaxEvaluations() {
return maxEvaluations;
}
/** {@inheritDoc} */
public int getEvaluations() {
return objectiveEvaluations;
}
/** {@inheritDoc} */
public int getJacobianEvaluations() {
return jacobianEvaluations;
}
/** {@inheritDoc} */
public void setConvergenceChecker(VectorialConvergenceChecker checker) {
this.checker = checker;
}
/** {@inheritDoc} */
public VectorialConvergenceChecker getConvergenceChecker() {
return checker;
}
/** Increment the iterations counter by 1.
* @exception OptimizationException if the maximal number
* of iterations is exceeded
*/
protected void incrementIterationsCounter()
throws OptimizationException {
if (++iterations > maxIterations) {
throw new OptimizationException(new MaxIterationsExceededException(maxIterations));
}
}
/**
* Update the jacobian matrix.
* @exception FunctionEvaluationException if the function jacobian
* cannot be evaluated or its dimension doesn't match problem dimension
*/
protected void updateJacobian() throws FunctionEvaluationException {
++jacobianEvaluations;
jacobian = jF.value(point);
if (jacobian.length != rows) {
throw new FunctionEvaluationException(point, "dimension mismatch {0} != {1}",
jacobian.length, rows);
}
for (int i = 0; i < rows; i++) {
final double[] ji = jacobian[i];
final double factor = -Math.sqrt(weights[i]);
for (int j = 0; j < cols; ++j) {
ji[j] *= factor;
}
}
}
/**
* Update the residuals array and cost function value.
* @exception FunctionEvaluationException if the function cannot be evaluated
* or its dimension doesn't match problem dimension or maximal number of
* of evaluations is exceeded
*/
protected void updateResidualsAndCost()
throws FunctionEvaluationException {
if (++objectiveEvaluations > maxEvaluations) {
throw new FunctionEvaluationException(new MaxEvaluationsExceededException(maxEvaluations),
point);
}
objective = f.value(point);
if (objective.length != rows) {
throw new FunctionEvaluationException(point, "dimension mismatch {0} != {1}",
objective.length, rows);
}
cost = 0;
for (int i = 0, index = 0; i < rows; i++, index += cols) {
final double residual = target[i] - objective[i];
residuals[i] = residual;
cost += weights[i] * residual * residual;
}
cost = Math.sqrt(cost);
}
/**
* Get the Root Mean Square value.
* Get the Root Mean Square value, i.e. the root of the arithmetic
* mean of the square of all weighted residuals. This is related to the
* criterion that is minimized by the optimizer as follows: if
* <em>c</em> if the criterion, and <em>n</em> is the number of
* measurements, then the RMS is <em>sqrt (c/n)</em>.
*
* @return RMS value
*/
public double getRMS() {
double criterion = 0;
for (int i = 0; i < rows; ++i) {
final double residual = residuals[i];
criterion += weights[i] * residual * residual;
}
return Math.sqrt(criterion / rows);
}
/**
* Get the Chi-Square value.
* @return chi-square value
*/
public double getChiSquare() {
double chiSquare = 0;
for (int i = 0; i < rows; ++i) {
final double residual = residuals[i];
chiSquare += residual * residual / weights[i];
}
return chiSquare;
}
/**
* Get the covariance matrix of optimized parameters.
* @return covariance matrix
* @exception FunctionEvaluationException if the function jacobian cannot
* be evaluated
* @exception OptimizationException if the covariance matrix
* cannot be computed (singular problem)
*/
public double[][] getCovariances()
throws FunctionEvaluationException, OptimizationException {
// set up the jacobian
updateJacobian();
// compute transpose(J).J, avoiding building big intermediate matrices
double[][] jTj = new double[cols][cols];
for (int i = 0; i < cols; ++i) {
for (int j = i; j < cols; ++j) {
double sum = 0;
for (int k = 0; k < rows; ++k) {
sum += jacobian[k][i] * jacobian[k][j];
}
jTj[i][j] = sum;
jTj[j][i] = sum;
}
}
try {
// compute the covariances matrix
RealMatrix inverse =
new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
return inverse.getData();
} catch (InvalidMatrixException ime) {
throw new OptimizationException("unable to compute covariances: singular problem");
}
}
/**
* Guess the errors in optimized parameters.
* <p>Guessing is covariance-based, it only gives rough order of magnitude.</p>
* @return errors in optimized parameters
* @exception FunctionEvaluationException if the function jacobian cannot b evaluated
* @exception OptimizationException if the covariances matrix cannot be computed
* or the number of degrees of freedom is not positive (number of measurements
* lesser or equal to number of parameters)
*/
public double[] guessParametersErrors()
throws FunctionEvaluationException, OptimizationException {
if (rows <= cols) {
throw new OptimizationException(
"no degrees of freedom ({0} measurements, {1} parameters)",
rows, cols);
}
double[] errors = new double[cols];
final double c = Math.sqrt(getChiSquare() / (rows - cols));
double[][] covar = getCovariances();
for (int i = 0; i < errors.length; ++i) {
errors[i] = Math.sqrt(covar[i][i]) * c;
}
return errors;
}
/** {@inheritDoc} */
public VectorialPointValuePair optimize(final DifferentiableMultivariateVectorialFunction f,
final double[] target, final double[] weights,
final double[] startPoint)
throws FunctionEvaluationException, OptimizationException, IllegalArgumentException {
if (target.length != weights.length) {
throw new OptimizationException("dimension mismatch {0} != {1}",
target.length, weights.length);
}
// reset counters
iterations = 0;
objectiveEvaluations = 0;
jacobianEvaluations = 0;
// store least squares problem characteristics
this.f = f;
jF = f.jacobian();
this.target = target.clone();
this.weights = weights.clone();
this.point = startPoint.clone();
this.residuals = new double[target.length];
// arrays shared with the o