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.linear.InvalidMatrixException;
23 import org.apache.commons.math.linear.LUDecompositionImpl;
24 import org.apache.commons.math.linear.MatrixUtils;
25 import org.apache.commons.math.linear.RealMatrix;
26 import org.apache.commons.math.linear.RealVector;
27 import org.apache.commons.math.linear.ArrayRealVector;
28
29 /**
30 * This class implements a solver for estimation problems.
31 *
32 * <p>This class solves estimation problems using a weighted least
33 * squares criterion on the measurement residuals. It uses a
34 * Gauss-Newton algorithm.</p>
35 *
36 * @version $Revision: 783702 $ $Date: 2009-06-11 04:54:02 -0400 (Thu, 11 Jun 2009) $
37 * @since 1.2
38 * @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
39 * been deprecated and replaced by package org.apache.commons.math.optimization.general
40 *
41 */
42 @Deprecated
43 public class GaussNewtonEstimator extends AbstractEstimator implements Serializable {
44
45 /** Serializable version identifier */
46 private static final long serialVersionUID = 5485001826076289109L;
47
48 /** Default threshold for cost steady state detection. */
49 private static final double DEFAULT_STEADY_STATE_THRESHOLD = 1.0e-6;
50
51 /** Default threshold for cost convergence. */
52 private static final double DEFAULT_CONVERGENCE = 1.0e-6;
53
54 /** Threshold for cost steady state detection. */
55 private double steadyStateThreshold;
56
57 /** Threshold for cost convergence. */
58 private double convergence;
59
60 /** Simple constructor with default settings.
61 * <p>
62 * The estimator is built with default values for all settings.
63 * </p>
64 * @see #DEFAULT_STEADY_STATE_THRESHOLD
65 * @see #DEFAULT_CONVERGENCE
66 * @see AbstractEstimator#DEFAULT_MAX_COST_EVALUATIONS
67 */
68 public GaussNewtonEstimator() {
69 this.steadyStateThreshold = DEFAULT_STEADY_STATE_THRESHOLD;
70 this.convergence = DEFAULT_CONVERGENCE;
71 }
72
73 /**
74 * Simple constructor.
75 *
76 * <p>This constructor builds an estimator and stores its convergence
77 * characteristics.</p>
78 *
79 * <p>An estimator is considered to have converged whenever either
80 * the criterion goes below a physical threshold under which
81 * improvements are considered useless or when the algorithm is
82 * unable to improve it (even if it is still high). The first
83 * condition that is met stops the iterations.</p>
84 *
85 * <p>The fact an estimator has converged does not mean that the
86 * model accurately fits the measurements. It only means no better
87 * solution can be found, it does not mean this one is good. Such an
88 * analysis is left to the caller.</p>
89 *
90 * <p>If neither conditions are fulfilled before a given number of
91 * iterations, the algorithm is considered to have failed and an
92 * {@link EstimationException} is thrown.</p>
93 *
94 * @param maxCostEval maximal number of cost evaluations allowed
95 * @param convergence criterion threshold below which we do not need
96 * to improve the criterion anymore
97 * @param steadyStateThreshold steady state detection threshold, the
98 * problem has converged has reached a steady state if
99 * <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 }