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 package org.apache.commons.math3.optim.nonlinear.vector.jacobian;
018
019 import org.apache.commons.math3.exception.DimensionMismatchException;
020 import org.apache.commons.math3.exception.TooManyEvaluationsException;
021 import org.apache.commons.math3.linear.ArrayRealVector;
022 import org.apache.commons.math3.linear.RealMatrix;
023 import org.apache.commons.math3.linear.DiagonalMatrix;
024 import org.apache.commons.math3.linear.DecompositionSolver;
025 import org.apache.commons.math3.linear.MatrixUtils;
026 import org.apache.commons.math3.linear.QRDecomposition;
027 import org.apache.commons.math3.linear.EigenDecomposition;
028 import org.apache.commons.math3.optim.OptimizationData;
029 import org.apache.commons.math3.optim.ConvergenceChecker;
030 import org.apache.commons.math3.optim.PointVectorValuePair;
031 import org.apache.commons.math3.optim.nonlinear.vector.Weight;
032 import org.apache.commons.math3.optim.nonlinear.vector.JacobianMultivariateVectorOptimizer;
033 import org.apache.commons.math3.util.FastMath;
034
035 /**
036 * Base class for implementing least-squares optimizers.
037 * It provides methods for error estimation.
038 *
039 * @version $Id$
040 * @since 3.1
041 */
042 public abstract class AbstractLeastSquaresOptimizer
043 extends JacobianMultivariateVectorOptimizer {
044 /** Square-root of the weight matrix. */
045 private RealMatrix weightMatrixSqrt;
046 /** Cost value (square root of the sum of the residuals). */
047 private double cost;
048
049 /**
050 * @param checker Convergence checker.
051 */
052 protected AbstractLeastSquaresOptimizer(ConvergenceChecker<PointVectorValuePair> checker) {
053 super(checker);
054 }
055
056 /**
057 * Computes the weighted Jacobian matrix.
058 *
059 * @param params Model parameters at which to compute the Jacobian.
060 * @return the weighted Jacobian: W<sup>1/2</sup> J.
061 * @throws DimensionMismatchException if the Jacobian dimension does not
062 * match problem dimension.
063 */
064 protected RealMatrix computeWeightedJacobian(double[] params) {
065 return weightMatrixSqrt.multiply(MatrixUtils.createRealMatrix(computeJacobian(params)));
066 }
067
068 /**
069 * Computes the cost.
070 *
071 * @param residuals Residuals.
072 * @return the cost.
073 * @see #computeResiduals(double[])
074 */
075 protected double computeCost(double[] residuals) {
076 final ArrayRealVector r = new ArrayRealVector(residuals);
077 return FastMath.sqrt(r.dotProduct(getWeight().operate(r)));
078 }
079
080 /**
081 * Gets the root-mean-square (RMS) value.
082 *
083 * The RMS the root of the arithmetic mean of the square of all weighted
084 * residuals.
085 * This is related to the criterion that is minimized by the optimizer
086 * as follows: If <em>c</em> if the criterion, and <em>n</em> is the
087 * number of measurements, then the RMS is <em>sqrt (c/n)</em>.
088 *
089 * @return the RMS value.
090 */
091 public double getRMS() {
092 return FastMath.sqrt(getChiSquare() / getTargetSize());
093 }
094
095 /**
096 * Get a Chi-Square-like value assuming the N residuals follow N
097 * distinct normal distributions centered on 0 and whose variances are
098 * the reciprocal of the weights.
099 * @return chi-square value
100 */
101 public double getChiSquare() {
102 return cost * cost;
103 }
104
105 /**
106 * Gets the square-root of the weight matrix.
107 *
108 * @return the square-root of the weight matrix.
109 */
110 public RealMatrix getWeightSquareRoot() {
111 return weightMatrixSqrt.copy();
112 }
113
114 /**
115 * Sets the cost.
116 *
117 * @param cost Cost value.
118 */
119 protected void setCost(double cost) {
120 this.cost = cost;
121 }
122
123 /**
124 * Get the covariance matrix of the optimized parameters.
125 * <br/>
126 * Note that this operation involves the inversion of the
127 * <code>J<sup>T</sup>J</code> matrix, where {@code J} is the
128 * Jacobian matrix.
129 * The {@code threshold} parameter is a way for the caller to specify
130 * that the result of this computation should be considered meaningless,
131 * and thus trigger an exception.
132 *
133 * @param params Model parameters.
134 * @param threshold Singularity threshold.
135 * @return the covariance matrix.
136 * @throws org.apache.commons.math3.linear.SingularMatrixException
137 * if the covariance matrix cannot be computed (singular problem).
138 */
139 public double[][] computeCovariances(double[] params,
140 double threshold) {
141 // Set up the Jacobian.
142 final RealMatrix j = computeWeightedJacobian(params);
143
144 // Compute transpose(J)J.
145 final RealMatrix jTj = j.transpose().multiply(j);
146
147 // Compute the covariances matrix.
148 final DecompositionSolver solver
149 = new QRDecomposition(jTj, threshold).getSolver();
150 return solver.getInverse().getData();
151 }
152
153 /**
154 * Computes an estimate of the standard deviation of the parameters. The
155 * returned values are the square root of the diagonal coefficients of the
156 * covariance matrix, {@code sd(a[i]) ~= sqrt(C[i][i])}, where {@code a[i]}
157 * is the optimized value of the {@code i}-th parameter, and {@code C} is
158 * the covariance matrix.
159 *
160 * @param params Model parameters.
161 * @param covarianceSingularityThreshold Singularity threshold (see
162 * {@link #computeCovariances(double[],double) computeCovariances}).
163 * @return an estimate of the standard deviation of the optimized parameters
164 * @throws org.apache.commons.math3.linear.SingularMatrixException
165 * if the covariance matrix cannot be computed.
166 */
167 public double[] computeSigma(double[] params,
168 double covarianceSingularityThreshold) {
169 final int nC = params.length;
170 final double[] sig = new double[nC];
171 final double[][] cov = computeCovariances(params, covarianceSingularityThreshold);
172 for (int i = 0; i < nC; ++i) {
173 sig[i] = FastMath.sqrt(cov[i][i]);
174 }
175 return sig;
176 }
177
178 /**
179 * {@inheritDoc}
180 *
181 * @param optData Optimization data. The following data will be looked for:
182 * <ul>
183 * <li>{@link org.apache.commons.math3.optim.MaxEval}</li>
184 * <li>{@link org.apache.commons.math3.optim.InitialGuess}</li>
185 * <li>{@link org.apache.commons.math3.optim.SimpleBounds}</li>
186 * <li>{@link org.apache.commons.math3.optim.nonlinear.vector.Target}</li>
187 * <li>{@link org.apache.commons.math3.optim.nonlinear.vector.Weight}</li>
188 * <li>{@link org.apache.commons.math3.optim.nonlinear.vector.ModelFunction}</li>
189 * <li>{@link org.apache.commons.math3.optim.nonlinear.vector.ModelFunctionJacobian}</li>
190 * </ul>
191 * @return {@inheritDoc}
192 * @throws TooManyEvaluationsException if the maximal number of
193 * evaluations is exceeded.
194 * @throws DimensionMismatchException if the initial guess, target, and weight
195 * arguments have inconsistent dimensions.
196 */
197 @Override
198 public PointVectorValuePair optimize(OptimizationData... optData)
199 throws TooManyEvaluationsException {
200 // Retrieve settings.
201 parseOptimizationData(optData);
202 // Set up base class and perform computation.
203 return super.optimize(optData);
204 }
205
206 /**
207 * Computes the residuals.
208 * The residual is the difference between the observed (target)
209 * values and the model (objective function) value.
210 * There is one residual for each element of the vector-valued
211 * function.
212 *
213 * @param objectiveValue Value of the the objective function. This is
214 * the value returned from a call to
215 * {@link #computeObjectiveValue(double[]) computeObjectiveValue}
216 * (whose array argument contains the model parameters).
217 * @return the residuals.
218 * @throws DimensionMismatchException if {@code params} has a wrong
219 * length.
220 */
221 protected double[] computeResiduals(double[] objectiveValue) {
222 final double[] target = getTarget();
223 if (objectiveValue.length != target.length) {
224 throw new DimensionMismatchException(target.length,
225 objectiveValue.length);
226 }
227
228 final double[] residuals = new double[target.length];
229 for (int i = 0; i < target.length; i++) {
230 residuals[i] = target[i] - objectiveValue[i];
231 }
232
233 return residuals;
234 }
235
236 /**
237 * Scans the list of (required and optional) optimization data that
238 * characterize the problem.
239 * If the weight matrix is specified, the {@link #weightMatrixSqrt}
240 * field is recomputed.
241 *
242 * @param optData Optimization data. The following data will be looked for:
243 * <ul>
244 * <li>{@link Weight}</li>
245 * </ul>
246 */
247 private void parseOptimizationData(OptimizationData... optData) {
248 // The existing values (as set by the previous call) are reused if
249 // not provided in the argument list.
250 for (OptimizationData data : optData) {
251 if (data instanceof Weight) {
252 weightMatrixSqrt = squareRoot(((Weight) data).getWeight());
253 // If more data must be parsed, this statement _must_ be
254 // changed to "continue".
255 break;
256 }
257 }
258 }
259
260 /**
261 * Computes the square-root of the weight matrix.
262 *
263 * @param m Symmetric, positive-definite (weight) matrix.
264 * @return the square-root of the weight matrix.
265 */
266 private RealMatrix squareRoot(RealMatrix m) {
267 if (m instanceof DiagonalMatrix) {
268 final int dim = m.getRowDimension();
269 final RealMatrix sqrtM = new DiagonalMatrix(dim);
270 for (int i = 0; i < dim; i++) {
271 sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i)));
272 }
273 return sqrtM;
274 } else {
275 final EigenDecomposition dec = new EigenDecomposition(m);
276 return dec.getSquareRoot();
277 }
278 }
279 }