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