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 }