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