View Javadoc

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>) &lt;
100      * J<sub>n</sub> &times 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>) &lt;
125      * J<sub>n</sub> &times 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 }