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: 825919 $ $Date: 2009-10-16 10:51:55 -0400 (Fri, 16 Oct 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 /** 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 /** Target value for the objective functions at optimum. */
067 protected double[] targetValues;
068
069 /** Weight for the least squares cost computation. */
070 protected double[] residualsWeights;
071
072 /** Current point. */
073 protected double[] point;
074
075 /** Current objective function value. */
076 protected double[] objective;
077
078 /** Current residuals. */
079 protected double[] residuals;
080
081 /** Cost value (square root of the sum of the residuals). */
082 protected double cost;
083
084 /** Maximal number of iterations allowed. */
085 private int maxIterations;
086
087 /** Number of iterations already performed. */
088 private int iterations;
089
090 /** Maximal number of evaluations allowed. */
091 private int maxEvaluations;
092
093 /** Number of evaluations already performed. */
094 private int objectiveEvaluations;
095
096 /** Number of jacobian evaluations. */
097 private int jacobianEvaluations;
098
099 /** Objective function. */
100 private DifferentiableMultivariateVectorialFunction function;
101
102 /** Objective function derivatives. */
103 private MultivariateMatrixFunction jF;
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 convergenceChecker) {
152 this.checker = convergenceChecker;
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(residualsWeights[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 = function.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 int index = 0;
212 for (int i = 0; i < rows; i++) {
213 final double residual = targetValues[i] - objective[i];
214 residuals[i] = residual;
215 cost += residualsWeights[i] * residual * residual;
216 index += cols;
217 }
218 cost = Math.sqrt(cost);
219
220 }
221
222 /**
223 * Get the Root Mean Square value.
224 * Get the Root Mean Square value, i.e. the root of the arithmetic
225 * mean of the square of all weighted residuals. This is related to the
226 * criterion that is minimized by the optimizer as follows: if
227 * <em>c</em> if the criterion, and <em>n</em> is the number of
228 * measurements, then the RMS is <em>sqrt (c/n)</em>.
229 *
230 * @return RMS value
231 */
232 public double getRMS() {
233 double criterion = 0;
234 for (int i = 0; i < rows; ++i) {
235 final double residual = residuals[i];
236 criterion += residualsWeights[i] * residual * residual;
237 }
238 return Math.sqrt(criterion / rows);
239 }
240
241 /**
242 * Get the Chi-Square value.
243 * @return chi-square value
244 */
245 public double getChiSquare() {
246 double chiSquare = 0;
247 for (int i = 0; i < rows; ++i) {
248 final double residual = residuals[i];
249 chiSquare += residual * residual / residualsWeights[i];
250 }
251 return chiSquare;
252 }
253
254 /**
255 * Get the covariance matrix of optimized parameters.
256 * @return covariance matrix
257 * @exception FunctionEvaluationException if the function jacobian cannot
258 * be evaluated
259 * @exception OptimizationException if the covariance matrix
260 * cannot be computed (singular problem)
261 */
262 public double[][] getCovariances()
263 throws FunctionEvaluationException, OptimizationException {
264
265 // set up the jacobian
266 updateJacobian();
267
268 // compute transpose(J).J, avoiding building big intermediate matrices
269 double[][] jTj = new double[cols][cols];
270 for (int i = 0; i < cols; ++i) {
271 for (int j = i; j < cols; ++j) {
272 double sum = 0;
273 for (int k = 0; k < rows; ++k) {
274 sum += jacobian[k][i] * jacobian[k][j];
275 }
276 jTj[i][j] = sum;
277 jTj[j][i] = sum;
278 }
279 }
280
281 try {
282 // compute the covariances matrix
283 RealMatrix inverse =
284 new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
285 return inverse.getData();
286 } catch (InvalidMatrixException ime) {
287 throw new OptimizationException("unable to compute covariances: singular problem");
288 }
289
290 }
291
292 /**
293 * Guess the errors in optimized parameters.
294 * <p>Guessing is covariance-based, it only gives rough order of magnitude.</p>
295 * @return errors in optimized parameters
296 * @exception FunctionEvaluationException if the function jacobian cannot b evaluated
297 * @exception OptimizationException if the covariances matrix cannot be computed
298 * or the number of degrees of freedom is not positive (number of measurements
299 * lesser or equal to number of parameters)
300 */
301 public double[] guessParametersErrors()
302 throws FunctionEvaluationException, OptimizationException {
303 if (rows <= cols) {
304 throw new OptimizationException(
305 "no degrees of freedom ({0} measurements, {1} parameters)",
306 rows, cols);
307 }
308 double[] errors = new double[cols];
309 final double c = Math.sqrt(getChiSquare() / (rows - cols));
310 double[][] covar = getCovariances();
311 for (int i = 0; i < errors.length; ++i) {
312 errors[i] = Math.sqrt(covar[i][i]) * c;
313 }
314 return errors;
315 }
316
317 /** {@inheritDoc} */
318 public VectorialPointValuePair optimize(final DifferentiableMultivariateVectorialFunction f,
319 final double[] target, final double[] weights,
320 final double[] startPoint)
321 throws FunctionEvaluationException, OptimizationException, IllegalArgumentException {
322
323 if (target.length != weights.length) {
324 throw new OptimizationException("dimension mismatch {0} != {1}",
325 target.length, weights.length);
326 }
327
328 // reset counters
329 iterations = 0;
330 objectiveEvaluations = 0;
331 jacobianEvaluations = 0;
332
333 // store least squares problem characteristics
334 function = f;
335 jF = f.jacobian();
336 targetValues = target.clone();
337 residualsWeights = weights.clone();
338 this.point = startPoint.clone();
339 this.residuals = new double[target.length];
340
341 // arrays shared with the other private methods
342 rows = target.length;
343 cols = point.length;
344 jacobian = new double[rows][cols];
345
346 cost = Double.POSITIVE_INFINITY;
347
348 return doOptimize();
349
350 }
351
352 /** Perform the bulk of optimization algorithm.
353 * @return the point/value pair giving the optimal value for objective function
354 * @exception FunctionEvaluationException if the objective function throws one during
355 * the search
356 * @exception OptimizationException if the algorithm failed to converge
357 * @exception IllegalArgumentException if the start point dimension is wrong
358 */
359 protected abstract VectorialPointValuePair doOptimize()
360 throws FunctionEvaluationException, OptimizationException, IllegalArgumentException;
361
362 }