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.exception.util.LocalizedFormats;
026    import org.apache.commons.math.linear.InvalidMatrixException;
027    import org.apache.commons.math.linear.LUDecompositionImpl;
028    import org.apache.commons.math.linear.MatrixUtils;
029    import org.apache.commons.math.linear.RealMatrix;
030    import org.apache.commons.math.optimization.OptimizationException;
031    import org.apache.commons.math.optimization.SimpleVectorialValueChecker;
032    import org.apache.commons.math.optimization.VectorialConvergenceChecker;
033    import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer;
034    import org.apache.commons.math.optimization.VectorialPointValuePair;
035    import org.apache.commons.math.util.FastMath;
036    
037    /**
038     * Base class for implementing least squares optimizers.
039     * <p>This base class handles the boilerplate methods associated to thresholds
040     * settings, jacobian and error estimation.</p>
041     * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 f??vr. 2011) $
042     * @since 1.2
043     *
044     */
045    public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMultivariateVectorialOptimizer {
046    
047        /** Default maximal number of iterations allowed. */
048        public static final int DEFAULT_MAX_ITERATIONS = 100;
049    
050        /** Convergence checker. */
051        protected VectorialConvergenceChecker checker;
052    
053        /**
054         * Jacobian matrix.
055         * <p>This matrix is in canonical form just after the calls to
056         * {@link #updateJacobian()}, but may be modified by the solver
057         * in the derived class (the {@link LevenbergMarquardtOptimizer
058         * Levenberg-Marquardt optimizer} does this).</p>
059         */
060        protected double[][] jacobian;
061    
062        /** Number of columns of the jacobian matrix. */
063        protected int cols;
064    
065        /** Number of rows of the jacobian matrix. */
066        protected int rows;
067    
068        /**
069         * Target value for the objective functions at optimum.
070         * @since 2.1
071         */
072        protected double[] targetValues;
073    
074        /**
075         * Weight for the least squares cost computation.
076         * @since 2.1
077         */
078        protected double[] residualsWeights;
079    
080        /** Current point. */
081        protected double[] point;
082    
083        /** Current objective function value. */
084        protected double[] objective;
085    
086        /** Current residuals. */
087        protected double[] residuals;
088    
089        /** Weighted Jacobian */
090        protected double[][] wjacobian;
091    
092        /** Weighted residuals */
093        protected double[] wresiduals;
094    
095        /** Cost value (square root of the sum of the residuals). */
096        protected double cost;
097    
098        /** Maximal number of iterations allowed. */
099        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    }