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.optimization;
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    /** This class converts {@link MultivariateVectorFunction vectorial
026     * objective functions} to {@link MultivariateFunction scalar objective functions}
027     * when the goal is to minimize them.
028     * <p>
029     * This class is mostly used when the vectorial objective function represents
030     * a theoretical result computed from a point set applied to a model and
031     * the models point must be adjusted to fit the theoretical result to some
032     * reference observations. The observations may be obtained for example from
033     * physical measurements whether the model is built from theoretical
034     * considerations.
035     * </p>
036     * <p>
037     * This class computes a possibly weighted squared sum of the residuals, which is
038     * a scalar value. The residuals are the difference between the theoretical model
039     * (i.e. the output of the vectorial objective function) and the observations. The
040     * class implements the {@link MultivariateFunction} interface and can therefore be
041     * minimized by any optimizer supporting scalar objectives functions.This is one way
042     * to perform a least square estimation. There are other ways to do this without using
043     * this converter, as some optimization algorithms directly support vectorial objective
044     * functions.
045     * </p>
046     * <p>
047     * This class support combination of residuals with or without weights and correlations.
048     * </p>
049      *
050     * @see MultivariateFunction
051     * @see MultivariateVectorFunction
052     * @version $Id: LeastSquaresConverter.java 1422230 2012-12-15 12:11:13Z erans $
053     * @deprecated As of 3.1 (to be removed in 4.0).
054     * @since 2.0
055     */
056    
057    @Deprecated
058    public class LeastSquaresConverter implements MultivariateFunction {
059    
060        /** Underlying vectorial function. */
061        private final MultivariateVectorFunction 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 MultivariateVectorFunction 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 DimensionMismatchException}
103         * will be 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 DimensionMismatchException if the observations vector and the weights
109         * vector dimensions do not match (objective function dimension is checked only when
110         * the {@link #value(double[])} method is called)
111         */
112        public LeastSquaresConverter(final MultivariateVectorFunction function,
113                                     final double[] observations, 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        /** Build a simple converter for correlated residuals with the specific weights.
124         * <p>
125         * The scalar objective function value is computed as:
126         * <pre>
127         * objective = y<sup>T</sup>y with y = scale&times;(observation-objective)
128         * </pre>
129         * </p>
130         * <p>
131         * The array computed by the objective function, the observations array and the
132         * the scaling matrix must have consistent sizes or a {@link DimensionMismatchException}
133         * will be triggered while computing the scalar objective.
134         * </p>
135         * @param function vectorial residuals function to wrap
136         * @param observations observations to be compared to objective function to compute residuals
137         * @param scale scaling matrix
138         * @throws DimensionMismatchException if the observations vector and the scale
139         * matrix dimensions do not match (objective function dimension is checked only when
140         * the {@link #value(double[])} method is called)
141         */
142        public LeastSquaresConverter(final MultivariateVectorFunction function,
143                                     final double[] observations, final RealMatrix scale) {
144            if (observations.length != scale.getColumnDimension()) {
145                throw new DimensionMismatchException(observations.length, scale.getColumnDimension());
146            }
147            this.function     = function;
148            this.observations = observations.clone();
149            this.weights      = null;
150            this.scale        = scale.copy();
151        }
152    
153        /** {@inheritDoc} */
154        public double value(final double[] point) {
155            // compute residuals
156            final double[] residuals = function.value(point);
157            if (residuals.length != observations.length) {
158                throw new DimensionMismatchException(residuals.length, observations.length);
159            }
160            for (int i = 0; i < residuals.length; ++i) {
161                residuals[i] -= observations[i];
162            }
163    
164            // compute sum of squares
165            double sumSquares = 0;
166            if (weights != null) {
167                for (int i = 0; i < residuals.length; ++i) {
168                    final double ri = residuals[i];
169                    sumSquares +=  weights[i] * ri * ri;
170                }
171            } else if (scale != null) {
172                for (final double yi : scale.operate(residuals)) {
173                    sumSquares += yi * yi;
174                }
175            } else {
176                for (final double ri : residuals) {
177                    sumSquares += ri * ri;
178                }
179            }
180    
181            return sumSquares;
182        }
183    }