/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.optimization.general;
import org.apache.commons.math3.analysis.DifferentiableMultivariateVectorFunction;
import org.apache.commons.math3.analysis.FunctionUtils;
import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
import org.apache.commons.math3.analysis.differentiation.MultivariateDifferentiableVectorFunction;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.NumberIsTooSmallException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.linear.ArrayRealVector;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.DiagonalMatrix;
import org.apache.commons.math3.linear.DecompositionSolver;
import org.apache.commons.math3.linear.MatrixUtils;
import org.apache.commons.math3.linear.QRDecomposition;
import org.apache.commons.math3.linear.EigenDecomposition;
import org.apache.commons.math3.optimization.OptimizationData;
import org.apache.commons.math3.optimization.InitialGuess;
import org.apache.commons.math3.optimization.Target;
import org.apache.commons.math3.optimization.Weight;
import org.apache.commons.math3.optimization.ConvergenceChecker;
import org.apache.commons.math3.optimization.DifferentiableMultivariateVectorOptimizer;
import org.apache.commons.math3.optimization.PointVectorValuePair;
import org.apache.commons.math3.optimization.direct.BaseAbstractMultivariateVectorOptimizer;
import org.apache.commons.math3.util.FastMath;
/**
* Base class for implementing least squares optimizers.
* It handles the boilerplate methods associated to thresholds settings,
* Jacobian and error estimation.
* <br/>
* This class constructs the Jacobian matrix of the function argument in method
* {@link BaseAbstractMultivariateVectorOptimizer#optimize(int,
* org.apache.commons.math3.analysis.MultivariateVectorFunction,OptimizationData[])
* optimize} and assumes that the rows of that matrix iterate on the model
* functions while the columns iterate on the parameters; thus, the numbers
* of rows is equal to the dimension of the
* {@link org.apache.commons.math3.optimization.Target Target} while
* the number of columns is equal to the dimension of the
* {@link org.apache.commons.math3.optimization.InitialGuess InitialGuess}.
*
* @deprecated As of 3.1 (to be removed in 4.0).
* @since 1.2
*/
@Deprecated
public abstract class AbstractLeastSquaresOptimizer
extends BaseAbstractMultivariateVectorOptimizer<DifferentiableMultivariateVectorFunction>
implements DifferentiableMultivariateVectorOptimizer {
/**
* Singularity threshold (cf. {@link #getCovariances(double)}).
* @deprecated As of 3.1.
*/
@Deprecated
private static final double DEFAULT_SINGULARITY_THRESHOLD = 1e-14;
/**
* Jacobian matrix of the weighted residuals.
* 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).
* @deprecated As of 3.1. To be removed in 4.0. Please use
* {@link #computeWeightedJacobian(double[])} instead.
*/
@Deprecated
protected double[][] weightedResidualJacobian;
/** Number of columns of the jacobian matrix.
* @deprecated As of 3.1.
*/
@Deprecated
protected int cols;
/** Number of rows of the jacobian matrix.
* @deprecated As of 3.1.
*/
@Deprecated
protected int rows;
/** Current point.
* @deprecated As of 3.1.
*/
@Deprecated
protected double[] point;
/** Current objective function value.
* @deprecated As of 3.1.
*/
@Deprecated
protected double[] objective;
/** Weighted residuals
* @deprecated As of 3.1.
*/
@Deprecated
protected double[] weightedResiduals;
/** Cost value (square root of the sum of the residuals).
* @deprecated As of 3.1. Field to become "private" in 4.0.
* Please use {@link #setCost(double)}.
*/
@Deprecated
protected double cost;
/** Objective function derivatives. */
private MultivariateDifferentiableVectorFunction jF;
/** Number of evaluations of the Jacobian. */
private int jacobianEvaluations;
/** Square-root of the weight matrix. */
private RealMatrix weightMatrixSqrt;
/**
* Simple constructor with default settings.
* The convergence check is set to a {@link
* org.apache.commons.math3.optimization.SimpleVectorValueChecker}.
* @deprecated See {@link org.apache.commons.math3.optimization.SimpleValueChecker#SimpleValueChecker()}
*/
@Deprecated
protected AbstractLeastSquaresOptimizer() {}
/**
* @param checker Convergence checker.
*/
protected AbstractLeastSquaresOptimizer(ConvergenceChecker<PointVectorValuePair> checker) {
super(checker);
}
/**
* @return the number of evaluations of the Jacobian function.
*/
public int getJacobianEvaluations() {
return jacobianEvaluations;
}
/**
* Update the jacobian matrix.
*
* @throws DimensionMismatchException if the Jacobian dimension does not
* match problem dimension.
* @deprecated As of 3.1. Please use {@link #computeWeightedJacobian(double[])}
* instead.
*/
@Deprecated
protected void updateJacobian() {
final RealMatrix weightedJacobian = computeWeightedJacobian(point);
weightedResidualJacobian = weightedJacobian.scalarMultiply(-1).getData();
}
/**
* Computes the Jacobian matrix.
*
* @param params Model parameters at which to compute the Jacobian.
* @return the weighted Jacobian: W<sup>1/2</sup> J.
* @throws DimensionMismatchException if the Jacobian dimension does not
* match problem dimension.
* @since 3.1
*/
protected RealMatrix computeWeightedJacobian(double[] params) {
++jacobianEvaluations;
final DerivativeStructure[] dsPoint = new DerivativeStructure[params.length];
final int nC = params.length;
for (int i = 0; i < nC; ++i) {
dsPoint[i] = new DerivativeStructure(nC, 1, i, params[i]);
}
final DerivativeStructure[] dsValue = jF.value(dsPoint);
final int nR = getTarget().length;
if (dsValue.length != nR) {
throw new DimensionMismatchException(dsValue.length, nR);
}
final double[][] jacobianData = new double[nR][nC];
for (int i = 0; i < nR; ++i) {
int[] orders = new int[nC];
for (int j = 0; j < nC; ++j) {
orders[j] = 1;
jacobianData[i][j] = dsValue[i].getPartialDerivative(orders);
orders[j] = 0;
}
}
return weightMatrixSqrt.multiply(MatrixUtils.createRealMatrix(jacobianData));
}
/**
* Update the residuals array and cost function value.
* @throws DimensionMismatchException if the dimension does not match the
* problem dimension.
* @throws org.apache.commons.math3.exception.TooManyEvaluationsException
* if the maximal number of evaluations is exceeded.
* @deprecated As of 3.1. Please use {@link #computeResiduals(double[])},
* {@link #computeObjectiveValue(double[])}, {@link #computeCost(double[])}
* and {@link #setCost(double)} instead.
*/
@Deprecated
protected void updateResidualsAndCost() {
objective = computeObjectiveValue(point);
final double[] res = computeResiduals(objective);
// Compute cost.
cost = computeCost(res);
// Compute weighted residuals.
final ArrayRealVector residuals = new ArrayRealVector(res);
weightedResiduals = weightMatrixSqrt.operate(residuals).toArray();
}
/**
* Computes the cost.
*
* @param residuals Residuals.
* @return the cost.
* @see #computeResiduals(double[])
* @since 3.1
*/
protected double computeCost(double[] residuals) {
final ArrayRealVector r = new ArrayRealVector(residuals);
return FastMath.sqrt(r.dotProduct(getWeight().operate(r)));
}
/**
* 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;
}
/**
* Gets the square-root of the weight matrix.
*
* @return the square-root of the weight matrix.
* @since 3.1
*/
public RealMatrix getWeightSquareRoot() {
return weightMatrixSqrt.copy();
}
/**
* Sets the cost.
*
* @param cost Cost value.
* @since 3.1
*/
protected void setCost(double cost) {
this.cost = cost;
}
/**
* Get the covariance matrix of the optimized parameters.
*
* @return the covariance matrix.
* @throws org.apache.commons.math3.linear.SingularMatrixException
* if the covariance matrix cannot be computed (singular problem).
* @see #getCovariances(double)
* @deprecated As of 3.1. Please use {@link #computeCovariances(double[],double)}
* instead.
*/
@Deprecated
public double[][] getCovariances() {
return getCovariances(DEFAULT_SINGULARITY_THRESHOLD);
}
/**
* Get the covariance matrix of the optimized parameters.
* <br/>
* Note that this operation involves the inversion of the
* <code>J<sup>T</sup>J</code> matrix, where {@code J} is the
* Jacobian matrix.
* The {@code threshold} parameter is a way for the caller to specify
* that the result of this computation should be considered meaningless,
* and thus trigger an exception.
*
* @param threshold Singularity threshold.
* @return the covariance matrix.
* @throws org.apache.commons.math3.linear.SingularMatrixException
* if the covariance matrix cannot be computed (singular problem).
* @deprecated As of 3.1. Please use {@link #computeCovariances(double[],double)}
* instead.
*/
@Deprecated
public double[][] getCovariances(double threshold) {
return computeCovariances(point, threshold);
}
/**
* Get the covariance matrix of the optimized parameters.
* <br/>
* Note that this operation involves the inversion of the
* <code>J<sup>T</sup>J</code> matrix, where {@code J} is the
* Jacobian matrix.
* The {@code threshold} parameter is a way for the caller to specify
* that the result of this computation should be considered meaningless,
* and thus trigger an exception.
*
* @param params Model parameters.
* @param threshold Singularity threshold.
* @return the covariance matrix.
* @throws org.apache.commons.math3.linear.SingularMatrixException
* if the covariance matrix cannot be computed (singular problem).
* @since 3.1
*/
public double[][] computeCovariances(double[] params,
double threshold) {
// Set up the Jacobian.
final RealMatrix j = computeWeightedJacobian(params);
// Compute transpose(J)J.
final RealMatrix jTj = j.transpose().multiply(j);
// Compute the covariances matrix.
final DecompositionSolver solver
= new QRDecomposition(jTj, threshold).getSolver();
return solver.getInverse().getData();
}
/**
* <p>
* Returns an estimate of the standard deviation of each parameter. The
* returned values are the so-called (asymptotic) standard errors on the
* parameters, defined as {@code sd(a[i]) = sqrt(S / (n - m) * C[i][i])},
* where {@code a[i]} is the optimized value of the {@code i}-th parameter,
* {@code S} is the minimized value of the sum of squares objective function
* (as returned by {@link #getChiSquare()}), {@code n} is the number of
* observations, {@code m} is the number of parameters and {@code C} is the
* covariance matrix.
* </p>
* <p>
* See also
* <a href="http://en.wikipedia.org/wiki/Least_squares">Wikipedia</a>,
* or
* <a href="http://mathworld.wolfram.com/LeastSquaresFitting.html">MathWorld</a>,
* equations (34) and (35) for a particular case.
* </p>
*
* @return an estimate of the standard deviation of the optimized parameters
* @throws org.apache.commons.math3.linear.SingularMatrixException
* if the covariance matrix cannot be computed.
* @throws NumberIsTooSmallException if the number of degrees of freedom is not
* positive, i.e. the number of measurements is less or equal to the number of
* parameters.
* @deprecated as of version 3.1, {@link #computeSigma(double[],double)} should be used
* instead. It should be emphasized that {@code guessParametersErrors} and
* {@code computeSigma} are <em>not</em> strictly equivalent.
*/
@Deprecated
public double[] guessParametersErrors() {
if (rows <= cols) {
throw new NumberIsTooSmallException(LocalizedFormats.NO_DEGREES_OF_FREEDOM,
rows, cols, false);
}
double[] errors = new double[cols];
final double c = FastMath.sqrt(getChiSquare() / (rows - cols));
double[][] covar = computeCovariances(point, 1e-14);
for (int i = 0; i < errors.length; ++i) {
errors[i] = FastMath.sqrt(covar[i][i]) * c;
}
return errors;
}
/**
* Computes an estimate of the standard deviation of the parameters. The
* returned values are the square root of the diagonal coefficients of the
* covariance matrix, {@code sd(a[i]) ~= sqrt(C[i][i])}, where {@code a[i]}
* is the optimized value of the {@code i}-th parameter, and {@code C} is
* the covariance matrix.
*
* @param params Model parameters.
* @param covarianceSingularityThreshold Singularity threshold (see
* {@link #computeCovariances(double[],double) computeCovariances}).
* @return an estimate of the standard deviation of the optimized parameters
* @throws org.apache.commons.math3.linear.SingularMatrixException
* if the covariance matrix cannot be computed.
* @since 3.1
*/
public double[] computeSigma(double[] params,
double covarianceSingularityThreshold) {
final int nC = params.length;
final double[] sig = new double[nC];
final double[][] cov = computeCovariances(params, covarianceSingularityThreshold);
for (int i = 0; i < nC; ++i) {
sig[i] = FastMath.sqrt(cov[i][i]);
}
return sig;
}
/** {@inheritDoc}
* @deprecated As of 3.1. Please use
* {@link BaseAbstractMultivariateVectorOptimizer#optimize(int,
* org.apache.commons.math3.analysis.MultivariateVectorFunction,OptimizationData[])
* optimize(int,MultivariateDifferentiableVectorFunction,OptimizationData...)}
* instead.
*/
@Override
@Deprecated
public PointVectorValuePair optimize(int maxEval,
final DifferentiableMultivariateVectorFunction f,
final double[] target, final double[] weights,
final double[] startPoint) {
return optimizeInternal(maxEval,
FunctionUtils.toMultivariateDifferentiableVectorFunction(f),
new Target(target),
new Weight(weights),
new InitialGuess(startPoint));
}
/**
* Optimize an objective function.
* Optimization is considered to be a weighted least-squares minimization.
* The cost function to be minimized is
* <code>∑weight<sub>i</sub>(objective<sub>i</sub> - target<sub>i</sub>)<sup>2</sup></code>
*
* @param f Objective function.
* @param target Target value for the objective functions at optimum.
* @param weights Weights for the least squares cost computation.
* @param startPoint Start point for optimization.
* @return the point/value pair giving the optimal value for objective
* function.
* @param maxEval Maximum number of function evaluations.
* @throws org.apache.commons.math3.exception.DimensionMismatchException
* if the start point dimension is wrong.
* @throws org.apache.commons.math3.exception.TooManyEvaluationsException
* if the maximal number of evaluations is exceeded.
* @throws org.apache.commons.math3.exception.NullArgumentException if
* any argument is {@code null}.
* @deprecated As of 3.1. Please use
* {@link BaseAbstractMultivariateVectorOptimizer#optimize(int,
* org.apache.commons.math3.analysis.MultivariateVectorFunction,OptimizationData[])
* optimize(int,MultivariateDifferentiableVectorFunction,OptimizationData...)}
* instead.
*/
@Deprecated
public PointVectorValuePair optimize(final int maxEval,
final MultivariateDifferentiableVectorFunction f,
final double[] target, final double[] weights,
final double[] startPoint) {
return optimizeInternal(maxEval, f,
new Target(target),
new Weight(weights),
new InitialGuess(startPoint));
}
/**
* Optimize an objective function.
* Optimization is considered to be a weighted least-squares minimization.
* The cost function to be minimized is
* <code>∑weight<sub>i</sub>(objective<sub>i</sub> - target<sub>i</sub>)<sup>2</sup></code>
*
* @param maxEval Allowed number of evaluations of the objective function.
* @param f Objective function.
* @param optData Optimization data. The following data will be looked for:
* <ul>
* <li>{@link Target}</li>
* <li>{@link Weight}</li>
* <li>{@link InitialGuess}</li>
* </ul>
* @return the point/value pair giving the optimal value of the objective
* function.
* @throws org.apache.commons.math3.exception.TooManyEvaluationsException if
* the maximal number of evaluations is exceeded.
* @throws DimensionMismatchException if the target, and weight arguments
* have inconsistent dimensions.
* @see BaseAbstractMultivariateVectorOptimizer#optimizeInternal(int,
* org.apache.commons.math3.analysis.MultivariateVectorFunction,OptimizationData[])
* @since 3.1
* @deprecated As of 3.1. Override is necessary only until this class's generic
* argument is changed to {@code MultivariateDifferentiableVectorFunction}.
*/
@Deprecated
protected PointVectorValuePair optimizeInternal(final int maxEval,
final MultivariateDifferentiableVectorFunction f,
OptimizationData... optData) {
// XXX Conversion will be removed when the generic argument of the
// base class becomes "MultivariateDifferentiableVectorFunction".
return super.optimizeInternal(maxEval, FunctionUtils.toDifferentiableMultivariateVectorFunction(f), optData);
}
/** {@inheritDoc} */
@Override
protected void setUp() {
super.setUp();
// Reset counter.
jacobianEvaluations = 0;
// Square-root of the weight matrix.
weightMatrixSqrt = squareRoot(getWeight());
// Store least squares problem characteristics.
// XXX The conversion won't be necessary when the generic argument of
// the base class becomes "MultivariateDifferentiableVectorFunction".
// XXX "jF" is not strictly necessary anymore but is currently more
// efficient than converting the value returned from "getObjectiveFunction()"
// every time it is used.
jF = FunctionUtils.toMultivariateDifferentiableVectorFunction((DifferentiableMultivariateVectorFunction) getObjectiveFunction());
// Arrays shared with "private" and "protected" methods.
point = getStartPoint();
rows = getTarget().length;
cols = point.length;
}
/**
* Computes the residuals.
* The residual is the difference between the observed (target)
* values and the model (objective function) value.
* There is one residual for each element of the vector-valued
* function.
*
* @param objectiveValue Value of the the objective function. This is
* the value returned from a call to
* {@link #computeObjectiveValue(double[]) computeObjectiveValue}
* (whose array argument contains the model parameters).
* @return the residuals.
* @throws DimensionMismatchException if {@code params} has a wrong
* length.
* @since 3.1
*/
protected double[] computeResiduals(double[] objectiveValue) {
final double[] target = getTarget();
if (objectiveValue.length != target.length) {
throw new DimensionMismatchException(target.length,
objectiveValue.length);
}
final double[] residuals = new double[target.length];
for (int i = 0; i < target.length; i++) {
residuals[i] = target[i] - objectiveValue[i];
}
return residuals;
}
/**
* Computes the square-root of the weight matrix.
*
* @param m Symmetric, positive-definite (weight) matrix.
* @return the square-root of the weight matrix.
*/
private RealMatrix squareRoot(RealMatrix m) {
if (m instanceof DiagonalMatrix) {
final int dim = m.getRowDimension();
final RealMatrix sqrtM = new DiagonalMatrix(dim);
for (int i = 0; i < dim; i++) {
sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
}
return sqrtM;
} else {
final EigenDecomposition dec = new EigenDecomposition(m);
return dec.getSquareRoot();
}
}
}