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 018 package org.apache.commons.math.optimization.general; 019 020 import org.apache.commons.math.FunctionEvaluationException; 021 import org.apache.commons.math.MaxEvaluationsExceededException; 022 import org.apache.commons.math.MaxIterationsExceededException; 023 import org.apache.commons.math.analysis.DifferentiableMultivariateVectorialFunction; 024 import org.apache.commons.math.analysis.MultivariateMatrixFunction; 025 import org.apache.commons.math.linear.InvalidMatrixException; 026 import org.apache.commons.math.linear.LUDecompositionImpl; 027 import org.apache.commons.math.linear.MatrixUtils; 028 import org.apache.commons.math.linear.RealMatrix; 029 import org.apache.commons.math.optimization.OptimizationException; 030 import org.apache.commons.math.optimization.SimpleVectorialValueChecker; 031 import org.apache.commons.math.optimization.VectorialConvergenceChecker; 032 import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer; 033 import org.apache.commons.math.optimization.VectorialPointValuePair; 034 035 /** 036 * Base class for implementing least squares optimizers. 037 * <p>This base class handles the boilerplate methods associated to thresholds 038 * settings, jacobian and error estimation.</p> 039 * @version $Revision: 786466 $ $Date: 2009-06-19 08:03:14 -0400 (Fri, 19 Jun 2009) $ 040 * @since 1.2 041 * 042 */ 043 public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMultivariateVectorialOptimizer { 044 045 /** Default maximal number of iterations allowed. */ 046 public static final int DEFAULT_MAX_ITERATIONS = 100; 047 048 /** Maximal number of iterations allowed. */ 049 private int maxIterations; 050 051 /** Number of iterations already performed. */ 052 private int iterations; 053 054 /** Maximal number of evaluations allowed. */ 055 private int maxEvaluations; 056 057 /** Number of evaluations already performed. */ 058 private int objectiveEvaluations; 059 060 /** Number of jacobian evaluations. */ 061 private int jacobianEvaluations; 062 063 /** Convergence checker. */ 064 protected VectorialConvergenceChecker checker; 065 066 /** 067 * Jacobian matrix. 068 * <p>This matrix is in canonical form just after the calls to 069 * {@link #updateJacobian()}, but may be modified by the solver 070 * in the derived class (the {@link LevenbergMarquardtOptimizer 071 * Levenberg-Marquardt optimizer} does this).</p> 072 */ 073 protected double[][] jacobian; 074 075 /** Number of columns of the jacobian matrix. */ 076 protected int cols; 077 078 /** Number of rows of the jacobian matrix. */ 079 protected int rows; 080 081 /** Objective function. */ 082 private DifferentiableMultivariateVectorialFunction f; 083 084 /** Objective function derivatives. */ 085 private MultivariateMatrixFunction jF; 086 087 /** Target value for the objective functions at optimum. */ 088 protected double[] target; 089 090 /** Weight for the least squares cost computation. */ 091 protected double[] weights; 092 093 /** Current point. */ 094 protected double[] point; 095 096 /** Current objective function value. */ 097 protected double[] objective; 098 099 /** Current residuals. */ 100 protected double[] residuals; 101 102 /** Cost value (square root of the sum of the residuals). */ 103 protected double cost; 104 105 /** Simple constructor with default settings. 106 * <p>The convergence check is set to a {@link SimpleVectorialValueChecker} 107 * and the maximal number of evaluation is set to its default value.</p> 108 */ 109 protected AbstractLeastSquaresOptimizer() { 110 setConvergenceChecker(new SimpleVectorialValueChecker()); 111 setMaxIterations(DEFAULT_MAX_ITERATIONS); 112 setMaxEvaluations(Integer.MAX_VALUE); 113 } 114 115 /** {@inheritDoc} */ 116 public void setMaxIterations(int maxIterations) { 117 this.maxIterations = maxIterations; 118 } 119 120 /** {@inheritDoc} */ 121 public int getMaxIterations() { 122 return maxIterations; 123 } 124 125 /** {@inheritDoc} */ 126 public int getIterations() { 127 return iterations; 128 } 129 130 /** {@inheritDoc} */ 131 public void setMaxEvaluations(int maxEvaluations) { 132 this.maxEvaluations = maxEvaluations; 133 } 134 135 /** {@inheritDoc} */ 136 public int getMaxEvaluations() { 137 return maxEvaluations; 138 } 139 140 /** {@inheritDoc} */ 141 public int getEvaluations() { 142 return objectiveEvaluations; 143 } 144 145 /** {@inheritDoc} */ 146 public int getJacobianEvaluations() { 147 return jacobianEvaluations; 148 } 149 150 /** {@inheritDoc} */ 151 public void setConvergenceChecker(VectorialConvergenceChecker checker) { 152 this.checker = checker; 153 } 154 155 /** {@inheritDoc} */ 156 public VectorialConvergenceChecker getConvergenceChecker() { 157 return checker; 158 } 159 160 /** Increment the iterations counter by 1. 161 * @exception OptimizationException if the maximal number 162 * of iterations is exceeded 163 */ 164 protected void incrementIterationsCounter() 165 throws OptimizationException { 166 if (++iterations > maxIterations) { 167 throw new OptimizationException(new MaxIterationsExceededException(maxIterations)); 168 } 169 } 170 171 /** 172 * Update the jacobian matrix. 173 * @exception FunctionEvaluationException if the function jacobian 174 * cannot be evaluated or its dimension doesn't match problem dimension 175 */ 176 protected void updateJacobian() throws FunctionEvaluationException { 177 ++jacobianEvaluations; 178 jacobian = jF.value(point); 179 if (jacobian.length != rows) { 180 throw new FunctionEvaluationException(point, "dimension mismatch {0} != {1}", 181 jacobian.length, rows); 182 } 183 for (int i = 0; i < rows; i++) { 184 final double[] ji = jacobian[i]; 185 final double factor = -Math.sqrt(weights[i]); 186 for (int j = 0; j < cols; ++j) { 187 ji[j] *= factor; 188 } 189 } 190 } 191 192 /** 193 * Update the residuals array and cost function value. 194 * @exception FunctionEvaluationException if the function cannot be evaluated 195 * or its dimension doesn't match problem dimension or maximal number of 196 * of evaluations is exceeded 197 */ 198 protected void updateResidualsAndCost() 199 throws FunctionEvaluationException { 200 201 if (++objectiveEvaluations > maxEvaluations) { 202 throw new FunctionEvaluationException(new MaxEvaluationsExceededException(maxEvaluations), 203 point); 204 } 205 objective = f.value(point); 206 if (objective.length != rows) { 207 throw new FunctionEvaluationException(point, "dimension mismatch {0} != {1}", 208 objective.length, rows); 209 } 210 cost = 0; 211 for (int i = 0, index = 0; i < rows; i++, index += cols) { 212 final double residual = target[i] - objective[i]; 213 residuals[i] = residual; 214 cost += weights[i] * residual * residual; 215 } 216 cost = Math.sqrt(cost); 217 218 } 219 220 /** 221 * Get the Root Mean Square value. 222 * Get the Root Mean Square value, i.e. the root of the arithmetic 223 * mean of the square of all weighted residuals. This is related to the 224 * criterion that is minimized by the optimizer as follows: if 225 * <em>c</em> if the criterion, and <em>n</em> is the number of 226 * measurements, then the RMS is <em>sqrt (c/n)</em>. 227 * 228 * @return RMS value 229 */ 230 public double getRMS() { 231 double criterion = 0; 232 for (int i = 0; i < rows; ++i) { 233 final double residual = residuals[i]; 234 criterion += weights[i] * residual * residual; 235 } 236 return Math.sqrt(criterion / rows); 237 } 238 239 /** 240 * Get the Chi-Square value. 241 * @return chi-square value 242 */ 243 public double getChiSquare() { 244 double chiSquare = 0; 245 for (int i = 0; i < rows; ++i) { 246 final double residual = residuals[i]; 247 chiSquare += residual * residual / weights[i]; 248 } 249 return chiSquare; 250 } 251 252 /** 253 * Get the covariance matrix of optimized parameters. 254 * @return covariance matrix 255 * @exception FunctionEvaluationException if the function jacobian cannot 256 * be evaluated 257 * @exception OptimizationException if the covariance matrix 258 * cannot be computed (singular problem) 259 */ 260 public double[][] getCovariances() 261 throws FunctionEvaluationException, OptimizationException { 262 263 // set up the jacobian 264 updateJacobian(); 265 266 // compute transpose(J).J, avoiding building big intermediate matrices 267 double[][] jTj = new double[cols][cols]; 268 for (int i = 0; i < cols; ++i) { 269 for (int j = i; j < cols; ++j) { 270 double sum = 0; 271 for (int k = 0; k < rows; ++k) { 272 sum += jacobian[k][i] * jacobian[k][j]; 273 } 274 jTj[i][j] = sum; 275 jTj[j][i] = sum; 276 } 277 } 278 279 try { 280 // compute the covariances matrix 281 RealMatrix inverse = 282 new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse(); 283 return inverse.getData(); 284 } catch (InvalidMatrixException ime) { 285 throw new OptimizationException("unable to compute covariances: singular problem"); 286 } 287 288 } 289 290 /** 291 * Guess the errors in optimized parameters. 292 * <p>Guessing is covariance-based, it only gives rough order of magnitude.</p> 293 * @return errors in optimized parameters 294 * @exception FunctionEvaluationException if the function jacobian cannot b evaluated 295 * @exception OptimizationException if the covariances matrix cannot be computed 296 * or the number of degrees of freedom is not positive (number of measurements 297 * lesser or equal to number of parameters) 298 */ 299 public double[] guessParametersErrors() 300 throws FunctionEvaluationException, OptimizationException { 301 if (rows <= cols) { 302 throw new OptimizationException( 303 "no degrees of freedom ({0} measurements, {1} parameters)", 304 rows, cols); 305 } 306 double[] errors = new double[cols]; 307 final double c = Math.sqrt(getChiSquare() / (rows - cols)); 308 double[][] covar = getCovariances(); 309 for (int i = 0; i < errors.length; ++i) { 310 errors[i] = Math.sqrt(covar[i][i]) * c; 311 } 312 return errors; 313 } 314 315 /** {@inheritDoc} */ 316 public VectorialPointValuePair optimize(final DifferentiableMultivariateVectorialFunction f, 317 final double[] target, final double[] weights, 318 final double[] startPoint) 319 throws FunctionEvaluationException, OptimizationException, IllegalArgumentException { 320 321 if (target.length != weights.length) { 322 throw new OptimizationException("dimension mismatch {0} != {1}", 323 target.length, weights.length); 324 } 325 326 // reset counters 327 iterations = 0; 328 objectiveEvaluations = 0; 329 jacobianEvaluations = 0; 330 331 // store least squares problem characteristics 332 this.f = f; 333 jF = f.jacobian(); 334 this.target = target.clone(); 335 this.weights = weights.clone(); 336 this.point = startPoint.clone(); 337 this.residuals = new double[target.length]; 338 339 // arrays shared with the other private methods 340 rows = target.length; 341 cols = point.length; 342 jacobian = new double[rows][cols]; 343 344 cost = Double.POSITIVE_INFINITY; 345 346 return doOptimize(); 347 348 } 349 350 /** Perform the bulk of optimization algorithm. 351 * @return the point/value pair giving the optimal value for objective function 352 * @exception FunctionEvaluationException if the objective function throws one during 353 * the search 354 * @exception OptimizationException if the algorithm failed to converge 355 * @exception IllegalArgumentException if the start point dimension is wrong 356 */ 357 abstract protected VectorialPointValuePair doOptimize() 358 throws FunctionEvaluationException, OptimizationException, IllegalArgumentException; 359 360 }