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
18package org.apache.commons.math.optimization;
19
20import org.apache.commons.math.FunctionEvaluationException;
21import org.apache.commons.math.MathRuntimeException;
22import org.apache.commons.math.analysis.MultivariateRealFunction;
23import org.apache.commons.math.analysis.MultivariateVectorialFunction;
24import org.apache.commons.math.exception.util.LocalizedFormats;
25import org.apache.commons.math.linear.RealMatrix;
26
27/** This class converts {@link MultivariateVectorialFunction vectorial
28 * objective functions} to {@link MultivariateRealFunction scalar objective functions}
29 * when the goal is to minimize them.
30 * <p>
31 * This class is mostly used when the vectorial objective function represents
32 * a theoretical result computed from a point set applied to a model and
33 * the models point must be adjusted to fit the theoretical result to some
34 * reference observations. The observations may be obtained for example from
35 * physical measurements whether the model is built from theoretical
36 * considerations.
37 * </p>
38 * <p>
39 * This class computes a possibly weighted squared sum of the residuals, which is
40 * a scalar value. The residuals are the difference between the theoretical model
41 * (i.e. the output of the vectorial objective function) and the observations. The
42 * class implements the {@link MultivariateRealFunction} interface and can therefore be
43 * minimized by any optimizer supporting scalar objectives functions.This is one way
44 * to perform a least square estimation. There are other ways to do this without using
45 * this converter, as some optimization algorithms directly support vectorial objective
46 * functions.
47 * </p>
48 * <p>
49 * This class support combination of residuals with or without weights and correlations.
50 * </p>
51  *
52 * @see MultivariateRealFunction
53 * @see MultivariateVectorialFunction
54 * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 févr. 2011) $
55 * @since 2.0
56 */
57
58public class LeastSquaresConverter implements MultivariateRealFunction {
59
60    /** Underlying vectorial function. */
61    private final MultivariateVectorialFunction function;
62
63    /** Observations to be compared to objective function to compute residuals. */
64    private final double[] observations;
65
66    /** Optional weights for the residuals. */
67    private final double[] weights;
68
69    /** Optional scaling matrix (weight and correlations) for the residuals. */
70    private final RealMatrix scale;
71
72    /** Build a simple converter for uncorrelated residuals with the same weight.
73     * @param function vectorial residuals function to wrap
74     * @param observations observations to be compared to objective function to compute residuals
75     */
76    public LeastSquaresConverter(final MultivariateVectorialFunction function,
77                                 final double[] observations) {
78        this.function     = function;
79        this.observations = observations.clone();
80        this.weights      = null;
81        this.scale        = null;
82    }
83
84    /** Build a simple converter for uncorrelated residuals with the specific weights.
85     * <p>
86     * The scalar objective function value is computed as:
87     * <pre>
88     * objective = &sum;weight<sub>i</sub>(observation<sub>i</sub>-objective<sub>i</sub>)<sup>2</sup>
89     * </pre>
90     * </p>
91     * <p>
92     * Weights can be used for example to combine residuals with different standard
93     * deviations. As an example, consider a residuals array in which even elements
94     * are angular measurements in degrees with a 0.01&deg; standard deviation and
95     * odd elements are distance measurements in meters with a 15m standard deviation.
96     * In this case, the weights array should be initialized with value
97     * 1.0/(0.01<sup>2</sup>) in the even elements and 1.0/(15.0<sup>2</sup>) in the
98     * odd elements (i.e. reciprocals of variances).
99     * </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}
194