1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 18 package org.apache.commons.math.estimation; 19 20 import java.io.Serializable; 21 22 import org.apache.commons.math.exception.util.LocalizedFormats; 23 import org.apache.commons.math.linear.InvalidMatrixException; 24 import org.apache.commons.math.linear.LUDecompositionImpl; 25 import org.apache.commons.math.linear.MatrixUtils; 26 import org.apache.commons.math.linear.RealMatrix; 27 import org.apache.commons.math.linear.RealVector; 28 import org.apache.commons.math.linear.ArrayRealVector; 29 import org.apache.commons.math.util.FastMath; 30 31 /** 32 * This class implements a solver for estimation problems. 33 * 34 * <p>This class solves estimation problems using a weighted least 35 * squares criterion on the measurement residuals. It uses a 36 * Gauss-Newton algorithm.</p> 37 * 38 * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 aot 2010) $ 39 * @since 1.2 40 * @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has 41 * been deprecated and replaced by package org.apache.commons.math.optimization.general 42 * 43 */ 44 @Deprecated 45 public class GaussNewtonEstimator extends AbstractEstimator implements Serializable { 46 47 /** Serializable version identifier */ 48 private static final long serialVersionUID = 5485001826076289109L; 49 50 /** Default threshold for cost steady state detection. */ 51 private static final double DEFAULT_STEADY_STATE_THRESHOLD = 1.0e-6; 52 53 /** Default threshold for cost convergence. */ 54 private static final double DEFAULT_CONVERGENCE = 1.0e-6; 55 56 /** Threshold for cost steady state detection. */ 57 private double steadyStateThreshold; 58 59 /** Threshold for cost convergence. */ 60 private double convergence; 61 62 /** Simple constructor with default settings. 63 * <p> 64 * The estimator is built with default values for all settings. 65 * </p> 66 * @see #DEFAULT_STEADY_STATE_THRESHOLD 67 * @see #DEFAULT_CONVERGENCE 68 * @see AbstractEstimator#DEFAULT_MAX_COST_EVALUATIONS 69 */ 70 public GaussNewtonEstimator() { 71 this.steadyStateThreshold = DEFAULT_STEADY_STATE_THRESHOLD; 72 this.convergence = DEFAULT_CONVERGENCE; 73 } 74 75 /** 76 * Simple constructor. 77 * 78 * <p>This constructor builds an estimator and stores its convergence 79 * characteristics.</p> 80 * 81 * <p>An estimator is considered to have converged whenever either 82 * the criterion goes below a physical threshold under which 83 * improvements are considered useless or when the algorithm is 84 * unable to improve it (even if it is still high). The first 85 * condition that is met stops the iterations.</p> 86 * 87 * <p>The fact an estimator has converged does not mean that the 88 * model accurately fits the measurements. It only means no better 89 * solution can be found, it does not mean this one is good. Such an 90 * analysis is left to the caller.</p> 91 * 92 * <p>If neither conditions are fulfilled before a given number of 93 * iterations, the algorithm is considered to have failed and an 94 * {@link EstimationException} is thrown.</p> 95 * 96 * @param maxCostEval maximal number of cost evaluations allowed 97 * @param convergence criterion threshold below which we do not need 98 * to improve the criterion anymore 99 * @param steadyStateThreshold steady state detection threshold, the 100 * problem has converged has reached a steady state if 101 * <code>FastMath.abs(J<sub>n</sub> - J<sub>n-1</sub>) < 102 * J<sub>n</sub> × convergence</code>, where <code>J<sub>n</sub></code> 103 * and <code>J<sub>n-1</sub></code> are the current and preceding criterion 104 * values (square sum of the weighted residuals of considered measurements). 105 */ 106 public GaussNewtonEstimator(final int maxCostEval, final double convergence, 107 final double steadyStateThreshold) { 108 setMaxCostEval(maxCostEval); 109 this.steadyStateThreshold = steadyStateThreshold; 110 this.convergence = convergence; 111 } 112 113 /** 114 * Set the convergence criterion threshold. 115 * @param convergence criterion threshold below which we do not need 116 * to improve the criterion anymore 117 */ 118 public void setConvergence(final double convergence) { 119 this.convergence = convergence; 120 } 121 122 /** 123 * Set the steady state detection threshold. 124 * <p> 125 * The problem has converged has reached a steady state if 126 * <code>FastMath.abs(J<sub>n</sub> - J<sub>n-1</sub>) < 127 * J<sub>n</sub> × convergence</code>, where <code>J<sub>n</sub></code> 128 * and <code>J<sub>n-1</sub></code> are the current and preceding criterion 129 * values (square sum of the weighted residuals of considered measurements). 130 * </p> 131 * @param steadyStateThreshold steady state detection threshold 132 */ 133 public void setSteadyStateThreshold(final double steadyStateThreshold) { 134 this.steadyStateThreshold = steadyStateThreshold; 135 } 136 137 /** 138 * Solve an estimation problem using a least squares criterion. 139 * 140 * <p>This method set the unbound parameters of the given problem 141 * starting from their current values through several iterations. At 142 * each step, the unbound parameters are changed in order to 143 * minimize a weighted least square criterion based on the 144 * measurements of the problem.</p> 145 * 146 * <p>The iterations are stopped either when the criterion goes 147 * below a physical threshold under which improvement are considered 148 * useless or when the algorithm is unable to improve it (even if it 149 * is still high). The first condition that is met stops the 150 * iterations. If the convergence it not reached before the maximum 151 * number of iterations, an {@link EstimationException} is 152 * thrown.</p> 153 * 154 * @param problem estimation problem to solve 155 * @exception EstimationException if the problem cannot be solved 156 * 157 * @see EstimationProblem 158 * 159 */ 160 @Override 161 public void estimate(EstimationProblem problem) 162 throws EstimationException { 163 164 initializeEstimate(problem); 165 166 // work matrices 167 double[] grad = new double[parameters.length]; 168 ArrayRealVector bDecrement = new ArrayRealVector(parameters.length); 169 double[] bDecrementData = bDecrement.getDataRef(); 170 RealMatrix wGradGradT = MatrixUtils.createRealMatrix(parameters.length, parameters.length); 171 172 // iterate until convergence is reached 173 double previous = Double.POSITIVE_INFINITY; 174 do { 175 176 // build the linear problem 177 incrementJacobianEvaluationsCounter(); 178 RealVector b = new ArrayRealVector(parameters.length); 179 RealMatrix a = MatrixUtils.createRealMatrix(parameters.length, parameters.length); 180 for (int i = 0; i < measurements.length; ++i) { 181 if (! measurements [i].isIgnored()) { 182 183 double weight = measurements[i].getWeight(); 184 double residual = measurements[i].getResidual(); 185 186 // compute the normal equation 187 for (int j = 0; j < parameters.length; ++j) { 188 grad[j] = measurements[i].getPartial(parameters[j]); 189 bDecrementData[j] = weight * residual * grad[j]; 190 } 191 192 // build the contribution matrix for measurement i 193 for (int k = 0; k < parameters.length; ++k) { 194 double gk = grad[k]; 195 for (int l = 0; l < parameters.length; ++l) { 196 wGradGradT.setEntry(k, l, weight * gk * grad[l]); 197 } 198 } 199 200 // update the matrices 201 a = a.add(wGradGradT); 202 b = b.add(bDecrement); 203 204 } 205 } 206 207 try { 208 209 // solve the linearized least squares problem 210 RealVector dX = new LUDecompositionImpl(a).getSolver().solve(b); 211 212 // update the estimated parameters 213 for (int i = 0; i < parameters.length; ++i) { 214 parameters[i].setEstimate(parameters[i].getEstimate() + dX.getEntry(i)); 215 } 216 217 } catch(InvalidMatrixException e) { 218 throw new EstimationException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM); 219 } 220 221 222 previous = cost; 223 updateResidualsAndCost(); 224 225 } while ((getCostEvaluations() < 2) || 226 (FastMath.abs(previous - cost) > (cost * steadyStateThreshold) && 227 (FastMath.abs(cost) > convergence))); 228 229 } 230 231 } 232