001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 package org.apache.commons.math3.optim.nonlinear.vector.jacobian; 018 019 import org.apache.commons.math3.exception.DimensionMismatchException; 020 import org.apache.commons.math3.exception.TooManyEvaluationsException; 021 import org.apache.commons.math3.linear.ArrayRealVector; 022 import org.apache.commons.math3.linear.RealMatrix; 023 import org.apache.commons.math3.linear.DecompositionSolver; 024 import org.apache.commons.math3.linear.MatrixUtils; 025 import org.apache.commons.math3.linear.QRDecomposition; 026 import org.apache.commons.math3.linear.EigenDecomposition; 027 import org.apache.commons.math3.optim.OptimizationData; 028 import org.apache.commons.math3.optim.ConvergenceChecker; 029 import org.apache.commons.math3.optim.PointVectorValuePair; 030 import org.apache.commons.math3.optim.nonlinear.vector.Weight; 031 import org.apache.commons.math3.optim.nonlinear.vector.JacobianMultivariateVectorOptimizer; 032 import org.apache.commons.math3.util.FastMath; 033 034 /** 035 * Base class for implementing least-squares optimizers. 036 * It provides methods for error estimation. 037 * 038 * @version $Id$ 039 * @since 3.1 040 */ 041 public abstract class AbstractLeastSquaresOptimizer 042 extends JacobianMultivariateVectorOptimizer { 043 /** Square-root of the weight matrix. */ 044 private RealMatrix weightMatrixSqrt; 045 /** Cost value (square root of the sum of the residuals). */ 046 private double cost; 047 048 /** 049 * @param checker Convergence checker. 050 */ 051 protected AbstractLeastSquaresOptimizer(ConvergenceChecker<PointVectorValuePair> checker) { 052 super(checker); 053 } 054 055 /** 056 * Computes the weighted Jacobian matrix. 057 * 058 * @param params Model parameters at which to compute the Jacobian. 059 * @return the weighted Jacobian: W<sup>1/2</sup> J. 060 * @throws DimensionMismatchException if the Jacobian dimension does not 061 * match problem dimension. 062 */ 063 protected RealMatrix computeWeightedJacobian(double[] params) { 064 return weightMatrixSqrt.multiply(MatrixUtils.createRealMatrix(computeJacobian(params))); 065 } 066 067 /** 068 * Computes the cost. 069 * 070 * @param residuals Residuals. 071 * @return the cost. 072 * @see #computeResiduals(double[]) 073 */ 074 protected double computeCost(double[] residuals) { 075 final ArrayRealVector r = new ArrayRealVector(residuals); 076 return FastMath.sqrt(r.dotProduct(getWeight().operate(r))); 077 } 078 079 /** 080 * Gets the root-mean-square (RMS) value. 081 * 082 * The RMS the root of the arithmetic mean of the square of all weighted 083 * residuals. 084 * This is related to the criterion that is minimized by the optimizer 085 * as follows: If <em>c</em> if the criterion, and <em>n</em> is the 086 * number of measurements, then the RMS is <em>sqrt (c/n)</em>. 087 * 088 * @return the RMS value. 089 */ 090 public double getRMS() { 091 return FastMath.sqrt(getChiSquare() / getTargetSize()); 092 } 093 094 /** 095 * Get a Chi-Square-like value assuming the N residuals follow N 096 * distinct normal distributions centered on 0 and whose variances are 097 * the reciprocal of the weights. 098 * @return chi-square value 099 */ 100 public double getChiSquare() { 101 return cost * cost; 102 } 103 104 /** 105 * Gets the square-root of the weight matrix. 106 * 107 * @return the square-root of the weight matrix. 108 */ 109 public RealMatrix getWeightSquareRoot() { 110 return weightMatrixSqrt.copy(); 111 } 112 113 /** 114 * Sets the cost. 115 * 116 * @param cost Cost value. 117 */ 118 protected void setCost(double cost) { 119 this.cost = cost; 120 } 121 122 /** 123 * Get the covariance matrix of the optimized parameters. 124 * <br/> 125 * Note that this operation involves the inversion of the 126 * <code>J<sup>T</sup>J</code> matrix, where {@code J} is the 127 * Jacobian matrix. 128 * The {@code threshold} parameter is a way for the caller to specify 129 * that the result of this computation should be considered meaningless, 130 * and thus trigger an exception. 131 * 132 * @param params Model parameters. 133 * @param threshold Singularity threshold. 134 * @return the covariance matrix. 135 * @throws org.apache.commons.math3.linear.SingularMatrixException 136 * if the covariance matrix cannot be computed (singular problem). 137 */ 138 public double[][] computeCovariances(double[] params, 139 double threshold) { 140 // Set up the Jacobian. 141 final RealMatrix j = computeWeightedJacobian(params); 142 143 // Compute transpose(J)J. 144 final RealMatrix jTj = j.transpose().multiply(j); 145 146 // Compute the covariances matrix. 147 final DecompositionSolver solver 148 = new QRDecomposition(jTj, threshold).getSolver(); 149 return solver.getInverse().getData(); 150 } 151 152 /** 153 * Computes an estimate of the standard deviation of the parameters. The 154 * returned values are the square root of the diagonal coefficients of the 155 * covariance matrix, {@code sd(a[i]) ~= sqrt(C[i][i])}, where {@code a[i]} 156 * is the optimized value of the {@code i}-th parameter, and {@code C} is 157 * the covariance matrix. 158 * 159 * @param params Model parameters. 160 * @param covarianceSingularityThreshold Singularity threshold (see 161 * {@link #computeCovariances(double[],double) computeCovariances}). 162 * @return an estimate of the standard deviation of the optimized parameters 163 * @throws org.apache.commons.math3.linear.SingularMatrixException 164 * if the covariance matrix cannot be computed. 165 */ 166 public double[] computeSigma(double[] params, 167 double covarianceSingularityThreshold) { 168 final int nC = params.length; 169 final double[] sig = new double[nC]; 170 final double[][] cov = computeCovariances(params, covarianceSingularityThreshold); 171 for (int i = 0; i < nC; ++i) { 172 sig[i] = FastMath.sqrt(cov[i][i]); 173 } 174 return sig; 175 } 176 177 /** 178 * {@inheritDoc} 179 * 180 * @param optData Optimization data. The following data will be looked for: 181 * <ul> 182 * <li>{@link org.apache.commons.math3.optim.MaxEval}</li> 183 * <li>{@link org.apache.commons.math3.optim.InitialGuess}</li> 184 * <li>{@link org.apache.commons.math3.optim.SimpleBounds}</li> 185 * <li>{@link org.apache.commons.math3.optim.nonlinear.vector.Target}</li> 186 * <li>{@link org.apache.commons.math3.optim.nonlinear.vector.Weight}</li> 187 * <li>{@link org.apache.commons.math3.optim.nonlinear.vector.ModelFunction}</li> 188 * <li>{@link org.apache.commons.math3.optim.nonlinear.vector.ModelFunctionJacobian}</li> 189 * </ul> 190 * @return {@inheritDoc} 191 * @throws TooManyEvaluationsException if the maximal number of 192 * evaluations is exceeded. 193 * @throws DimensionMismatchException if the initial guess, target, and weight 194 * arguments have inconsistent dimensions. 195 */ 196 @Override 197 public PointVectorValuePair optimize(OptimizationData... optData) 198 throws TooManyEvaluationsException { 199 // Retrieve settings. 200 parseOptimizationData(optData); 201 // Set up base class and perform computation. 202 return super.optimize(optData); 203 } 204 205 /** 206 * Computes the residuals. 207 * The residual is the difference between the observed (target) 208 * values and the model (objective function) value. 209 * There is one residual for each element of the vector-valued 210 * function. 211 * 212 * @param objectiveValue Value of the the objective function. This is 213 * the value returned from a call to 214 * {@link #computeObjectiveValue(double[]) computeObjectiveValue} 215 * (whose array argument contains the model parameters). 216 * @return the residuals. 217 * @throws DimensionMismatchException if {@code params} has a wrong 218 * length. 219 */ 220 protected double[] computeResiduals(double[] objectiveValue) { 221 final double[] target = getTarget(); 222 if (objectiveValue.length != target.length) { 223 throw new DimensionMismatchException(target.length, 224 objectiveValue.length); 225 } 226 227 final double[] residuals = new double[target.length]; 228 for (int i = 0; i < target.length; i++) { 229 residuals[i] = target[i] - objectiveValue[i]; 230 } 231 232 return residuals; 233 } 234 235 /** 236 * Scans the list of (required and optional) optimization data that 237 * characterize the problem. 238 * If the weight matrix is specified, the {@link #weightMatrixSqrt} 239 * field is recomputed. 240 * 241 * @param optData Optimization data. The following data will be looked for: 242 * <ul> 243 * <li>{@link Weight}</li> 244 * </ul> 245 */ 246 private void parseOptimizationData(OptimizationData... optData) { 247 // The existing values (as set by the previous call) are reused if 248 // not provided in the argument list. 249 for (OptimizationData data : optData) { 250 if (data instanceof Weight) { 251 weightMatrixSqrt = squareRoot(((Weight) data).getWeight()); 252 // If more data must be parsed, this statement _must_ be 253 // changed to "continue". 254 break; 255 } 256 } 257 } 258 259 /** 260 * Computes the square-root of the weight matrix. 261 * 262 * @param m Symmetric, positive-definite (weight) matrix. 263 * @return the square-root of the weight matrix. 264 */ 265 private RealMatrix squareRoot(RealMatrix m) { 266 final EigenDecomposition dec = new EigenDecomposition(m); 267 return dec.getSquareRoot(); 268 } 269 }