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    *      https://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  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  
23  package org.hipparchus.optim.nonlinear.scalar;
24  
25  import org.hipparchus.analysis.MultivariateFunction;
26  import org.hipparchus.analysis.MultivariateVectorFunction;
27  import org.hipparchus.exception.LocalizedCoreFormats;
28  import org.hipparchus.exception.MathIllegalArgumentException;
29  import org.hipparchus.linear.RealMatrix;
30  
31  /**
32   * This class converts
33   * {@link MultivariateVectorFunction vectorial objective functions} to
34   * {@link MultivariateFunction scalar objective functions}
35   * when the goal is to minimize them.
36   * <br>
37   * This class is mostly used when the vectorial objective function represents
38   * a theoretical result computed from a point set applied to a model and
39   * the models point must be adjusted to fit the theoretical result to some
40   * reference observations. The observations may be obtained for example from
41   * physical measurements whether the model is built from theoretical
42   * considerations.
43   * <br>
44   * This class computes a possibly weighted squared sum of the residuals, which is
45   * a scalar value. The residuals are the difference between the theoretical model
46   * (i.e. the output of the vectorial objective function) and the observations. The
47   * class implements the {@link MultivariateFunction} interface and can therefore be
48   * minimized by any optimizer supporting scalar objectives functions.This is one way
49   * to perform a least square estimation. There are other ways to do this without using
50   * this converter, as some optimization algorithms directly support vectorial objective
51   * functions.
52   * <br>
53   * This class support combination of residuals with or without weights and correlations.
54    *
55   * @see MultivariateFunction
56   * @see MultivariateVectorFunction
57   */
58  
59  public class LeastSquaresConverter implements MultivariateFunction {
60      /** Underlying vectorial function. */
61      private final MultivariateVectorFunction function;
62      /** Observations to be compared to objective function to compute residuals. */
63      private final double[] observations;
64      /** Optional weights for the residuals. */
65      private final double[] weights;
66      /** Optional scaling matrix (weight and correlations) for the residuals. */
67      private final RealMatrix scale;
68  
69      /**
70       * Builds a simple converter for uncorrelated residuals with identical
71       * weights.
72       *
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 MultivariateVectorFunction 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      /**
85       * Builds a simple converter for uncorrelated residuals with the
86       * specified weights.
87       * <p>
88       * The scalar objective function value is computed as:
89       * objective = &sum;weight<sub>i</sub>(observation<sub>i</sub>-objective<sub>i</sub>)<sup>2</sup>
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 MathIllegalArgumentException}
103      * will be triggered while computing the scalar objective.
104      * </p>
105      *
106      * @param function vectorial residuals function to wrap
107      * @param observations observations to be compared to objective function to compute residuals
108      * @param weights weights to apply to the residuals
109      * @throws MathIllegalArgumentException if the observations vector and the weights
110      * vector dimensions do not match (objective function dimension is checked only when
111      * the {@link #value(double[])} method is called)
112      */
113     public LeastSquaresConverter(final MultivariateVectorFunction function,
114                                  final double[] observations,
115                                  final double[] weights) {
116         if (observations.length != weights.length) {
117             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
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     /**
127      * Builds a simple converter for correlated residuals with the
128      * specified weights.
129      * <p>
130      * The scalar objective function value is computed as:
131      * objective = y<sup>T</sup>y with y = scale&times;(observation-objective)
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 MathIllegalArgumentException}
136      * will be triggered while computing the scalar objective.
137      * </p>
138      *
139      * @param function vectorial residuals function to wrap
140      * @param observations observations to be compared to objective function to compute residuals
141      * @param scale scaling matrix
142      * @throws MathIllegalArgumentException if the observations vector and the scale
143      * matrix dimensions do not match (objective function dimension is checked only when
144      * the {@link #value(double[])} method is called)
145      */
146     public LeastSquaresConverter(final MultivariateVectorFunction function,
147                                  final double[] observations,
148                                  final RealMatrix scale) {
149         if (observations.length != scale.getColumnDimension()) {
150             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
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     @Override
161     public double value(final double[] point) {
162         // compute residuals
163         final double[] residuals = function.value(point);
164         if (residuals.length != observations.length) {
165             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
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 }