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.estimation; 019 020 import java.io.Serializable; 021 022 import org.apache.commons.math.linear.InvalidMatrixException; 023 import org.apache.commons.math.linear.LUDecompositionImpl; 024 import org.apache.commons.math.linear.MatrixUtils; 025 import org.apache.commons.math.linear.RealMatrix; 026 import org.apache.commons.math.linear.RealVector; 027 import org.apache.commons.math.linear.ArrayRealVector; 028 029 /** 030 * This class implements a solver for estimation problems. 031 * 032 * <p>This class solves estimation problems using a weighted least 033 * squares criterion on the measurement residuals. It uses a 034 * Gauss-Newton algorithm.</p> 035 * 036 * @version $Revision: 783702 $ $Date: 2009-06-11 04:54:02 -0400 (Thu, 11 Jun 2009) $ 037 * @since 1.2 038 * @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has 039 * been deprecated and replaced by package org.apache.commons.math.optimization.general 040 * 041 */ 042 @Deprecated 043 public class GaussNewtonEstimator extends AbstractEstimator implements Serializable { 044 045 /** Serializable version identifier */ 046 private static final long serialVersionUID = 5485001826076289109L; 047 048 /** Default threshold for cost steady state detection. */ 049 private static final double DEFAULT_STEADY_STATE_THRESHOLD = 1.0e-6; 050 051 /** Default threshold for cost convergence. */ 052 private static final double DEFAULT_CONVERGENCE = 1.0e-6; 053 054 /** Threshold for cost steady state detection. */ 055 private double steadyStateThreshold; 056 057 /** Threshold for cost convergence. */ 058 private double convergence; 059 060 /** Simple constructor with default settings. 061 * <p> 062 * The estimator is built with default values for all settings. 063 * </p> 064 * @see #DEFAULT_STEADY_STATE_THRESHOLD 065 * @see #DEFAULT_CONVERGENCE 066 * @see AbstractEstimator#DEFAULT_MAX_COST_EVALUATIONS 067 */ 068 public GaussNewtonEstimator() { 069 this.steadyStateThreshold = DEFAULT_STEADY_STATE_THRESHOLD; 070 this.convergence = DEFAULT_CONVERGENCE; 071 } 072 073 /** 074 * Simple constructor. 075 * 076 * <p>This constructor builds an estimator and stores its convergence 077 * characteristics.</p> 078 * 079 * <p>An estimator is considered to have converged whenever either 080 * the criterion goes below a physical threshold under which 081 * improvements are considered useless or when the algorithm is 082 * unable to improve it (even if it is still high). The first 083 * condition that is met stops the iterations.</p> 084 * 085 * <p>The fact an estimator has converged does not mean that the 086 * model accurately fits the measurements. It only means no better 087 * solution can be found, it does not mean this one is good. Such an 088 * analysis is left to the caller.</p> 089 * 090 * <p>If neither conditions are fulfilled before a given number of 091 * iterations, the algorithm is considered to have failed and an 092 * {@link EstimationException} is thrown.</p> 093 * 094 * @param maxCostEval maximal number of cost evaluations allowed 095 * @param convergence criterion threshold below which we do not need 096 * to improve the criterion anymore 097 * @param steadyStateThreshold steady state detection threshold, the 098 * problem has converged has reached a steady state if 099 * <code>Math.abs(J<sub>n</sub> - J<sub>n-1</sub>) < 100 * J<sub>n</sub> × convergence</code>, where <code>J<sub>n</sub></code> 101 * and <code>J<sub>n-1</sub></code> are the current and preceding criterion 102 * values (square sum of the weighted residuals of considered measurements). 103 */ 104 public GaussNewtonEstimator(final int maxCostEval, final double convergence, 105 final double steadyStateThreshold) { 106 setMaxCostEval(maxCostEval); 107 this.steadyStateThreshold = steadyStateThreshold; 108 this.convergence = convergence; 109 } 110 111 /** 112 * Set the convergence criterion threshold. 113 * @param convergence criterion threshold below which we do not need 114 * to improve the criterion anymore 115 */ 116 public void setConvergence(final double convergence) { 117 this.convergence = convergence; 118 } 119 120 /** 121 * Set the steady state detection threshold. 122 * <p> 123 * The problem has converged has reached a steady state if 124 * <code>Math.abs(J<sub>n</sub> - J<sub>n-1</sub>) < 125 * J<sub>n</sub> × convergence</code>, where <code>J<sub>n</sub></code> 126 * and <code>J<sub>n-1</sub></code> are the current and preceding criterion 127 * values (square sum of the weighted residuals of considered measurements). 128 * </p> 129 * @param steadyStateThreshold steady state detection threshold 130 */ 131 public void setSteadyStateThreshold(final double steadyStateThreshold) { 132 this.steadyStateThreshold = steadyStateThreshold; 133 } 134 135 /** 136 * Solve an estimation problem using a least squares criterion. 137 * 138 * <p>This method set the unbound parameters of the given problem 139 * starting from their current values through several iterations. At 140 * each step, the unbound parameters are changed in order to 141 * minimize a weighted least square criterion based on the 142 * measurements of the problem.</p> 143 * 144 * <p>The iterations are stopped either when the criterion goes 145 * below a physical threshold under which improvement are considered 146 * useless or when the algorithm is unable to improve it (even if it 147 * is still high). The first condition that is met stops the 148 * iterations. If the convergence it not reached before the maximum 149 * number of iterations, an {@link EstimationException} is 150 * thrown.</p> 151 * 152 * @param problem estimation problem to solve 153 * @exception EstimationException if the problem cannot be solved 154 * 155 * @see EstimationProblem 156 * 157 */ 158 @Override 159 public void estimate(EstimationProblem problem) 160 throws EstimationException { 161 162 initializeEstimate(problem); 163 164 // work matrices 165 double[] grad = new double[parameters.length]; 166 ArrayRealVector bDecrement = new ArrayRealVector(parameters.length); 167 double[] bDecrementData = bDecrement.getDataRef(); 168 RealMatrix wGradGradT = MatrixUtils.createRealMatrix(parameters.length, parameters.length); 169 170 // iterate until convergence is reached 171 double previous = Double.POSITIVE_INFINITY; 172 do { 173 174 // build the linear problem 175 incrementJacobianEvaluationsCounter(); 176 RealVector b = new ArrayRealVector(parameters.length); 177 RealMatrix a = MatrixUtils.createRealMatrix(parameters.length, parameters.length); 178 for (int i = 0; i < measurements.length; ++i) { 179 if (! measurements [i].isIgnored()) { 180 181 double weight = measurements[i].getWeight(); 182 double residual = measurements[i].getResidual(); 183 184 // compute the normal equation 185 for (int j = 0; j < parameters.length; ++j) { 186 grad[j] = measurements[i].getPartial(parameters[j]); 187 bDecrementData[j] = weight * residual * grad[j]; 188 } 189 190 // build the contribution matrix for measurement i 191 for (int k = 0; k < parameters.length; ++k) { 192 double gk = grad[k]; 193 for (int l = 0; l < parameters.length; ++l) { 194 wGradGradT.setEntry(k, l, weight * gk * grad[l]); 195 } 196 } 197 198 // update the matrices 199 a = a.add(wGradGradT); 200 b = b.add(bDecrement); 201 202 } 203 } 204 205 try { 206 207 // solve the linearized least squares problem 208 RealVector dX = new LUDecompositionImpl(a).getSolver().solve(b); 209 210 // update the estimated parameters 211 for (int i = 0; i < parameters.length; ++i) { 212 parameters[i].setEstimate(parameters[i].getEstimate() + dX.getEntry(i)); 213 } 214 215 } catch(InvalidMatrixException e) { 216 throw new EstimationException("unable to solve: singular problem"); 217 } 218 219 220 previous = cost; 221 updateResidualsAndCost(); 222 223 } while ((getCostEvaluations() < 2) || 224 (Math.abs(previous - cost) > (cost * steadyStateThreshold) && 225 (Math.abs(cost) > convergence))); 226 227 } 228 229 }