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 = ∑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° 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×(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 }