/*
* 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.exception.util.LocalizedFormats;
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;
import org.apache.commons.math.util.FastMath;
/**
* 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: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $
* @since 1.2
*
*/
public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMultivariateVectorialOptimizer {
/** Default maximal number of iterations allowed. */
public static final int DEFAULT_MAX_ITERATIONS = 100;
/** 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;
/**
* Target value for the objective functions at optimum.
* @since 2.1
*/
protected double[] targetValues;
/**
* Weight for the least squares cost computation.
* @since 2.1
*/
protected double[] residualsWeights;
/** Current point. */
protected double[] point;
/** Current objective function value. */
protected double[] objective;
/** Current residuals. */
protected double[] residuals;
/** Weighted Jacobian */
protected double[][] wjacobian;
/** Weighted residuals */
protected double[] wresiduals;
/** Cost value (square root of the sum of the residuals). */
protected double cost;
/** 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;
/** Objective function. */
private DifferentiableMultivariateVectorialFunction function;
/** Objective function derivatives. */
private MultivariateMatrixFunction jF;
/** 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 convergenceChecker) {
this.checker = convergenceChecker;
}
/** {@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, LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
jacobian.length, rows);
}
for (int i = 0; i < rows; i++) {
final double[] ji = jacobian[i];
double wi = FastMath.sqrt(residualsWeights[i]);
for (int j = 0; j < cols; ++j) {
ji[j] *= -1.0;
wjacobian[i][j] = ji[j]*wi;
}
}
}
/**
* 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 = function.value(point);
if (objective.length != rows) {
throw new FunctionEvaluationException(point, LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
objective.length, rows);
}
cost = 0;
int index = 0;
for (int i = 0; i < rows; i++) {
final double residual = targetValues[i] - objective[i];
residuals[i] = residual;
wresiduals[i]= residual*FastMath.sqrt(residualsWeights[i]);
cost += residualsWeights[i] * residual * residual;
index += cols;
}
cost = FastMath.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() {
return FastMath.sqrt(getChiSquare() / rows);
}
/**
* Get a Chi-Square-like value assuming the N residuals follow N
* distinct normal distributions centered on 0 and whose variances are
* the reciprocal of the weights.
* @return chi-square value
*/
public double getChiSquare() {
return cost*cost;
}
/**
* 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 += wjacobian[k][i] * wjacobian[k][j];
}
jTj[i][j] = sum;
jTj[j][i] = sum;
}
}
try {
// compute the covariance matrix
RealMatrix inverse =
new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
return inverse.getData();
} catch (InvalidMatrixException ime) {
throw new OptimizationException(LocalizedFormats.UNABLE_TO_COMPUTE_COVARIANCE_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(
LocalizedFormats.NO_DEGREES_OF_FREEDOM,
rows, cols);
}
double[] errors = new double[cols];
final double c = FastMath.sqrt(getChiSquare() / (rows - cols));
double[][] covar = getCovariances();
for (int i = 0; i < errors.length; ++i) {
errors[i] = FastMath.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(LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
target.length, weights.length);
}
// reset counters
iterations = 0;
objectiveEvaluations = 0;
jacobianEvaluations = 0;
// store least squares problem characteristics
function = f;
jF = f.jacobian();
targetValues = target.clone();
residualsWeights = weights.clone();
this.point = startPoint.clone();
this.residuals = new double[target.length];
// arrays shared with the other private methods
rows = target.length;
cols = point.length;
jacobian = new double[rows][cols];
wjacobian = new double[rows][cols];
wresiduals = new double[rows];
cost = Double.POSITIVE_INFINITY;
return doOptimize();
}
/** Perform the bulk of optimization algorithm.
* @return the point/value pair giving the optimal value for objective function
* @exception FunctionEvaluationException if the objective function throws one during
* the search
* @exception OptimizationException if the algorithm failed to converge
* @exception IllegalArgumentException if the start point dimension is wrong
*/
protected abstract VectorialPointValuePair doOptimize()
throws FunctionEvaluationException, OptimizationException, IllegalArgumentException;
}