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.exception.util.LocalizedFormats; 023 import org.apache.commons.math.linear.InvalidMatrixException; 024 import org.apache.commons.math.linear.LUDecompositionImpl; 025 import org.apache.commons.math.linear.MatrixUtils; 026 import org.apache.commons.math.linear.RealMatrix; 027 import org.apache.commons.math.linear.RealVector; 028 import org.apache.commons.math.linear.ArrayRealVector; 029 import org.apache.commons.math.util.FastMath; 030 031 /** 032 * This class implements a solver for estimation problems. 033 * 034 * <p>This class solves estimation problems using a weighted least 035 * squares criterion on the measurement residuals. It uses a 036 * Gauss-Newton algorithm.</p> 037 * 038 * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 ao??t 2010) $ 039 * @since 1.2 040 * @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has 041 * been deprecated and replaced by package org.apache.commons.math.optimization.general 042 * 043 */ 044 @Deprecated 045 public class GaussNewtonEstimator extends AbstractEstimator implements Serializable { 046 047 /** Serializable version identifier */ 048 private static final long serialVersionUID = 5485001826076289109L; 049 050 /** Default threshold for cost steady state detection. */ 051 private static final double DEFAULT_STEADY_STATE_THRESHOLD = 1.0e-6; 052 053 /** Default threshold for cost convergence. */ 054 private static final double DEFAULT_CONVERGENCE = 1.0e-6; 055 056 /** Threshold for cost steady state detection. */ 057 private double steadyStateThreshold; 058 059 /** Threshold for cost convergence. */ 060 private double convergence; 061 062 /** Simple constructor with default settings. 063 * <p> 064 * The estimator is built with default values for all settings. 065 * </p> 066 * @see #DEFAULT_STEADY_STATE_THRESHOLD 067 * @see #DEFAULT_CONVERGENCE 068 * @see AbstractEstimator#DEFAULT_MAX_COST_EVALUATIONS 069 */ 070 public GaussNewtonEstimator() { 071 this.steadyStateThreshold = DEFAULT_STEADY_STATE_THRESHOLD; 072 this.convergence = DEFAULT_CONVERGENCE; 073 } 074 075 /** 076 * Simple constructor. 077 * 078 * <p>This constructor builds an estimator and stores its convergence 079 * characteristics.</p> 080 * 081 * <p>An estimator is considered to have converged whenever either 082 * the criterion goes below a physical threshold under which 083 * improvements are considered useless or when the algorithm is 084 * unable to improve it (even if it is still high). The first 085 * condition that is met stops the iterations.</p> 086 * 087 * <p>The fact an estimator has converged does not mean that the 088 * model accurately fits the measurements. It only means no better 089 * solution can be found, it does not mean this one is good. Such an 090 * analysis is left to the caller.</p> 091 * 092 * <p>If neither conditions are fulfilled before a given number of 093 * iterations, the algorithm is considered to have failed and an 094 * {@link EstimationException} is thrown.</p> 095 * 096 * @param maxCostEval maximal number of cost evaluations allowed 097 * @param convergence criterion threshold below which we do not need 098 * to improve the criterion anymore 099 * @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 }