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.math.optimization.general;
019
020 import org.apache.commons.math.exception.NumberIsTooSmallException;
021 import org.apache.commons.math.exception.DimensionMismatchException;
022 import org.apache.commons.math.analysis.DifferentiableMultivariateVectorialFunction;
023 import org.apache.commons.math.analysis.MultivariateMatrixFunction;
024 import org.apache.commons.math.exception.util.LocalizedFormats;
025 import org.apache.commons.math.linear.LUDecomposition;
026 import org.apache.commons.math.linear.DecompositionSolver;
027 import org.apache.commons.math.linear.MatrixUtils;
028 import org.apache.commons.math.optimization.ConvergenceChecker;
029 import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer;
030 import org.apache.commons.math.optimization.VectorialPointValuePair;
031 import org.apache.commons.math.optimization.direct.BaseAbstractVectorialOptimizer;
032 import org.apache.commons.math.util.FastMath;
033
034 /**
035 * Base class for implementing least squares optimizers.
036 * It handles the boilerplate methods associated to thresholds settings,
037 * jacobian and error estimation.
038 * <br/>
039 * This class uses the {@link DifferentiableMultivariateVectorialFunction#jacobian()}
040 * of the function argument in method
041 * {@link #optimize(int,DifferentiableMultivariateVectorialFunction,double[],double[],double[])
042 * optimize} and assumes that, in the matrix returned by the
043 * {@link MultivariateMatrixFunction#value(double[]) value} method, the rows
044 * iterate on the model functions while the columns iterate on the parameters; thus,
045 * the numbers of rows is equal to the length of the {@code target} array while the
046 * number of columns is equal to the length of the {@code startPoint} array.
047 *
048 * @version $Id: AbstractLeastSquaresOptimizer.java 1175099 2011-09-24 04:28:36Z celestin $
049 * @since 1.2
050 */
051 public abstract class AbstractLeastSquaresOptimizer
052 extends BaseAbstractVectorialOptimizer<DifferentiableMultivariateVectorialFunction>
053 implements DifferentiableMultivariateVectorialOptimizer {
054 /** Singularity threshold (cf. {@link #getCovariances(double)}). */
055 private static final double DEFAULT_SINGULARITY_THRESHOLD = 1e-14;
056 /**
057 * Jacobian matrix of the weighted residuals.
058 * This matrix is in canonical form just after the calls to
059 * {@link #updateJacobian()}, but may be modified by the solver
060 * in the derived class (the {@link LevenbergMarquardtOptimizer
061 * Levenberg-Marquardt optimizer} does this).
062 */
063 protected double[][] weightedResidualJacobian;
064 /** Number of columns of the jacobian matrix. */
065 protected int cols;
066 /** Number of rows of the jacobian matrix. */
067 protected int rows;
068 /** Current point. */
069 protected double[] point;
070 /** Current objective function value. */
071 protected double[] objective;
072 /** Current residuals. */
073 protected double[] residuals;
074 /** Weighted residuals */
075 protected double[] weightedResiduals;
076 /** Cost value (square root of the sum of the residuals). */
077 protected double cost;
078 /** Objective function derivatives. */
079 private MultivariateMatrixFunction jF;
080 /** Number of evaluations of the Jacobian. */
081 private int jacobianEvaluations;
082
083 /**
084 * Simple constructor with default settings.
085 * The convergence check is set to a {@link
086 * org.apache.commons.math.optimization.SimpleVectorialValueChecker}.
087 */
088 protected AbstractLeastSquaresOptimizer() {}
089 /**
090 * @param checker Convergence checker.
091 */
092 protected AbstractLeastSquaresOptimizer(ConvergenceChecker<VectorialPointValuePair> checker) {
093 super(checker);
094 }
095
096 /**
097 * @return the number of evaluations of the Jacobian function.
098 */
099 public int getJacobianEvaluations() {
100 return jacobianEvaluations;
101 }
102
103 /**
104 * Update the jacobian matrix.
105 *
106 * @throws DimensionMismatchException if the Jacobian dimension does not
107 * match problem dimension.
108 */
109 protected void updateJacobian() {
110 ++jacobianEvaluations;
111 weightedResidualJacobian = jF.value(point);
112 if (weightedResidualJacobian.length != rows) {
113 throw new DimensionMismatchException(weightedResidualJacobian.length, rows);
114 }
115
116 final double[] residualsWeights = getWeightRef();
117
118 for (int i = 0; i < rows; i++) {
119 final double[] ji = weightedResidualJacobian[i];
120 double wi = FastMath.sqrt(residualsWeights[i]);
121 for (int j = 0; j < cols; ++j) {
122 //ji[j] *= -1.0;
123 weightedResidualJacobian[i][j] = -ji[j]*wi;
124 }
125 }
126 }
127
128 /**
129 * Update the residuals array and cost function value.
130 * @throws DimensionMismatchException if the dimension does not match the
131 * problem dimension.
132 * @throws org.apache.commons.math.exception.TooManyEvaluationsException
133 * if the maximal number of evaluations is exceeded.
134 */
135 protected void updateResidualsAndCost() {
136 objective = computeObjectiveValue(point);
137 if (objective.length != rows) {
138 throw new DimensionMismatchException(objective.length, rows);
139 }
140
141 final double[] targetValues = getTargetRef();
142 final double[] residualsWeights = getWeightRef();
143
144 cost = 0;
145 for (int i = 0; i < rows; i++) {
146 final double residual = targetValues[i] - objective[i];
147 weightedResiduals[i]= residual*FastMath.sqrt(residualsWeights[i]);
148 cost += residualsWeights[i] * residual * residual;
149 }
150 cost = FastMath.sqrt(cost);
151 }
152
153 /**
154 * Get the Root Mean Square value.
155 * Get the Root Mean Square value, i.e. the root of the arithmetic
156 * mean of the square of all weighted residuals. This is related to the
157 * criterion that is minimized by the optimizer as follows: if
158 * <em>c</em> if the criterion, and <em>n</em> is the number of
159 * measurements, then the RMS is <em>sqrt (c/n)</em>.
160 *
161 * @return RMS value
162 */
163 public double getRMS() {
164 return FastMath.sqrt(getChiSquare() / rows);
165 }
166
167 /**
168 * Get a Chi-Square-like value assuming the N residuals follow N
169 * distinct normal distributions centered on 0 and whose variances are
170 * the reciprocal of the weights.
171 * @return chi-square value
172 */
173 public double getChiSquare() {
174 return cost * cost;
175 }
176
177 /**
178 * Get the covariance matrix of the optimized parameters.
179 *
180 * @return the covariance matrix.
181 * @throws org.apache.commons.math.linear.SingularMatrixException
182 * if the covariance matrix cannot be computed (singular problem).
183 */
184 public double[][] getCovariances() {
185 return getCovariances(DEFAULT_SINGULARITY_THRESHOLD);
186 }
187
188 /**
189 * Get the covariance matrix of the optimized parameters.
190 *
191 * @param threshold Singularity threshold.
192 * @return the covariance matrix.
193 * @throws org.apache.commons.math.linear.SingularMatrixException
194 * if the covariance matrix cannot be computed (singular problem).
195 */
196 public double[][] getCovariances(double threshold) {
197 // Set up the jacobian.
198 updateJacobian();
199
200 // Compute transpose(J)J, without building intermediate matrices.
201 double[][] jTj = new double[cols][cols];
202 for (int i = 0; i < cols; ++i) {
203 for (int j = i; j < cols; ++j) {
204 double sum = 0;
205 for (int k = 0; k < rows; ++k) {
206 sum += weightedResidualJacobian[k][i] * weightedResidualJacobian[k][j];
207 }
208 jTj[i][j] = sum;
209 jTj[j][i] = sum;
210 }
211 }
212
213 // Compute the covariances matrix.
214 final DecompositionSolver solver
215 = new LUDecomposition(MatrixUtils.createRealMatrix(jTj), threshold).getSolver();
216 return solver.getInverse().getData();
217 }
218
219 /**
220 * Guess the errors in optimized parameters.
221 * Guessing is covariance-based: It only gives a rough order of magnitude.
222 *
223 * @return errors in optimized parameters
224 * @throws org.apache.commons.math.linear.SingularMatrixException
225 * if the covariances matrix cannot be computed.
226 * @throws NumberIsTooSmallException if the number of degrees of freedom is not
227 * positive, i.e. the number of measurements is less or equal to the number of
228 * parameters.
229 */
230 public double[] guessParametersErrors() {
231 if (rows <= cols) {
232 throw new NumberIsTooSmallException(LocalizedFormats.NO_DEGREES_OF_FREEDOM,
233 rows, cols, false);
234 }
235 double[] errors = new double[cols];
236 final double c = FastMath.sqrt(getChiSquare() / (rows - cols));
237 double[][] covar = getCovariances();
238 for (int i = 0; i < errors.length; ++i) {
239 errors[i] = FastMath.sqrt(covar[i][i]) * c;
240 }
241 return errors;
242 }
243
244 /** {@inheritDoc} */
245 @Override
246 public VectorialPointValuePair optimize(int maxEval,
247 final DifferentiableMultivariateVectorialFunction f,
248 final double[] target, final double[] weights,
249 final double[] startPoint) {
250 // Reset counter.
251 jacobianEvaluations = 0;
252
253 // Store least squares problem characteristics.
254 jF = f.jacobian();
255 this.residuals = new double[target.length];
256
257 // Arrays shared with the other private methods.
258 point = startPoint.clone();
259 rows = target.length;
260 cols = point.length;
261
262 weightedResidualJacobian = new double[rows][cols];
263 this.weightedResiduals = new double[rows];
264
265 cost = Double.POSITIVE_INFINITY;
266
267 return super.optimize(maxEval, f, target, weights, startPoint);
268 }
269 }