View Javadoc

1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.apache.commons.math.optimization.general;
19  
20  import org.apache.commons.math.FunctionEvaluationException;
21  import org.apache.commons.math.MaxEvaluationsExceededException;
22  import org.apache.commons.math.MaxIterationsExceededException;
23  import org.apache.commons.math.analysis.DifferentiableMultivariateVectorialFunction;
24  import org.apache.commons.math.analysis.MultivariateMatrixFunction;
25  import org.apache.commons.math.linear.InvalidMatrixException;
26  import org.apache.commons.math.linear.LUDecompositionImpl;
27  import org.apache.commons.math.linear.MatrixUtils;
28  import org.apache.commons.math.linear.RealMatrix;
29  import org.apache.commons.math.optimization.OptimizationException;
30  import org.apache.commons.math.optimization.SimpleVectorialValueChecker;
31  import org.apache.commons.math.optimization.VectorialConvergenceChecker;
32  import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer;
33  import org.apache.commons.math.optimization.VectorialPointValuePair;
34  
35  /**
36   * Base class for implementing least squares optimizers.
37   * <p>This base class handles the boilerplate methods associated to thresholds
38   * settings, jacobian and error estimation.</p>
39   * @version $Revision: 786466 $ $Date: 2009-06-19 08:03:14 -0400 (Fri, 19 Jun 2009) $
40   * @since 1.2
41   *
42   */
43  public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMultivariateVectorialOptimizer {
44  
45      /** Default maximal number of iterations allowed. */
46      public static final int DEFAULT_MAX_ITERATIONS = 100;
47  
48      /** Maximal number of iterations allowed. */
49      private int maxIterations;
50  
51      /** Number of iterations already performed. */
52      private int iterations;
53  
54      /** Maximal number of evaluations allowed. */
55      private int maxEvaluations;
56  
57      /** Number of evaluations already performed. */
58      private int objectiveEvaluations;
59  
60      /** Number of jacobian evaluations. */
61      private int jacobianEvaluations;
62  
63      /** Convergence checker. */
64      protected VectorialConvergenceChecker checker;
65  
66      /** 
67       * Jacobian matrix.
68       * <p>This matrix is in canonical form just after the calls to
69       * {@link #updateJacobian()}, but may be modified by the solver
70       * in the derived class (the {@link LevenbergMarquardtOptimizer
71       * Levenberg-Marquardt optimizer} does this).</p>
72       */
73      protected double[][] jacobian;
74  
75      /** Number of columns of the jacobian matrix. */
76      protected int cols;
77  
78      /** Number of rows of the jacobian matrix. */
79      protected int rows;
80  
81      /** Objective function. */
82      private DifferentiableMultivariateVectorialFunction f;
83  
84      /** Objective function derivatives. */
85      private MultivariateMatrixFunction jF;
86  
87      /** Target value for the objective functions at optimum. */
88      protected double[] target;
89  
90      /** Weight for the least squares cost computation. */
91      protected double[] weights;
92  
93      /** Current point. */
94      protected double[] point;
95  
96      /** Current objective function value. */
97      protected double[] objective;
98  
99      /** 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 }