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
18package org.apache.commons.math.optimization.general;
19
20import org.apache.commons.math.FunctionEvaluationException;
21import org.apache.commons.math.MaxEvaluationsExceededException;
22import org.apache.commons.math.MaxIterationsExceededException;
23import org.apache.commons.math.analysis.DifferentiableMultivariateVectorialFunction;
24import org.apache.commons.math.analysis.MultivariateMatrixFunction;
25import org.apache.commons.math.exception.util.LocalizedFormats;
26import org.apache.commons.math.linear.InvalidMatrixException;
27import org.apache.commons.math.linear.LUDecompositionImpl;
28import org.apache.commons.math.linear.MatrixUtils;
29import org.apache.commons.math.linear.RealMatrix;
30import org.apache.commons.math.optimization.OptimizationException;
31import org.apache.commons.math.optimization.SimpleVectorialValueChecker;
32import org.apache.commons.math.optimization.VectorialConvergenceChecker;
33import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer;
34import org.apache.commons.math.optimization.VectorialPointValuePair;
35import org.apache.commons.math.util.FastMath;
36
37/**
38 * Base class for implementing least squares optimizers.
39 * <p>This base class handles the boilerplate methods associated to thresholds
40 * settings, jacobian and error estimation.</p>
41 * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $
42 * @since 1.2
43 *
44 */
45public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMultivariateVectorialOptimizer {
46
47    /** Default maximal number of iterations allowed. */
48    public static final int DEFAULT_MAX_ITERATIONS = 100;
49
50    /** Convergence checker. */
51    protected VectorialConvergenceChecker checker;
52
53    /**
54     * Jacobian matrix.
55     * <p>This matrix is in canonical form just after the calls to
56     * {@link #updateJacobian()}, but may be modified by the solver
57     * in the derived class (the {@link LevenbergMarquardtOptimizer
58     * Levenberg-Marquardt optimizer} does this).</p>
59     */
60    protected double[][] jacobian;
61
62    /** Number of columns of the jacobian matrix. */
63    protected int cols;
64
65    /** Number of rows of the jacobian matrix. */
66    protected int rows;
67
68    /**
69     * Target value for the objective functions at optimum.
70     * @since 2.1
71     */
72    protected double[] targetValues;
73
74    /**
75     * Weight for the least squares cost computation.
76     * @since 2.1
77     */
78    protected double[] residualsWeights;
79
80    /** Current point. */
81    protected double[] point;
82
83    /** Current objective function value. */
84    protected double[] objective;
85
86    /** Current residuals. */
87    protected double[] residuals;
88
89    /** Weighted Jacobian */
90    protected double[][] wjacobian;
91
92    /** Weighted residuals */
93    protected double[] wresiduals;
94
95    /** Cost value (square root of the sum of the residuals). */
96    protected double cost;
97
98    /** Maximal number of iterations allowed. */
99    private int maxIterations;
100
101    /** Number of iterations already performed. */
102    private int iterations;
103
104    /** Maximal number of evaluations allowed. */
105    private int maxEvaluations;
106
107    /** Number of evaluations already performed. */
108    private int objectiveEvaluations;
109
110    /** Number of jacobian evaluations. */
111    private int jacobianEvaluations;
112
113    /** Objective function. */
114    private DifferentiableMultivariateVectorialFunction function;
115
116    /** Objective function derivatives. */
117    private MultivariateMatrixFunction jF;
118
119    /** Simple constructor with default settings.
120     * <p>The convergence check is set to a {@link SimpleVectorialValueChecker}
121     * and the maximal number of evaluation is set to its default value.</p>
122     */
123    protected AbstractLeastSquaresOptimizer() {
124        setConvergenceChecker(new SimpleVectorialValueChecker());
125        setMaxIterations(DEFAULT_MAX_ITERATIONS);
126        setMaxEvaluations(Integer.MAX_VALUE);
127    }
128
129    /** {@inheritDoc} */
130    public void setMaxIterations(int maxIterations) {
131        this.maxIterations = maxIterations;
132    }
133
134    /** {@inheritDoc} */
135    public int getMaxIterations() {
136        return maxIterations;
137    }
138
139    /** {@inheritDoc} */
140    public int getIterations() {
141        return iterations;
142    }
143
144    /** {@inheritDoc} */
145    public void setMaxEvaluations(int maxEvaluations) {
146        this.maxEvaluations = maxEvaluations;
147    }
148
149    /** {@inheritDoc} */
150    public int getMaxEvaluations() {
151        return maxEvaluations;
152    }
153
154    /** {@inheritDoc} */
155    public int getEvaluations() {
156        return objectiveEvaluations;
157    }
158
159    /** {@inheritDoc} */
160    public int getJacobianEvaluations() {
161        return jacobianEvaluations;
162    }
163
164    /** {@inheritDoc} */
165    public void setConvergenceChecker(VectorialConvergenceChecker convergenceChecker) {
166        this.checker = convergenceChecker;
167    }
168
169    /** {@inheritDoc} */
170    public VectorialConvergenceChecker getConvergenceChecker() {
171        return checker;
172    }
173
174    /** Increment the iterations counter by 1.
175     * @exception OptimizationException if the maximal number
176     * of iterations is exceeded
177     */
178    protected void incrementIterationsCounter()
179        throws OptimizationException {
180        if (++iterations > maxIterations) {
181            throw new OptimizationException(new MaxIterationsExceededException(maxIterations));
182        }
183    }
184
185    /**
186     * Update the jacobian matrix.
187     * @exception FunctionEvaluationException if the function jacobian
188     * cannot be evaluated or its dimension doesn't match problem dimension
189     */
190    protected void updateJacobian() throws FunctionEvaluationException {
191        ++jacobianEvaluations;
192        jacobian = jF.value(point);
193        if (jacobian.length != rows) {
194            throw new FunctionEvaluationException(point, LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
195                                                  jacobian.length, rows);
196        }
197        for (int i = 0; i < rows; i++) {
198            final double[] ji = jacobian[i];
199            double wi = FastMath.sqrt(residualsWeights[i]);
200            for (int j = 0; j < cols; ++j) {
201                ji[j] *=  -1.0;
202                wjacobian[i][j] = ji[j]*wi;
203            }
204        }
205    }
206
207    /**
208     * Update the residuals array and cost function value.
209     * @exception FunctionEvaluationException if the function cannot be evaluated
210     * or its dimension doesn't match problem dimension or maximal number of
211     * of evaluations is exceeded
212     */
213    protected void updateResidualsAndCost()
214        throws FunctionEvaluationException {
215
216        if (++objectiveEvaluations > maxEvaluations) {
217            throw new FunctionEvaluationException(new MaxEvaluationsExceededException(maxEvaluations),
218                                                  point);
219        }
220        objective = function.value(point);
221        if (objective.length != rows) {
222            throw new FunctionEvaluationException(point, LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
223                                                  objective.length, rows);
224        }
225        cost = 0;
226        int index = 0;
227        for (int i = 0; i < rows; i++) {
228            final double residual = targetValues[i] - objective[i];
229            residuals[i] = residual;
230            wresiduals[i]= residual*FastMath.sqrt(residualsWeights[i]);
231            cost += residualsWeights[i] * residual * residual;
232            index += cols;
233        }
234        cost = FastMath.sqrt(cost);
235
236    }
237
238    /**
239     * Get the Root Mean Square value.
240     * Get the Root Mean Square value, i.e. the root of the arithmetic
241     * mean of the square of all weighted residuals. This is related to the
242     * criterion that is minimized by the optimizer as follows: if
243     * <em>c</em> if the criterion, and <em>n</em> is the number of
244     * measurements, then the RMS is <em>sqrt (c/n)</em>.
245     *
246     * @return RMS value
247     */
248    public double getRMS() {
249        return FastMath.sqrt(getChiSquare() / rows);
250    }
251
252    /**
253     * Get a Chi-Square-like value assuming the N residuals follow N
254     * distinct normal distributions centered on 0 and whose variances are
255     * the reciprocal of the weights.
256     * @return chi-square value
257     */
258    public double getChiSquare() {
259        return cost*cost;
260    }
261
262    /**
263     * Get the covariance matrix of optimized parameters.
264     * @return covariance matrix
265     * @exception FunctionEvaluationException if the function jacobian cannot
266     * be evaluated
267     * @exception OptimizationException if the covariance matrix
268     * cannot be computed (singular problem)
269     */
270    public double[][] getCovariances()
271        throws FunctionEvaluationException, OptimizationException {
272
273        // set up the jacobian
274        updateJacobian();
275
276        // compute transpose(J).J, avoiding building big intermediate matrices
277        double[][] jTj = new double[cols][cols];
278        for (int i = 0; i < cols; ++i) {
279            for (int j = i; j < cols; ++j) {
280                double sum = 0;
281                for (int k = 0; k < rows; ++k) {
282                    sum += wjacobian[k][i] * wjacobian[k][j];
283                }
284                jTj[i][j] = sum;
285                jTj[j][i] = sum;
286            }
287        }
288
289        try {
290            // compute the covariance matrix
291            RealMatrix inverse =
292                new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
293            return inverse.getData();
294        } catch (InvalidMatrixException ime) {
295            throw new OptimizationException(LocalizedFormats.UNABLE_TO_COMPUTE_COVARIANCE_SINGULAR_PROBLEM);
296        }
297
298    }
299
300    /**
301     * Guess the errors in optimized parameters.
302     * <p>Guessing is covariance-based, it only gives rough order of magnitude.</p>
303     * @return errors in optimized parameters
304     * @exception FunctionEvaluationException if the function jacobian cannot b evaluated
305     * @exception OptimizationException if the covariances matrix cannot be computed
306     * or the number of degrees of freedom is not positive (number of measurements
307     * lesser or equal to number of parameters)
308     */
309    public double[] guessParametersErrors()
310        throws FunctionEvaluationException, OptimizationException {
311        if (rows <= cols) {
312            throw new OptimizationException(
313                    LocalizedFormats.NO_DEGREES_OF_FREEDOM,
314                    rows, cols);
315        }
316        double[] errors = new double[cols];
317        final double c = FastMath.sqrt(getChiSquare() / (rows - cols));
318        double[][] covar = getCovariances();
319        for (int i = 0; i < errors.length; ++i) {
320            errors[i] = FastMath.sqrt(covar[i][i]) * c;
321        }
322        return errors;
323    }
324
325    /** {@inheritDoc} */
326    public VectorialPointValuePair optimize(final DifferentiableMultivariateVectorialFunction f,
327                                            final double[] target, final double[] weights,
328                                            final double[] startPoint)
329        throws FunctionEvaluationException, OptimizationException, IllegalArgumentException {
330
331        if (target.length != weights.length) {
332            throw new OptimizationException(LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
333                                            target.length, weights.length);
334        }
335
336        // reset counters
337        iterations           = 0;
338        objectiveEvaluations = 0;
339        jacobianEvaluations  = 0;
340
341        // store least squares problem characteristics
342        function         = f;
343        jF               = f.jacobian();
344        targetValues     = target.clone();
345        residualsWeights = weights.clone();
346        this.point       = startPoint.clone();
347        this.residuals   = new double[target.length];
348
349        // arrays shared with the other private methods
350        rows      = target.length;
351        cols      = point.length;
352        jacobian  = new double[rows][cols];
353
354        wjacobian = new double[rows][cols];
355        wresiduals = new double[rows];
356
357        cost = Double.POSITIVE_INFINITY;
358
359        return doOptimize();
360
361    }
362
363    /** Perform the bulk of optimization algorithm.
364     * @return the point/value pair giving the optimal value for objective function
365     * @exception FunctionEvaluationException if the objective function throws one during
366     * the search
367     * @exception OptimizationException if the algorithm failed to converge
368     * @exception IllegalArgumentException if the start point dimension is wrong
369     */
370    protected abstract VectorialPointValuePair doOptimize()
371        throws FunctionEvaluationException, OptimizationException, IllegalArgumentException;
372
373
374}
375