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