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.linear.RealMatrix;
025    
026    /** This class converts {@link MultivariateVectorialFunction vectorial
027     * objective functions} to {@link MultivariateRealFunction scalar objective functions}
028     * when the goal is to minimize them.
029     * <p>
030     * This class is mostly used when the vectorial objective function represents
031     * a theoretical result computed from a point set applied to a model and
032     * the models point must be adjusted to fit the theoretical result to some
033     * reference observations. The observations may be obtained for example from
034     * physical measurements whether the model is built from theoretical
035     * considerations.
036     * </p>
037     * <p>
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 MultivariateRealFunction} 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     * </p>
047     * <p>
048     * This class support combination of residuals with or without weights and correlations.
049     * </p>
050      *
051     * @see MultivariateRealFunction
052     * @see MultivariateVectorialFunction
053     * @version $Revision: 780674 $ $Date: 2009-06-01 11:10:55 -0400 (Mon, 01 Jun 2009) $
054     * @since 2.0
055     */
056    
057    public class LeastSquaresConverter implements MultivariateRealFunction {
058    
059        /** Underlying vectorial function. */
060        private final MultivariateVectorialFunction function;
061    
062        /** Observations to be compared to objective function to compute residuals. */
063        private final double[] observations;
064    
065        /** Optional weights for the residuals. */
066        private final double[] weights;
067    
068        /** Optional scaling matrix (weight and correlations) for the residuals. */
069        private final RealMatrix scale;
070    
071        /** Build a simple converter for uncorrelated residuals with the same weight.
072         * @param function vectorial residuals function to wrap
073         * @param observations observations to be compared to objective function to compute residuals
074         */
075        public LeastSquaresConverter(final MultivariateVectorialFunction function,
076                                     final double[] observations) {
077            this.function     = function;
078            this.observations = observations.clone();
079            this.weights      = null;
080            this.scale        = null;
081        }
082    
083        /** Build a simple converter for uncorrelated residuals with the specific weights.
084         * <p>
085         * The scalar objective function value is computed as:
086         * <pre>
087         * objective = &sum;weight<sub>i</sub>(observation<sub>i</sub>-objective<sub>i</sub>)<sup>2</sup>
088         * </pre>
089         * </p>
090         * <p>
091         * Weights can be used for example to combine residuals with different standard
092         * deviations. As an example, consider a residuals array in which even elements
093         * are angular measurements in degrees with a 0.01&deg; standard deviation and
094         * odd elements are distance measurements in meters with a 15m standard deviation.
095         * In this case, the weights array should be initialized with value
096         * 1.0/(0.01<sup>2</sup>) in the even elements and 1.0/(15.0<sup>2</sup>) in the
097         * odd elements (i.e. reciprocals of variances). 
098         * </p>
099         * <p>
100         * The array computed by the objective function, the observations array and the
101         * weights array must have consistent sizes or a {@link FunctionEvaluationException} will be
102         * triggered while computing the scalar objective.
103         * </p>
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         * @exception IllegalArgumentException if the observations vector and the weights
108         * vector dimensions don't match (objective function dimension is checked only when
109         * the {@link #value(double[])} method is called)
110         */
111        public LeastSquaresConverter(final MultivariateVectorialFunction function,
112                                     final double[] observations, final double[] weights)
113            throws IllegalArgumentException {
114            if (observations.length != weights.length) {
115                throw MathRuntimeException.createIllegalArgumentException(
116                        "dimension mismatch {0} != {1}",
117                        observations.length, weights.length);
118            }
119            this.function     = function;
120            this.observations = observations.clone();
121            this.weights      = weights.clone();
122            this.scale        = null;
123        }
124    
125        /** Build a simple converter for correlated residuals with the specific 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 FunctionEvaluationException}
135         * will be triggered while computing the scalar objective.
136         * </p>
137         * @param function vectorial residuals function to wrap
138         * @param observations observations to be compared to objective function to compute residuals
139         * @param scale scaling matrix
140         * @exception IllegalArgumentException if the observations vector and the scale
141         * matrix dimensions don't match (objective function dimension is checked only when
142         * the {@link #value(double[])} method is called)
143         */
144        public LeastSquaresConverter(final MultivariateVectorialFunction function,
145                                     final double[] observations, final RealMatrix scale)
146            throws IllegalArgumentException {
147            if (observations.length != scale.getColumnDimension()) {
148                throw MathRuntimeException.createIllegalArgumentException(
149                        "dimension mismatch {0} != {1}",
150                        observations.length, scale.getColumnDimension());
151            }
152            this.function     = function;
153            this.observations = observations.clone();
154            this.weights      = null;
155            this.scale        = scale.copy();
156        }
157    
158        /** {@inheritDoc} */
159        public double value(final double[] point) throws FunctionEvaluationException {
160    
161            // compute residuals
162            final double[] residuals = function.value(point);
163            if (residuals.length != observations.length) {
164                throw new FunctionEvaluationException(point, "dimension mismatch {0} != {1}",
165                                                      residuals.length, observations.length);
166            }
167            for (int i = 0; i < residuals.length; ++i) {
168                residuals[i] -= observations[i];
169            }
170    
171            // compute sum of squares
172            double sumSquares = 0;
173            if (weights != null) {
174                for (int i = 0; i < residuals.length; ++i) {
175                    final double ri = residuals[i];
176                    sumSquares +=  weights[i] * ri * ri;
177                }
178            } else if (scale != null) {
179                for (final double yi : scale.operate(residuals)) {
180                    sumSquares += yi * yi;
181                }
182            } else {
183                for (final double ri : residuals) {
184                    sumSquares += ri * ri;
185                }
186            }
187    
188            return sumSquares;
189    
190        }
191    
192    }