    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    }