View Javadoc

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