Home | History | Annotate | Download | only in general 
      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 
     18 package org.apache.commons.math.optimization.general;
     19 
     20 import org.apache.commons.math.FunctionEvaluationException;
     21 import org.apache.commons.math.MaxEvaluationsExceededException;
     22 import org.apache.commons.math.MaxIterationsExceededException;
     23 import org.apache.commons.math.analysis.DifferentiableMultivariateVectorialFunction;
     24 import org.apache.commons.math.analysis.MultivariateMatrixFunction;
     25 import org.apache.commons.math.exception.util.LocalizedFormats;
     26 import org.apache.commons.math.linear.InvalidMatrixException;
     27 import org.apache.commons.math.linear.LUDecompositionImpl;
     28 import org.apache.commons.math.linear.MatrixUtils;
     29 import org.apache.commons.math.linear.RealMatrix;
     30 import org.apache.commons.math.optimization.OptimizationException;
     31 import org.apache.commons.math.optimization.SimpleVectorialValueChecker;
     32 import org.apache.commons.math.optimization.VectorialConvergenceChecker;
     33 import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer;
     34 import org.apache.commons.math.optimization.VectorialPointValuePair;
     35 import 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 fvr. 2011) $
     42  * @since 1.2
     43  *
     44  */
     45 public 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