1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math.optimization.general;
19
20 import org.apache.commons.math.FunctionEvaluationException;
21 import org.apache.commons.math.MaxEvaluationsExceededException;
22 import org.apache.commons.math.MaxIterationsExceededException;
23 import org.apache.commons.math.analysis.DifferentiableMultivariateVectorialFunction;
24 import org.apache.commons.math.analysis.MultivariateMatrixFunction;
25 import org.apache.commons.math.linear.InvalidMatrixException;
26 import org.apache.commons.math.linear.LUDecompositionImpl;
27 import org.apache.commons.math.linear.MatrixUtils;
28 import org.apache.commons.math.linear.RealMatrix;
29 import org.apache.commons.math.optimization.OptimizationException;
30 import org.apache.commons.math.optimization.SimpleVectorialValueChecker;
31 import org.apache.commons.math.optimization.VectorialConvergenceChecker;
32 import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer;
33 import org.apache.commons.math.optimization.VectorialPointValuePair;
34
35
36
37
38
39
40
41
42
43 public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMultivariateVectorialOptimizer {
44
45
46 public static final int DEFAULT_MAX_ITERATIONS = 100;
47
48
49 private int maxIterations;
50
51
52 private int iterations;
53
54
55 private int maxEvaluations;
56
57
58 private int objectiveEvaluations;
59
60
61 private int jacobianEvaluations;
62
63
64 protected VectorialConvergenceChecker checker;
65
66
67
68
69
70
71
72
73 protected double[][] jacobian;
74
75
76 protected int cols;
77
78
79 protected int rows;
80
81
82 private DifferentiableMultivariateVectorialFunction f;
83
84
85 private MultivariateMatrixFunction jF;
86
87
88 protected double[] target;
89
90
91 protected double[] weights;
92
93
94 protected double[] point;
95
96
97 protected double[] objective;
98
99
100 protected double[] residuals;
101
102
103 protected double cost;
104
105
106
107
108
109 protected AbstractLeastSquaresOptimizer() {
110 setConvergenceChecker(new SimpleVectorialValueChecker());
111 setMaxIterations(DEFAULT_MAX_ITERATIONS);
112 setMaxEvaluations(Integer.MAX_VALUE);
113 }
114
115
116 public void setMaxIterations(int maxIterations) {
117 this.maxIterations = maxIterations;
118 }
119
120
121 public int getMaxIterations() {
122 return maxIterations;
123 }
124
125
126 public int getIterations() {
127 return iterations;
128 }
129
130
131 public void setMaxEvaluations(int maxEvaluations) {
132 this.maxEvaluations = maxEvaluations;
133 }
134
135
136 public int getMaxEvaluations() {
137 return maxEvaluations;
138 }
139
140
141 public int getEvaluations() {
142 return objectiveEvaluations;
143 }
144
145
146 public int getJacobianEvaluations() {
147 return jacobianEvaluations;
148 }
149
150
151 public void setConvergenceChecker(VectorialConvergenceChecker checker) {
152 this.checker = checker;
153 }
154
155
156 public VectorialConvergenceChecker getConvergenceChecker() {
157 return checker;
158 }
159
160
161
162
163
164 protected void incrementIterationsCounter()
165 throws OptimizationException {
166 if (++iterations > maxIterations) {
167 throw new OptimizationException(new MaxIterationsExceededException(maxIterations));
168 }
169 }
170
171
172
173
174
175
176 protected void updateJacobian() throws FunctionEvaluationException {
177 ++jacobianEvaluations;
178 jacobian = jF.value(point);
179 if (jacobian.length != rows) {
180 throw new FunctionEvaluationException(point, "dimension mismatch {0} != {1}",
181 jacobian.length, rows);
182 }
183 for (int i = 0; i < rows; i++) {
184 final double[] ji = jacobian[i];
185 final double factor = -Math.sqrt(weights[i]);
186 for (int j = 0; j < cols; ++j) {
187 ji[j] *= factor;
188 }
189 }
190 }
191
192
193
194
195
196
197
198 protected void updateResidualsAndCost()
199 throws FunctionEvaluationException {
200
201 if (++objectiveEvaluations > maxEvaluations) {
202 throw new FunctionEvaluationException(new MaxEvaluationsExceededException(maxEvaluations),
203 point);
204 }
205 objective = f.value(point);
206 if (objective.length != rows) {
207 throw new FunctionEvaluationException(point, "dimension mismatch {0} != {1}",
208 objective.length, rows);
209 }
210 cost = 0;
211 for (int i = 0, index = 0; i < rows; i++, index += cols) {
212 final double residual = target[i] - objective[i];
213 residuals[i] = residual;
214 cost += weights[i] * residual * residual;
215 }
216 cost = Math.sqrt(cost);
217
218 }
219
220
221
222
223
224
225
226
227
228
229
230 public double getRMS() {
231 double criterion = 0;
232 for (int i = 0; i < rows; ++i) {
233 final double residual = residuals[i];
234 criterion += weights[i] * residual * residual;
235 }
236 return Math.sqrt(criterion / rows);
237 }
238
239
240
241
242
243 public double getChiSquare() {
244 double chiSquare = 0;
245 for (int i = 0; i < rows; ++i) {
246 final double residual = residuals[i];
247 chiSquare += residual * residual / weights[i];
248 }
249 return chiSquare;
250 }
251
252
253
254
255
256
257
258
259
260 public double[][] getCovariances()
261 throws FunctionEvaluationException, OptimizationException {
262
263
264 updateJacobian();
265
266
267 double[][] jTj = new double[cols][cols];
268 for (int i = 0; i < cols; ++i) {
269 for (int j = i; j < cols; ++j) {
270 double sum = 0;
271 for (int k = 0; k < rows; ++k) {
272 sum += jacobian[k][i] * jacobian[k][j];
273 }
274 jTj[i][j] = sum;
275 jTj[j][i] = sum;
276 }
277 }
278
279 try {
280
281 RealMatrix inverse =
282 new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
283 return inverse.getData();
284 } catch (InvalidMatrixException ime) {
285 throw new OptimizationException("unable to compute covariances: singular problem");
286 }
287
288 }
289
290
291
292
293
294
295
296
297
298
299 public double[] guessParametersErrors()
300 throws FunctionEvaluationException, OptimizationException {
301 if (rows <= cols) {
302 throw new OptimizationException(
303 "no degrees of freedom ({0} measurements, {1} parameters)",
304 rows, cols);
305 }
306 double[] errors = new double[cols];
307 final double c = Math.sqrt(getChiSquare() / (rows - cols));
308 double[][] covar = getCovariances();
309 for (int i = 0; i < errors.length; ++i) {
310 errors[i] = Math.sqrt(covar[i][i]) * c;
311 }
312 return errors;
313 }
314
315
316 public VectorialPointValuePair optimize(final DifferentiableMultivariateVectorialFunction f,
317 final double[] target, final double[] weights,
318 final double[] startPoint)
319 throws FunctionEvaluationException, OptimizationException, IllegalArgumentException {
320
321 if (target.length != weights.length) {
322 throw new OptimizationException("dimension mismatch {0} != {1}",
323 target.length, weights.length);
324 }
325
326
327 iterations = 0;
328 objectiveEvaluations = 0;
329 jacobianEvaluations = 0;
330
331
332 this.f = f;
333 jF = f.jacobian();
334 this.target = target.clone();
335 this.weights = weights.clone();
336 this.point = startPoint.clone();
337 this.residuals = new double[target.length];
338
339
340 rows = target.length;
341 cols = point.length;
342 jacobian = new double[rows][cols];
343
344 cost = Double.POSITIVE_INFINITY;
345
346 return doOptimize();
347
348 }
349
350
351
352
353
354
355
356
357 abstract protected VectorialPointValuePair doOptimize()
358 throws FunctionEvaluationException, OptimizationException, IllegalArgumentException;
359
360 }