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.optimization.general;
19  
20  import org.apache.commons.math.FunctionEvaluationException;
21  import org.apache.commons.math.linear.BlockRealMatrix;
22  import org.apache.commons.math.linear.DecompositionSolver;
23  import org.apache.commons.math.linear.InvalidMatrixException;
24  import org.apache.commons.math.linear.LUDecompositionImpl;
25  import org.apache.commons.math.linear.QRDecompositionImpl;
26  import org.apache.commons.math.linear.RealMatrix;
27  import org.apache.commons.math.optimization.OptimizationException;
28  import org.apache.commons.math.optimization.SimpleVectorialValueChecker;
29  import org.apache.commons.math.optimization.VectorialPointValuePair;
30  
31  /** 
32   * Gauss-Newton least-squares solver.
33   * <p>
34   * This class solve a least-square problem by solving the normal equations
35   * of the linearized problem at each iteration. Either LU decomposition or
36   * QR decomposition can be used to solve the normal equations. LU decomposition
37   * is faster but QR decomposition is more robust for difficult problems.
38   * </p>
39   *
40   * @version $Revision: 795978 $ $Date: 2009-07-20 15:57:08 -0400 (Mon, 20 Jul 2009) $
41   * @since 2.0
42   *
43   */
44  
45  public class GaussNewtonOptimizer extends AbstractLeastSquaresOptimizer {
46  
47      /** Indicator for using LU decomposition. */
48      private final boolean useLU;
49  
50      /** Simple constructor with default settings.
51       * <p>The convergence check is set to a {@link SimpleVectorialValueChecker}
52       * and the maximal number of evaluation is set to
53       * {@link AbstractLeastSquaresOptimizer#DEFAULT_MAX_ITERATIONS}.
54       * @param useLU if true, the normal equations will be solved using LU
55       * decomposition, otherwise they will be solved using QR decomposition
56       */
57      public GaussNewtonOptimizer(final boolean useLU) {
58          this.useLU = useLU;
59      }
60  
61      /** {@inheritDoc} */
62      @Override
63      public VectorialPointValuePair doOptimize()
64          throws FunctionEvaluationException, OptimizationException, IllegalArgumentException {
65  
66          // iterate until convergence is reached
67          VectorialPointValuePair current = null;
68          for (boolean converged = false; !converged;) {
69  
70              incrementIterationsCounter();
71  
72              // evaluate the objective function and its jacobian
73              VectorialPointValuePair previous = current;
74              updateResidualsAndCost();
75              updateJacobian();
76              current = new VectorialPointValuePair(point, objective);
77  
78              // build the linear problem
79              final double[]   b = new double[cols];
80              final double[][] a = new double[cols][cols];
81              for (int i = 0; i < rows; ++i) {
82  
83                  final double[] grad   = jacobian[i];
84                  final double weight   = weights[i];
85                  final double residual = objective[i] - target[i];
86  
87                  // compute the normal equation
88                  final double wr = weight * residual;
89                  for (int j = 0; j < cols; ++j) {
90                      b[j] += wr * grad[j];
91                  }
92  
93                  // build the contribution matrix for measurement i
94                  for (int k = 0; k < cols; ++k) {
95                      double[] ak = a[k];
96                      double wgk = weight * grad[k];
97                      for (int l = 0; l < cols; ++l) {
98                          ak[l] += wgk * grad[l];
99                      }
100                 }
101 
102             }
103 
104             try {
105 
106                 // solve the linearized least squares problem
107                 RealMatrix mA = new BlockRealMatrix(a);
108                 DecompositionSolver solver = useLU ?
109                         new LUDecompositionImpl(mA).getSolver() :
110                         new QRDecompositionImpl(mA).getSolver();
111                 final double[] dX = solver.solve(b);
112 
113                 // update the estimated parameters
114                 for (int i = 0; i < cols; ++i) {
115                     point[i] += dX[i];
116                 }
117 
118             } catch(InvalidMatrixException e) {
119                 throw new OptimizationException("unable to solve: singular problem");
120             }
121 
122             // check convergence
123             if (previous != null) {
124                 converged = checker.converged(getIterations(), previous, current);
125             }
126 
127         }
128 
129         // we have converged
130         return current;
131 
132     }
133 
134 }