diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/AbstractLeastSquaresOptimizer.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/AbstractLeastSquaresOptimizer.java | 281 |
1 files changed, 281 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/AbstractLeastSquaresOptimizer.java b/src/main/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/AbstractLeastSquaresOptimizer.java new file mode 100644 index 0000000..67682eb --- /dev/null +++ b/src/main/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/AbstractLeastSquaresOptimizer.java @@ -0,0 +1,281 @@ +/* + * 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.optim.nonlinear.vector.jacobian; + +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.TooManyEvaluationsException; +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.optim.OptimizationData; +import org.apache.commons.math3.optim.ConvergenceChecker; +import org.apache.commons.math3.optim.PointVectorValuePair; +import org.apache.commons.math3.optim.nonlinear.vector.Weight; +import org.apache.commons.math3.optim.nonlinear.vector.JacobianMultivariateVectorOptimizer; +import org.apache.commons.math3.util.FastMath; + +/** + * Base class for implementing least-squares optimizers. + * It provides methods for error estimation. + * + * @since 3.1 + * @deprecated All classes and interfaces in this package are deprecated. + * The optimizers that were provided here were moved to the + * {@link org.apache.commons.math3.fitting.leastsquares} package + * (cf. MATH-1008). + */ +@Deprecated +public abstract class AbstractLeastSquaresOptimizer + extends JacobianMultivariateVectorOptimizer { + /** Square-root of the weight matrix. */ + private RealMatrix weightMatrixSqrt; + /** Cost value (square root of the sum of the residuals). */ + private double cost; + + /** + * @param checker Convergence checker. + */ + protected AbstractLeastSquaresOptimizer(ConvergenceChecker<PointVectorValuePair> checker) { + super(checker); + } + + /** + * Computes the weighted 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. + */ + protected RealMatrix computeWeightedJacobian(double[] params) { + return weightMatrixSqrt.multiply(MatrixUtils.createRealMatrix(computeJacobian(params))); + } + + /** + * Computes the cost. + * + * @param residuals Residuals. + * @return the cost. + * @see #computeResiduals(double[]) + */ + protected double computeCost(double[] residuals) { + final ArrayRealVector r = new ArrayRealVector(residuals); + return FastMath.sqrt(r.dotProduct(getWeight().operate(r))); + } + + /** + * Gets the root-mean-square (RMS) value. + * + * The RMS 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 the RMS value. + */ + public double getRMS() { + return FastMath.sqrt(getChiSquare() / getTargetSize()); + } + + /** + * 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. + */ + public RealMatrix getWeightSquareRoot() { + return weightMatrixSqrt.copy(); + } + + /** + * Sets the cost. + * + * @param cost Cost value. + */ + protected void setCost(double cost) { + this.cost = cost; + } + + /** + * 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). + */ + 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(); + } + + /** + * 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. + */ + 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} + * + * @param optData Optimization data. In addition to those documented in + * {@link JacobianMultivariateVectorOptimizer#parseOptimizationData(OptimizationData[]) + * JacobianMultivariateVectorOptimizer}, this method will register the following data: + * <ul> + * <li>{@link org.apache.commons.math3.optim.nonlinear.vector.Weight}</li> + * </ul> + * @return {@inheritDoc} + * @throws TooManyEvaluationsException if the maximal number of + * evaluations is exceeded. + * @throws DimensionMismatchException if the initial guess, target, and weight + * arguments have inconsistent dimensions. + */ + @Override + public PointVectorValuePair optimize(OptimizationData... optData) + throws TooManyEvaluationsException { + // Set up base class and perform computation. + return super.optimize(optData); + } + + /** + * 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. + */ + 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; + } + + /** + * Scans the list of (required and optional) optimization data that + * characterize the problem. + * If the weight matrix is specified, the {@link #weightMatrixSqrt} + * field is recomputed. + * + * @param optData Optimization data. The following data will be looked for: + * <ul> + * <li>{@link Weight}</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 Weight) { + weightMatrixSqrt = squareRoot(((Weight) data).getWeight()); + // If more data must be parsed, this statement _must_ be + // changed to "continue". + break; + } + } + } + + /** + * 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(); + } + } +} |