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