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: 786466 $ $Date: 2009-06-19 08:03:14 -0400 (Fri, 19 Jun 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        /** Maximal number of iterations allowed. */
049        private int maxIterations;
050    
051        /** Number of iterations already performed. */
052        private int iterations;
053    
054        /** Maximal number of evaluations allowed. */
055        private int maxEvaluations;
056    
057        /** Number of evaluations already performed. */
058        private int objectiveEvaluations;
059    
060        /** Number of jacobian evaluations. */
061        private int jacobianEvaluations;
062    
063        /** Convergence checker. */
064        protected VectorialConvergenceChecker checker;
065    
066        /** 
067         * Jacobian matrix.
068         * <p>This matrix is in canonical form just after the calls to
069         * {@link #updateJacobian()}, but may be modified by the solver
070         * in the derived class (the {@link LevenbergMarquardtOptimizer
071         * Levenberg-Marquardt optimizer} does this).</p>
072         */
073        protected double[][] jacobian;
074    
075        /** Number of columns of the jacobian matrix. */
076        protected int cols;
077    
078        /** Number of rows of the jacobian matrix. */
079        protected int rows;
080    
081        /** Objective function. */
082        private DifferentiableMultivariateVectorialFunction f;
083    
084        /** Objective function derivatives. */
085        private MultivariateMatrixFunction jF;
086    
087        /** Target value for the objective functions at optimum. */
088        protected double[] target;
089    
090        /** Weight for the least squares cost computation. */
091        protected double[] weights;
092    
093        /** Current point. */
094        protected double[] point;
095    
096        /** Current objective function value. */
097        protected double[] objective;
098    
099        /** 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    }