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.math3.optim.nonlinear.scalar;
19  
20  import org.apache.commons.math3.analysis.MultivariateFunction;
21  import org.apache.commons.math3.analysis.MultivariateVectorFunction;
22  import org.apache.commons.math3.exception.DimensionMismatchException;
23  import org.apache.commons.math3.linear.RealMatrix;
24  
25  /**
26   * This class converts
27   * {@link MultivariateVectorFunction vectorial objective functions} to
28   * {@link MultivariateFunction scalar objective functions}
29   * when the goal is to minimize them.
30   * <br/>
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   * <br/>
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 MultivariateFunction} 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   * <br/>
47   * This class support combination of residuals with or without weights and correlations.
48    *
49   * @see MultivariateFunction
50   * @see MultivariateVectorFunction
51   * @since 2.0
52   */
53  
54  public class LeastSquaresConverter implements MultivariateFunction {
55      /** Underlying vectorial function. */
56      private final MultivariateVectorFunction function;
57      /** Observations to be compared to objective function to compute residuals. */
58      private final double[] observations;
59      /** Optional weights for the residuals. */
60      private final double[] weights;
61      /** Optional scaling matrix (weight and correlations) for the residuals. */
62      private final RealMatrix scale;
63  
64      /**
65       * Builds a simple converter for uncorrelated residuals with identical
66       * weights.
67       *
68       * @param function vectorial residuals function to wrap
69       * @param observations observations to be compared to objective function to compute residuals
70       */
71      public LeastSquaresConverter(final MultivariateVectorFunction function,
72                                   final double[] observations) {
73          this.function     = function;
74          this.observations = observations.clone();
75          this.weights      = null;
76          this.scale        = null;
77      }
78  
79      /**
80       * Builds a simple converter for uncorrelated residuals with the
81       * specified weights.
82       * <p>
83       * The scalar objective function value is computed as:
84       * <pre>
85       * objective = &sum;weight<sub>i</sub>(observation<sub>i</sub>-objective<sub>i</sub>)<sup>2</sup>
86       * </pre>
87       * </p>
88       * <p>
89       * Weights can be used for example to combine residuals with different standard
90       * deviations. As an example, consider a residuals array in which even elements
91       * are angular measurements in degrees with a 0.01&deg; standard deviation and
92       * odd elements are distance measurements in meters with a 15m standard deviation.
93       * In this case, the weights array should be initialized with value
94       * 1.0/(0.01<sup>2</sup>) in the even elements and 1.0/(15.0<sup>2</sup>) in the
95       * odd elements (i.e. reciprocals of variances).
96       * </p>
97       * <p>
98       * The array computed by the objective function, the observations array and the
99       * weights array must have consistent sizes or a {@link DimensionMismatchException}
100      * will be triggered while computing the scalar objective.
101      * </p>
102      *
103      * @param function vectorial residuals function to wrap
104      * @param observations observations to be compared to objective function to compute residuals
105      * @param weights weights to apply to the residuals
106      * @throws DimensionMismatchException if the observations vector and the weights
107      * vector dimensions do not match (objective function dimension is checked only when
108      * the {@link #value(double[])} method is called)
109      */
110     public LeastSquaresConverter(final MultivariateVectorFunction function,
111                                  final double[] observations,
112                                  final double[] weights) {
113         if (observations.length != weights.length) {
114             throw new DimensionMismatchException(observations.length, weights.length);
115         }
116         this.function     = function;
117         this.observations = observations.clone();
118         this.weights      = weights.clone();
119         this.scale        = null;
120     }
121 
122     /**
123      * Builds a simple converter for correlated residuals with the
124      * specified weights.
125      * <p>
126      * The scalar objective function value is computed as:
127      * <pre>
128      * objective = y<sup>T</sup>y with y = scale&times;(observation-objective)
129      * </pre>
130      * </p>
131      * <p>
132      * The array computed by the objective function, the observations array and the
133      * the scaling matrix must have consistent sizes or a {@link DimensionMismatchException}
134      * will be triggered while computing the scalar objective.
135      * </p>
136      *
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      * @throws DimensionMismatchException if the observations vector and the scale
141      * matrix dimensions do not match (objective function dimension is checked only when
142      * the {@link #value(double[])} method is called)
143      */
144     public LeastSquaresConverter(final MultivariateVectorFunction function,
145                                  final double[] observations,
146                                  final RealMatrix scale) {
147         if (observations.length != scale.getColumnDimension()) {
148             throw new DimensionMismatchException(observations.length, scale.getColumnDimension());
149         }
150         this.function     = function;
151         this.observations = observations.clone();
152         this.weights      = null;
153         this.scale        = scale.copy();
154     }
155 
156     /** {@inheritDoc} */
157     public double value(final double[] point) {
158         // compute residuals
159         final double[] residuals = function.value(point);
160         if (residuals.length != observations.length) {
161             throw new DimensionMismatchException(residuals.length, observations.length);
162         }
163         for (int i = 0; i < residuals.length; ++i) {
164             residuals[i] -= observations[i];
165         }
166 
167         // compute sum of squares
168         double sumSquares = 0;
169         if (weights != null) {
170             for (int i = 0; i < residuals.length; ++i) {
171                 final double ri = residuals[i];
172                 sumSquares +=  weights[i] * ri * ri;
173             }
174         } else if (scale != null) {
175             for (final double yi : scale.operate(residuals)) {
176                 sumSquares += yi * yi;
177             }
178         } else {
179             for (final double ri : residuals) {
180                 sumSquares += ri * ri;
181             }
182         }
183 
184         return sumSquares;
185     }
186 }