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 }