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 }