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: 925812 $ $Date: 2010-03-21 11:49:31 -0400 (Sun, 21 Mar 2010) $
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        /**
067         * Target value for the objective functions at optimum.
068         * @since 2.1
069         */
070        protected double[] targetValues;
071    
072        /**
073         * Weight for the least squares cost computation.
074         * @since 2.1
075         */
076        protected double[] residualsWeights;
077    
078        /** Current point. */
079        protected double[] point;
080    
081        /** Current objective function value. */
082        protected double[] objective;
083    
084        /** Current residuals. */
085        protected double[] residuals;
086    
087        /** Cost value (square root of the sum of the residuals). */
088        protected double cost;
089    
090        /** Maximal number of iterations allowed. */
091        private int maxIterations;
092    
093        /** Number of iterations already performed. */
094        private int iterations;
095    
096        /** Maximal number of evaluations allowed. */
097        private int maxEvaluations;
098    
099        /** Number of evaluations already performed. */
100        private int objectiveEvaluations;
101    
102        /** Number of jacobian evaluations. */
103        private int jacobianEvaluations;
104    
105        /** Objective function. */
106        private DifferentiableMultivariateVectorialFunction function;
107    
108        /** Objective function derivatives. */
109        private MultivariateMatrixFunction jF;
110    
111        /** Simple constructor with default settings.
112         * <p>The convergence check is set to a {@link SimpleVectorialValueChecker}
113         * and the maximal number of evaluation is set to its default value.</p>
114         */
115        protected AbstractLeastSquaresOptimizer() {
116            setConvergenceChecker(new SimpleVectorialValueChecker());
117            setMaxIterations(DEFAULT_MAX_ITERATIONS);
118            setMaxEvaluations(Integer.MAX_VALUE);
119        }
120    
121        /** {@inheritDoc} */
122        public void setMaxIterations(int maxIterations) {
123            this.maxIterations = maxIterations;
124        }
125    
126        /** {@inheritDoc} */
127        public int getMaxIterations() {
128            return maxIterations;
129        }
130    
131        /** {@inheritDoc} */
132        public int getIterations() {
133            return iterations;
134        }
135    
136        /** {@inheritDoc} */
137        public void setMaxEvaluations(int maxEvaluations) {
138            this.maxEvaluations = maxEvaluations;
139        }
140    
141        /** {@inheritDoc} */
142        public int getMaxEvaluations() {
143            return maxEvaluations;
144        }
145    
146        /** {@inheritDoc} */
147        public int getEvaluations() {
148            return objectiveEvaluations;
149        }
150    
151        /** {@inheritDoc} */
152        public int getJacobianEvaluations() {
153            return jacobianEvaluations;
154        }
155    
156        /** {@inheritDoc} */
157        public void setConvergenceChecker(VectorialConvergenceChecker convergenceChecker) {
158            this.checker = convergenceChecker;
159        }
160    
161        /** {@inheritDoc} */
162        public VectorialConvergenceChecker getConvergenceChecker() {
163            return checker;
164        }
165    
166        /** Increment the iterations counter by 1.
167         * @exception OptimizationException if the maximal number
168         * of iterations is exceeded
169         */
170        protected void incrementIterationsCounter()
171            throws OptimizationException {
172            if (++iterations > maxIterations) {
173                throw new OptimizationException(new MaxIterationsExceededException(maxIterations));
174            }
175        }
176    
177        /**
178         * Update the jacobian matrix.
179         * @exception FunctionEvaluationException if the function jacobian
180         * cannot be evaluated or its dimension doesn't match problem dimension
181         */
182        protected void updateJacobian() throws FunctionEvaluationException {
183            ++jacobianEvaluations;
184            jacobian = jF.value(point);
185            if (jacobian.length != rows) {
186                throw new FunctionEvaluationException(point, "dimension mismatch {0} != {1}",
187                                                      jacobian.length, rows);
188            }
189            for (int i = 0; i < rows; i++) {
190                final double[] ji = jacobian[i];
191                final double factor = -Math.sqrt(residualsWeights[i]);
192                for (int j = 0; j < cols; ++j) {
193                    ji[j] *= factor;
194                }
195            }
196        }
197    
198        /**
199         * Update the residuals array and cost function value.
200         * @exception FunctionEvaluationException if the function cannot be evaluated
201         * or its dimension doesn't match problem dimension or maximal number of
202         * of evaluations is exceeded
203         */
204        protected void updateResidualsAndCost()
205            throws FunctionEvaluationException {
206    
207            if (++objectiveEvaluations > maxEvaluations) {
208                throw new FunctionEvaluationException(new MaxEvaluationsExceededException(maxEvaluations),
209                                                      point);
210            }
211            objective = function.value(point);
212            if (objective.length != rows) {
213                throw new FunctionEvaluationException(point, "dimension mismatch {0} != {1}",
214                                                      objective.length, rows);
215            }
216            cost = 0;
217            int index = 0;
218            for (int i = 0; i < rows; i++) {
219                final double residual = targetValues[i] - objective[i];
220                residuals[i] = residual;
221                cost += residualsWeights[i] * residual * residual;
222                index += cols;
223            }
224            cost = Math.sqrt(cost);
225    
226        }
227    
228        /**
229         * Get the Root Mean Square value.
230         * Get the Root Mean Square value, i.e. the root of the arithmetic
231         * mean of the square of all weighted residuals. This is related to the
232         * criterion that is minimized by the optimizer as follows: if
233         * <em>c</em> if the criterion, and <em>n</em> is the number of
234         * measurements, then the RMS is <em>sqrt (c/n)</em>.
235         *
236         * @return RMS value
237         */
238        public double getRMS() {
239            double criterion = 0;
240            for (int i = 0; i < rows; ++i) {
241                final double residual = residuals[i];
242                criterion += residualsWeights[i] * residual * residual;
243            }
244            return Math.sqrt(criterion / rows);
245        }
246    
247        /**
248         * Get the Chi-Square value.
249         * @return chi-square value
250         */
251        public double getChiSquare() {
252            double chiSquare = 0;
253            for (int i = 0; i < rows; ++i) {
254                final double residual = residuals[i];
255                chiSquare += residual * residual / residualsWeights[i];
256            }
257            return chiSquare;
258        }
259    
260        /**
261         * Get the covariance matrix of optimized parameters.
262         * @return covariance matrix
263         * @exception FunctionEvaluationException if the function jacobian cannot
264         * be evaluated
265         * @exception OptimizationException if the covariance matrix
266         * cannot be computed (singular problem)
267         */
268        public double[][] getCovariances()
269            throws FunctionEvaluationException, OptimizationException {
270    
271            // set up the jacobian
272            updateJacobian();
273    
274            // compute transpose(J).J, avoiding building big intermediate matrices
275            double[][] jTj = new double[cols][cols];
276            for (int i = 0; i < cols; ++i) {
277                for (int j = i; j < cols; ++j) {
278                    double sum = 0;
279                    for (int k = 0; k < rows; ++k) {
280                        sum += jacobian[k][i] * jacobian[k][j];
281                    }
282                    jTj[i][j] = sum;
283                    jTj[j][i] = sum;
284                }
285            }
286    
287            try {
288                // compute the covariances matrix
289                RealMatrix inverse =
290                    new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
291                return inverse.getData();
292            } catch (InvalidMatrixException ime) {
293                throw new OptimizationException("unable to compute covariances: singular problem");
294            }
295    
296        }
297    
298        /**
299         * Guess the errors in optimized parameters.
300         * <p>Guessing is covariance-based, it only gives rough order of magnitude.</p>
301         * @return errors in optimized parameters
302         * @exception FunctionEvaluationException if the function jacobian cannot b evaluated
303         * @exception OptimizationException if the covariances matrix cannot be computed
304         * or the number of degrees of freedom is not positive (number of measurements
305         * lesser or equal to number of parameters)
306         */
307        public double[] guessParametersErrors()
308            throws FunctionEvaluationException, OptimizationException {
309            if (rows <= cols) {
310                throw new OptimizationException(
311                        "no degrees of freedom ({0} measurements, {1} parameters)",
312                        rows, cols);
313            }
314            double[] errors = new double[cols];
315            final double c = Math.sqrt(getChiSquare() / (rows - cols));
316            double[][] covar = getCovariances();
317            for (int i = 0; i < errors.length; ++i) {
318                errors[i] = Math.sqrt(covar[i][i]) * c;
319            }
320            return errors;
321        }
322    
323        /** {@inheritDoc} */
324        public VectorialPointValuePair optimize(final DifferentiableMultivariateVectorialFunction f,
325                                                final double[] target, final double[] weights,
326                                                final double[] startPoint)
327            throws FunctionEvaluationException, OptimizationException, IllegalArgumentException {
328    
329            if (target.length != weights.length) {
330                throw new OptimizationException("dimension mismatch {0} != {1}",
331                                                target.length, weights.length);
332            }
333    
334            // reset counters
335            iterations           = 0;
336            objectiveEvaluations = 0;
337            jacobianEvaluations  = 0;
338    
339            // store least squares problem characteristics
340            function         = f;
341            jF               = f.jacobian();
342            targetValues     = target.clone();
343            residualsWeights = weights.clone();
344            this.point       = startPoint.clone();
345            this.residuals   = new double[target.length];
346    
347            // arrays shared with the other private methods
348            rows      = target.length;
349            cols      = point.length;
350            jacobian  = new double[rows][cols];
351    
352            cost = Double.POSITIVE_INFINITY;
353    
354            return doOptimize();
355    
356        }
357    
358        /** Perform the bulk of optimization algorithm.
359         * @return the point/value pair giving the optimal value for objective function
360         * @exception FunctionEvaluationException if the objective function throws one during
361         * the search
362         * @exception OptimizationException if the algorithm failed to converge
363         * @exception IllegalArgumentException if the start point dimension is wrong
364         */
365        protected abstract VectorialPointValuePair doOptimize()
366            throws FunctionEvaluationException, OptimizationException, IllegalArgumentException;
367    
368    }