001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018 package org.apache.commons.math3.optimization.general;
019
020 import org.apache.commons.math3.analysis.DifferentiableMultivariateVectorFunction;
021 import org.apache.commons.math3.analysis.FunctionUtils;
022 import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
023 import org.apache.commons.math3.analysis.differentiation.MultivariateDifferentiableVectorFunction;
024 import org.apache.commons.math3.exception.DimensionMismatchException;
025 import org.apache.commons.math3.exception.NumberIsTooSmallException;
026 import org.apache.commons.math3.exception.util.LocalizedFormats;
027 import org.apache.commons.math3.linear.ArrayRealVector;
028 import org.apache.commons.math3.linear.RealMatrix;
029 import org.apache.commons.math3.linear.DiagonalMatrix;
030 import org.apache.commons.math3.linear.DecompositionSolver;
031 import org.apache.commons.math3.linear.MatrixUtils;
032 import org.apache.commons.math3.linear.QRDecomposition;
033 import org.apache.commons.math3.linear.EigenDecomposition;
034 import org.apache.commons.math3.optimization.OptimizationData;
035 import org.apache.commons.math3.optimization.InitialGuess;
036 import org.apache.commons.math3.optimization.Target;
037 import org.apache.commons.math3.optimization.Weight;
038 import org.apache.commons.math3.optimization.ConvergenceChecker;
039 import org.apache.commons.math3.optimization.DifferentiableMultivariateVectorOptimizer;
040 import org.apache.commons.math3.optimization.PointVectorValuePair;
041 import org.apache.commons.math3.optimization.direct.BaseAbstractMultivariateVectorOptimizer;
042 import org.apache.commons.math3.util.FastMath;
043
044 /**
045 * Base class for implementing least squares optimizers.
046 * It handles the boilerplate methods associated to thresholds settings,
047 * Jacobian and error estimation.
048 * <br/>
049 * This class constructs the Jacobian matrix of the function argument in method
050 * {@link BaseAbstractMultivariateVectorOptimizer#optimize(int,MultivariateVectorFunction,OptimizationData[])
051 * optimize} and assumes that the rows of that matrix iterate on the model
052 * functions while the columns iterate on the parameters; thus, the numbers
053 * of rows is equal to the dimension of the
054 * {@link org.apache.commons.math3.optimization.Target Target} while
055 * the number of columns is equal to the dimension of the
056 * {@link org.apache.commons.math3.optimization.InitialGuess InitialGuess}.
057 *
058 * @version $Id: AbstractLeastSquaresOptimizer.java 1426759 2012-12-29 13:26:44Z erans $
059 * @deprecated As of 3.1 (to be removed in 4.0).
060 * @since 1.2
061 */
062 @Deprecated
063 public abstract class AbstractLeastSquaresOptimizer
064 extends BaseAbstractMultivariateVectorOptimizer<DifferentiableMultivariateVectorFunction>
065 implements DifferentiableMultivariateVectorOptimizer {
066 /**
067 * Singularity threshold (cf. {@link #getCovariances(double)}).
068 * @deprecated As of 3.1.
069 */
070 @Deprecated
071 private static final double DEFAULT_SINGULARITY_THRESHOLD = 1e-14;
072 /**
073 * Jacobian matrix of the weighted residuals.
074 * This matrix is in canonical form just after the calls to
075 * {@link #updateJacobian()}, but may be modified by the solver
076 * in the derived class (the {@link LevenbergMarquardtOptimizer
077 * Levenberg-Marquardt optimizer} does this).
078 * @deprecated As of 3.1. To be removed in 4.0. Please use
079 * {@link #computeWeightedJacobian(double[])} instead.
080 */
081 @Deprecated
082 protected double[][] weightedResidualJacobian;
083 /** Number of columns of the jacobian matrix.
084 * @deprecated As of 3.1.
085 */
086 @Deprecated
087 protected int cols;
088 /** Number of rows of the jacobian matrix.
089 * @deprecated As of 3.1.
090 */
091 @Deprecated
092 protected int rows;
093 /** Current point.
094 * @deprecated As of 3.1.
095 */
096 @Deprecated
097 protected double[] point;
098 /** Current objective function value.
099 * @deprecated As of 3.1.
100 */
101 @Deprecated
102 protected double[] objective;
103 /** Weighted residuals
104 * @deprecated As of 3.1.
105 */
106 @Deprecated
107 protected double[] weightedResiduals;
108 /** Cost value (square root of the sum of the residuals).
109 * @deprecated As of 3.1. Field to become "private" in 4.0.
110 * Please use {@link #setCost(double)}.
111 */
112 @Deprecated
113 protected double cost;
114 /** Objective function derivatives. */
115 private MultivariateDifferentiableVectorFunction jF;
116 /** Number of evaluations of the Jacobian. */
117 private int jacobianEvaluations;
118 /** Square-root of the weight matrix. */
119 private RealMatrix weightMatrixSqrt;
120
121 /**
122 * Simple constructor with default settings.
123 * The convergence check is set to a {@link
124 * org.apache.commons.math3.optimization.SimpleVectorValueChecker}.
125 * @deprecated See {@link org.apache.commons.math3.optimization.SimpleValueChecker#SimpleValueChecker()}
126 */
127 @Deprecated
128 protected AbstractLeastSquaresOptimizer() {}
129
130 /**
131 * @param checker Convergence checker.
132 */
133 protected AbstractLeastSquaresOptimizer(ConvergenceChecker<PointVectorValuePair> checker) {
134 super(checker);
135 }
136
137 /**
138 * @return the number of evaluations of the Jacobian function.
139 */
140 public int getJacobianEvaluations() {
141 return jacobianEvaluations;
142 }
143
144 /**
145 * Update the jacobian matrix.
146 *
147 * @throws DimensionMismatchException if the Jacobian dimension does not
148 * match problem dimension.
149 * @deprecated As of 3.1. Please use {@link #computeWeightedJacobian(double[])}
150 * instead.
151 */
152 @Deprecated
153 protected void updateJacobian() {
154 final RealMatrix weightedJacobian = computeWeightedJacobian(point);
155 weightedResidualJacobian = weightedJacobian.scalarMultiply(-1).getData();
156 }
157
158 /**
159 * Computes the Jacobian matrix.
160 *
161 * @param params Model parameters at which to compute the Jacobian.
162 * @return the weighted Jacobian: W<sup>1/2</sup> J.
163 * @throws DimensionMismatchException if the Jacobian dimension does not
164 * match problem dimension.
165 * @since 3.1
166 */
167 protected RealMatrix computeWeightedJacobian(double[] params) {
168 ++jacobianEvaluations;
169
170 final DerivativeStructure[] dsPoint = new DerivativeStructure[params.length];
171 final int nC = params.length;
172 for (int i = 0; i < nC; ++i) {
173 dsPoint[i] = new DerivativeStructure(nC, 1, i, params[i]);
174 }
175 final DerivativeStructure[] dsValue = jF.value(dsPoint);
176 final int nR = getTarget().length;
177 if (dsValue.length != nR) {
178 throw new DimensionMismatchException(dsValue.length, nR);
179 }
180 final double[][] jacobianData = new double[nR][nC];
181 for (int i = 0; i < nR; ++i) {
182 int[] orders = new int[nC];
183 for (int j = 0; j < nC; ++j) {
184 orders[j] = 1;
185 jacobianData[i][j] = dsValue[i].getPartialDerivative(orders);
186 orders[j] = 0;
187 }
188 }
189
190 return weightMatrixSqrt.multiply(MatrixUtils.createRealMatrix(jacobianData));
191 }
192
193 /**
194 * Update the residuals array and cost function value.
195 * @throws DimensionMismatchException if the dimension does not match the
196 * problem dimension.
197 * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
198 * if the maximal number of evaluations is exceeded.
199 * @deprecated As of 3.1. Please use {@link #computeResiduals(double[])},
200 * {@link #computeObjectiveValue(double[])}, {@link #computeCost(double[])}
201 * and {@link #setCost(double)} instead.
202 */
203 @Deprecated
204 protected void updateResidualsAndCost() {
205 objective = computeObjectiveValue(point);
206 final double[] res = computeResiduals(objective);
207
208 // Compute cost.
209 cost = computeCost(res);
210
211 // Compute weighted residuals.
212 final ArrayRealVector residuals = new ArrayRealVector(res);
213 weightedResiduals = weightMatrixSqrt.operate(residuals).toArray();
214 }
215
216 /**
217 * Computes the cost.
218 *
219 * @param residuals Residuals.
220 * @return the cost.
221 * @see #computeResiduals(double[])
222 * @since 3.1
223 */
224 protected double computeCost(double[] residuals) {
225 final ArrayRealVector r = new ArrayRealVector(residuals);
226 return FastMath.sqrt(r.dotProduct(getWeight().operate(r)));
227 }
228
229 /**
230 * Get the Root Mean Square value.
231 * Get the Root Mean Square value, i.e. the root of the arithmetic
232 * mean of the square of all weighted residuals. This is related to the
233 * criterion that is minimized by the optimizer as follows: if
234 * <em>c</em> if the criterion, and <em>n</em> is the number of
235 * measurements, then the RMS is <em>sqrt (c/n)</em>.
236 *
237 * @return RMS value
238 */
239 public double getRMS() {
240 return FastMath.sqrt(getChiSquare() / rows);
241 }
242
243 /**
244 * Get a Chi-Square-like value assuming the N residuals follow N
245 * distinct normal distributions centered on 0 and whose variances are
246 * the reciprocal of the weights.
247 * @return chi-square value
248 */
249 public double getChiSquare() {
250 return cost * cost;
251 }
252
253 /**
254 * Gets the square-root of the weight matrix.
255 *
256 * @return the square-root of the weight matrix.
257 * @since 3.1
258 */
259 public RealMatrix getWeightSquareRoot() {
260 return weightMatrixSqrt.copy();
261 }
262
263 /**
264 * Sets the cost.
265 *
266 * @param cost Cost value.
267 * @since 3.1
268 */
269 protected void setCost(double cost) {
270 this.cost = cost;
271 }
272
273 /**
274 * Get the covariance matrix of the optimized parameters.
275 *
276 * @return the covariance matrix.
277 * @throws org.apache.commons.math3.linear.SingularMatrixException
278 * if the covariance matrix cannot be computed (singular problem).
279 * @see #getCovariances(double)
280 * @deprecated As of 3.1. Please use {@link #computeCovariances(double[],double)}
281 * instead.
282 */
283 @Deprecated
284 public double[][] getCovariances() {
285 return getCovariances(DEFAULT_SINGULARITY_THRESHOLD);
286 }
287
288 /**
289 * Get the covariance matrix of the optimized parameters.
290 * <br/>
291 * Note that this operation involves the inversion of the
292 * <code>J<sup>T</sup>J</code> matrix, where {@code J} is the
293 * Jacobian matrix.
294 * The {@code threshold} parameter is a way for the caller to specify
295 * that the result of this computation should be considered meaningless,
296 * and thus trigger an exception.
297 *
298 * @param threshold Singularity threshold.
299 * @return the covariance matrix.
300 * @throws org.apache.commons.math3.linear.SingularMatrixException
301 * if the covariance matrix cannot be computed (singular problem).
302 * @deprecated As of 3.1. Please use {@link #computeCovariances(double[],double)}
303 * instead.
304 */
305 @Deprecated
306 public double[][] getCovariances(double threshold) {
307 return computeCovariances(point, threshold);
308 }
309
310 /**
311 * Get the covariance matrix of the optimized parameters.
312 * <br/>
313 * Note that this operation involves the inversion of the
314 * <code>J<sup>T</sup>J</code> matrix, where {@code J} is the
315 * Jacobian matrix.
316 * The {@code threshold} parameter is a way for the caller to specify
317 * that the result of this computation should be considered meaningless,
318 * and thus trigger an exception.
319 *
320 * @param params Model parameters.
321 * @param threshold Singularity threshold.
322 * @return the covariance matrix.
323 * @throws org.apache.commons.math3.linear.SingularMatrixException
324 * if the covariance matrix cannot be computed (singular problem).
325 * @since 3.1
326 */
327 public double[][] computeCovariances(double[] params,
328 double threshold) {
329 // Set up the Jacobian.
330 final RealMatrix j = computeWeightedJacobian(params);
331
332 // Compute transpose(J)J.
333 final RealMatrix jTj = j.transpose().multiply(j);
334
335 // Compute the covariances matrix.
336 final DecompositionSolver solver
337 = new QRDecomposition(jTj, threshold).getSolver();
338 return solver.getInverse().getData();
339 }
340
341 /**
342 * <p>
343 * Returns an estimate of the standard deviation of each parameter. The
344 * returned values are the so-called (asymptotic) standard errors on the
345 * parameters, defined as {@code sd(a[i]) = sqrt(S / (n - m) * C[i][i])},
346 * where {@code a[i]} is the optimized value of the {@code i}-th parameter,
347 * {@code S} is the minimized value of the sum of squares objective function
348 * (as returned by {@link #getChiSquare()}), {@code n} is the number of
349 * observations, {@code m} is the number of parameters and {@code C} is the
350 * covariance matrix.
351 * </p>
352 * <p>
353 * See also
354 * <a href="http://en.wikipedia.org/wiki/Least_squares">Wikipedia</a>,
355 * or
356 * <a href="http://mathworld.wolfram.com/LeastSquaresFitting.html">MathWorld</a>,
357 * equations (34) and (35) for a particular case.
358 * </p>
359 *
360 * @return an estimate of the standard deviation of the optimized parameters
361 * @throws org.apache.commons.math3.linear.SingularMatrixException
362 * if the covariance matrix cannot be computed.
363 * @throws NumberIsTooSmallException if the number of degrees of freedom is not
364 * positive, i.e. the number of measurements is less or equal to the number of
365 * parameters.
366 * @deprecated as of version 3.1, {@link #computeSigma(double[],double)} should be used
367 * instead. It should be emphasized that {@code guessParametersErrors} and
368 * {@code computeSigma} are <em>not</em> strictly equivalent.
369 */
370 @Deprecated
371 public double[] guessParametersErrors() {
372 if (rows <= cols) {
373 throw new NumberIsTooSmallException(LocalizedFormats.NO_DEGREES_OF_FREEDOM,
374 rows, cols, false);
375 }
376 double[] errors = new double[cols];
377 final double c = FastMath.sqrt(getChiSquare() / (rows - cols));
378 double[][] covar = computeCovariances(point, 1e-14);
379 for (int i = 0; i < errors.length; ++i) {
380 errors[i] = FastMath.sqrt(covar[i][i]) * c;
381 }
382 return errors;
383 }
384
385 /**
386 * Computes an estimate of the standard deviation of the parameters. The
387 * returned values are the square root of the diagonal coefficients of the
388 * covariance matrix, {@code sd(a[i]) ~= sqrt(C[i][i])}, where {@code a[i]}
389 * is the optimized value of the {@code i}-th parameter, and {@code C} is
390 * the covariance matrix.
391 *
392 * @param params Model parameters.
393 * @param covarianceSingularityThreshold Singularity threshold (see
394 * {@link #computeCovariances(double[],double) computeCovariances}).
395 * @return an estimate of the standard deviation of the optimized parameters
396 * @throws org.apache.commons.math3.linear.SingularMatrixException
397 * if the covariance matrix cannot be computed.
398 * @since 3.1
399 */
400 public double[] computeSigma(double[] params,
401 double covarianceSingularityThreshold) {
402 final int nC = params.length;
403 final double[] sig = new double[nC];
404 final double[][] cov = computeCovariances(params, covarianceSingularityThreshold);
405 for (int i = 0; i < nC; ++i) {
406 sig[i] = FastMath.sqrt(cov[i][i]);
407 }
408 return sig;
409 }
410
411 /** {@inheritDoc}
412 * @deprecated As of 3.1. Please use
413 * {@link BaseAbstractMultivariateVectorOptimizer#optimize(int,MultivariateVectorFunction,OptimizationData[])
414 * optimize(int,MultivariateDifferentiableVectorFunction,OptimizationData...)}
415 * instead.
416 */
417 @Override
418 @Deprecated
419 public PointVectorValuePair optimize(int maxEval,
420 final DifferentiableMultivariateVectorFunction f,
421 final double[] target, final double[] weights,
422 final double[] startPoint) {
423 return optimizeInternal(maxEval,
424 FunctionUtils.toMultivariateDifferentiableVectorFunction(f),
425 new Target(target),
426 new Weight(weights),
427 new InitialGuess(startPoint));
428 }
429
430 /**
431 * Optimize an objective function.
432 * Optimization is considered to be a weighted least-squares minimization.
433 * The cost function to be minimized is
434 * <code>∑weight<sub>i</sub>(objective<sub>i</sub> - target<sub>i</sub>)<sup>2</sup></code>
435 *
436 * @param f Objective function.
437 * @param target Target value for the objective functions at optimum.
438 * @param weights Weights for the least squares cost computation.
439 * @param startPoint Start point for optimization.
440 * @return the point/value pair giving the optimal value for objective
441 * function.
442 * @param maxEval Maximum number of function evaluations.
443 * @throws org.apache.commons.math3.exception.DimensionMismatchException
444 * if the start point dimension is wrong.
445 * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
446 * if the maximal number of evaluations is exceeded.
447 * @throws org.apache.commons.math3.exception.NullArgumentException if
448 * any argument is {@code null}.
449 * @deprecated As of 3.1. Please use
450 * {@link BaseAbstractMultivariateVectorOptimizer#optimize(int,MultivariateVectorFunction,OptimizationData[])
451 * optimize(int,MultivariateDifferentiableVectorFunction,OptimizationData...)}
452 * instead.
453 */
454 @Deprecated
455 public PointVectorValuePair optimize(final int maxEval,
456 final MultivariateDifferentiableVectorFunction f,
457 final double[] target, final double[] weights,
458 final double[] startPoint) {
459 return optimizeInternal(maxEval, f,
460 new Target(target),
461 new Weight(weights),
462 new InitialGuess(startPoint));
463 }
464
465 /**
466 * Optimize an objective function.
467 * Optimization is considered to be a weighted least-squares minimization.
468 * The cost function to be minimized is
469 * <code>∑weight<sub>i</sub>(objective<sub>i</sub> - target<sub>i</sub>)<sup>2</sup></code>
470 *
471 * @param maxEval Allowed number of evaluations of the objective function.
472 * @param f Objective function.
473 * @param optData Optimization data. The following data will be looked for:
474 * <ul>
475 * <li>{@link Target}</li>
476 * <li>{@link Weight}</li>
477 * <li>{@link InitialGuess}</li>
478 * </ul>
479 * @return the point/value pair giving the optimal value of the objective
480 * function.
481 * @throws org.apache.commons.math3.exception.TooManyEvaluationsException if
482 * the maximal number of evaluations is exceeded.
483 * @throws DimensionMismatchException if the target, and weight arguments
484 * have inconsistent dimensions.
485 * @see BaseAbstractMultivariateVectorOptimizer#optimizeInternal(int,MultivariateVectorFunction,OptimizationData[])
486 * @since 3.1
487 * @deprecated As of 3.1. Override is necessary only until this class's generic
488 * argument is changed to {@code MultivariateDifferentiableVectorFunction}.
489 */
490 @Deprecated
491 protected PointVectorValuePair optimizeInternal(final int maxEval,
492 final MultivariateDifferentiableVectorFunction f,
493 OptimizationData... optData) {
494 // XXX Conversion will be removed when the generic argument of the
495 // base class becomes "MultivariateDifferentiableVectorFunction".
496 return super.optimizeInternal(maxEval, FunctionUtils.toDifferentiableMultivariateVectorFunction(f), optData);
497 }
498
499 /** {@inheritDoc} */
500 @Override
501 protected void setUp() {
502 super.setUp();
503
504 // Reset counter.
505 jacobianEvaluations = 0;
506
507 // Square-root of the weight matrix.
508 weightMatrixSqrt = squareRoot(getWeight());
509
510 // Store least squares problem characteristics.
511 // XXX The conversion won't be necessary when the generic argument of
512 // the base class becomes "MultivariateDifferentiableVectorFunction".
513 // XXX "jF" is not strictly necessary anymore but is currently more
514 // efficient than converting the value returned from "getObjectiveFunction()"
515 // every time it is used.
516 jF = FunctionUtils.toMultivariateDifferentiableVectorFunction((DifferentiableMultivariateVectorFunction) getObjectiveFunction());
517
518 // Arrays shared with "private" and "protected" methods.
519 point = getStartPoint();
520 rows = getTarget().length;
521 cols = point.length;
522 }
523
524 /**
525 * Computes the residuals.
526 * The residual is the difference between the observed (target)
527 * values and the model (objective function) value.
528 * There is one residual for each element of the vector-valued
529 * function.
530 *
531 * @param objectiveValue Value of the the objective function. This is
532 * the value returned from a call to
533 * {@link #computeObjectiveValue(double[]) computeObjectiveValue}
534 * (whose array argument contains the model parameters).
535 * @return the residuals.
536 * @throws DimensionMismatchException if {@code params} has a wrong
537 * length.
538 * @since 3.1
539 */
540 protected double[] computeResiduals(double[] objectiveValue) {
541 final double[] target = getTarget();
542 if (objectiveValue.length != target.length) {
543 throw new DimensionMismatchException(target.length,
544 objectiveValue.length);
545 }
546
547 final double[] residuals = new double[target.length];
548 for (int i = 0; i < target.length; i++) {
549 residuals[i] = target[i] - objectiveValue[i];
550 }
551
552 return residuals;
553 }
554
555 /**
556 * Computes the square-root of the weight matrix.
557 *
558 * @param m Symmetric, positive-definite (weight) matrix.
559 * @return the square-root of the weight matrix.
560 */
561 private RealMatrix squareRoot(RealMatrix m) {
562 if (m instanceof DiagonalMatrix) {
563 final int dim = m.getRowDimension();
564 final RealMatrix sqrtM = new DiagonalMatrix(dim);
565 for (int i = 0; i < dim; i++) {
566 sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
567 }
568 return sqrtM;
569 } else {
570 final EigenDecomposition dec = new EigenDecomposition(m);
571 return dec.getSquareRoot();
572 }
573 }
574 }