    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;
    019
    020    import org.apache.commons.math.FunctionEvaluationException;
    021    import org.apache.commons.math.MathRuntimeException;
    022    import org.apache.commons.math.analysis.MultivariateRealFunction;
    023    import org.apache.commons.math.analysis.MultivariateVectorialFunction;
    024    import org.apache.commons.math.exception.util.LocalizedFormats;
    025    import org.apache.commons.math.linear.RealMatrix;
    026
    027    /** This class converts {@link MultivariateVectorialFunction vectorial
    028     * objective functions} to {@link MultivariateRealFunction scalar objective functions}
    029     * when the goal is to minimize them.
    030     * <p>
    031     * This class is mostly used when the vectorial objective function represents
    032     * a theoretical result computed from a point set applied to a model and
    033     * the models point must be adjusted to fit the theoretical result to some
    034     * reference observations. The observations may be obtained for example from
    035     * physical measurements whether the model is built from theoretical
    036     * considerations.
    037     * </p>
    038     * <p>
    039     * This class computes a possibly weighted squared sum of the residuals, which is
    040     * a scalar value. The residuals are the difference between the theoretical model
    041     * (i.e. the output of the vectorial objective function) and the observations. The
    042     * class implements the {@link MultivariateRealFunction} interface and can therefore be
    043     * minimized by any optimizer supporting scalar objectives functions.This is one way
    044     * to perform a least square estimation. There are other ways to do this without using
    045     * this converter, as some optimization algorithms directly support vectorial objective
    046     * functions.
    047     * </p>
    048     * <p>
    049     * This class support combination of residuals with or without weights and correlations.
    050     * </p>
    051      *
    052     * @see MultivariateRealFunction
    053     * @see MultivariateVectorialFunction
    054     * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $
    055     * @since 2.0
    056     */
    057
    058    public class LeastSquaresConverter implements MultivariateRealFunction {
    059
    060        /** Underlying vectorial function. */
    061        private final MultivariateVectorialFunction function;
    062
    063        /** Observations to be compared to objective function to compute residuals. */
    064        private final double[] observations;
    065
    066        /** Optional weights for the residuals. */
    067        private final double[] weights;
    068
    069        /** Optional scaling matrix (weight and correlations) for the residuals. */
    070        private final RealMatrix scale;
    071
    072        /** Build a simple converter for uncorrelated residuals with the same weight.
    073         * @param function vectorial residuals function to wrap
    074         * @param observations observations to be compared to objective function to compute residuals
    075         */
    076        public LeastSquaresConverter(final MultivariateVectorialFunction function,
    077                                     final double[] observations) {
    078            this.function     = function;
    079            this.observations = observations.clone();
    080            this.weights      = null;
    081            this.scale        = null;
    082        }
    083
    084        /** Build a simple converter for uncorrelated residuals with the specific weights.
    085         * <p>
    086         * The scalar objective function value is computed as:
    087         * <pre>
    088         * objective = &sum;weight<sub>i</sub>(observation<sub>i</sub>-objective<sub>i</sub>)<sup>2</sup>
    089         * </pre>
    090         * </p>
    091         * <p>
    092         * Weights can be used for example to combine residuals with different standard
    093         * deviations. As an example, consider a residuals array in which even elements
    094         * are angular measurements in degrees with a 0.01&deg; standard deviation and
    095         * odd elements are distance measurements in meters with a 15m standard deviation.
    096         * In this case, the weights array should be initialized with value
    097         * 1.0/(0.01<sup>2</sup>) in the even elements and 1.0/(15.0<sup>2</sup>) in the
    098         * odd elements (i.e. reciprocals of variances).
    099         * </p>
    100         * <p>
    101         * The array computed by the objective function, the observations array and the
    102         * weights array must have consistent sizes or a {@link FunctionEvaluationException} will be
    103         * triggered while computing the scalar objective.
    104         * </p>
    105         * @param function vectorial residuals function to wrap
    106         * @param observations observations to be compared to objective function to compute residuals
    107         * @param weights weights to apply to the residuals
    108         * @exception IllegalArgumentException if the observations vector and the weights
    109         * vector dimensions don't match (objective function dimension is checked only when
    110         * the {@link #value(double[])} method is called)
    111         */
    112        public LeastSquaresConverter(final MultivariateVectorialFunction function,
    113                                     final double[] observations, final double[] weights)
    114            throws IllegalArgumentException {
    115            if (observations.length != weights.length) {
    116                throw MathRuntimeException.createIllegalArgumentException(
    117                        LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
    118                        observations.length, weights.length);
    119            }
    120            this.function     = function;
    121            this.observations = observations.clone();
    122            this.weights      = weights.clone();
    123            this.scale        = null;
    124        }
    125
    126        /** Build a simple converter for correlated residuals with the specific weights.
    127         * <p>
    128         * The scalar objective function value is computed as:
    129         * <pre>
    130         * objective = y<sup>T</sup>y with y = scale&times;(observation-objective)
    131         * </pre>
    132         * </p>
    133         * <p>
    134         * The array computed by the objective function, the observations array and the
    135         * the scaling matrix must have consistent sizes or a {@link FunctionEvaluationException}
    136         * will be triggered while computing the scalar objective.
    137         * </p>
    138         * @param function vectorial residuals function to wrap
    139         * @param observations observations to be compared to objective function to compute residuals
    140         * @param scale scaling matrix
    141         * @exception IllegalArgumentException if the observations vector and the scale
    142         * matrix dimensions don't match (objective function dimension is checked only when
    143         * the {@link #value(double[])} method is called)
    144         */
    145        public LeastSquaresConverter(final MultivariateVectorialFunction function,
    146                                     final double[] observations, final RealMatrix scale)
    147            throws IllegalArgumentException {
    148            if (observations.length != scale.getColumnDimension()) {
    149                throw MathRuntimeException.createIllegalArgumentException(
    150                        LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
    151                        observations.length, scale.getColumnDimension());
    152            }
    153            this.function     = function;
    154            this.observations = observations.clone();
    155            this.weights      = null;
    156            this.scale        = scale.copy();
    157        }
    158
    159        /** {@inheritDoc} */
    160        public double value(final double[] point) throws FunctionEvaluationException {
    161
    162            // compute residuals
    163            final double[] residuals = function.value(point);
    164            if (residuals.length != observations.length) {
    165                throw new FunctionEvaluationException(point,LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
    166                                            residuals.length, observations.length);
    167            }
    168            for (int i = 0; i < residuals.length; ++i) {
    169                residuals[i] -= observations[i];
    170            }
    171
    172            // compute sum of squares
    173            double sumSquares = 0;
    174            if (weights != null) {
    175                for (int i = 0; i < residuals.length; ++i) {
    176                    final double ri = residuals[i];
    177                    sumSquares +=  weights[i] * ri * ri;
    178                }
    179            } else if (scale != null) {
    180                for (final double yi : scale.operate(residuals)) {
    181                    sumSquares += yi * yi;
    182                }
    183            } else {
    184                for (final double ri : residuals) {
    185                    sumSquares += ri * ri;
    186                }
    187            }
    188
    189            return sumSquares;
    190
    191        }
    192
    193    }