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  package org.apache.commons.math3.optim.nonlinear.vector.jacobian;
18  
19  import org.apache.commons.math3.exception.ConvergenceException;
20  import org.apache.commons.math3.exception.NullArgumentException;
21  import org.apache.commons.math3.exception.MathInternalError;
22  import org.apache.commons.math3.exception.MathUnsupportedOperationException;
23  import org.apache.commons.math3.exception.util.LocalizedFormats;
24  import org.apache.commons.math3.linear.ArrayRealVector;
25  import org.apache.commons.math3.linear.BlockRealMatrix;
26  import org.apache.commons.math3.linear.DecompositionSolver;
27  import org.apache.commons.math3.linear.LUDecomposition;
28  import org.apache.commons.math3.linear.QRDecomposition;
29  import org.apache.commons.math3.linear.RealMatrix;
30  import org.apache.commons.math3.linear.SingularMatrixException;
31  import org.apache.commons.math3.optim.ConvergenceChecker;
32  import org.apache.commons.math3.optim.PointVectorValuePair;
33  
34  /**
35   * Gauss-Newton least-squares solver.
36   * <br/>
37   * Constraints are not supported: the call to
38   * {@link #optimize(OptimizationData[]) optimize} will throw
39   * {@link MathUnsupportedOperationException} if bounds are passed to it.
40   *
41   * <p>
42   * This class solve a least-square problem by solving the normal equations
43   * of the linearized problem at each iteration. Either LU decomposition or
44   * QR decomposition can be used to solve the normal equations. LU decomposition
45   * is faster but QR decomposition is more robust for difficult problems.
46   * </p>
47   *
48   * @version $Id: GaussNewtonOptimizer.java 1515242 2013-08-18 23:27:29Z erans $
49   * @since 2.0
50   * @deprecated All classes and interfaces in this package are deprecated.
51   * The optimizers that were provided here were moved to the
52   * {@link org.apache.commons.math3.fitting.leastsquares} package
53   * (cf. MATH-1008).
54   */
55  @Deprecated
56  public class GaussNewtonOptimizer extends AbstractLeastSquaresOptimizer {
57      /** Indicator for using LU decomposition. */
58      private final boolean useLU;
59  
60      /**
61       * Simple constructor with default settings.
62       * The normal equations will be solved using LU decomposition.
63       *
64       * @param checker Convergence checker.
65       */
66      public GaussNewtonOptimizer(ConvergenceChecker<PointVectorValuePair> checker) {
67          this(true, checker);
68      }
69  
70      /**
71       * @param useLU If {@code true}, the normal equations will be solved
72       * using LU decomposition, otherwise they will be solved using QR
73       * decomposition.
74       * @param checker Convergence checker.
75       */
76      public GaussNewtonOptimizer(final boolean useLU,
77                                  ConvergenceChecker<PointVectorValuePair> checker) {
78          super(checker);
79          this.useLU = useLU;
80      }
81  
82      /** {@inheritDoc} */
83      @Override
84      public PointVectorValuePair doOptimize() {
85          checkParameters();
86  
87          final ConvergenceChecker<PointVectorValuePair> checker
88              = getConvergenceChecker();
89  
90          // Computation will be useless without a checker (see "for-loop").
91          if (checker == null) {
92              throw new NullArgumentException();
93          }
94  
95          final double[] targetValues = getTarget();
96          final int nR = targetValues.length; // Number of observed data.
97  
98          final RealMatrix weightMatrix = getWeight();
99          // Diagonal of the weight matrix.
100         final double[] residualsWeights = new double[nR];
101         for (int i = 0; i < nR; i++) {
102             residualsWeights[i] = weightMatrix.getEntry(i, i);
103         }
104 
105         final double[] currentPoint = getStartPoint();
106         final int nC = currentPoint.length;
107 
108         // iterate until convergence is reached
109         PointVectorValuePair current = null;
110         for (boolean converged = false; !converged;) {
111             incrementIterationCount();
112 
113             // evaluate the objective function and its jacobian
114             PointVectorValuePair previous = current;
115             // Value of the objective function at "currentPoint".
116             final double[] currentObjective = computeObjectiveValue(currentPoint);
117             final double[] currentResiduals = computeResiduals(currentObjective);
118             final RealMatrix weightedJacobian = computeWeightedJacobian(currentPoint);
119             current = new PointVectorValuePair(currentPoint, currentObjective);
120 
121             // build the linear problem
122             final double[]   b = new double[nC];
123             final double[][] a = new double[nC][nC];
124             for (int i = 0; i < nR; ++i) {
125 
126                 final double[] grad   = weightedJacobian.getRow(i);
127                 final double weight   = residualsWeights[i];
128                 final double residual = currentResiduals[i];
129 
130                 // compute the normal equation
131                 final double wr = weight * residual;
132                 for (int j = 0; j < nC; ++j) {
133                     b[j] += wr * grad[j];
134                 }
135 
136                 // build the contribution matrix for measurement i
137                 for (int k = 0; k < nC; ++k) {
138                     double[] ak = a[k];
139                     double wgk = weight * grad[k];
140                     for (int l = 0; l < nC; ++l) {
141                         ak[l] += wgk * grad[l];
142                     }
143                 }
144             }
145 
146             // Check convergence.
147             if (previous != null) {
148                 converged = checker.converged(getIterations(), previous, current);
149                 if (converged) {
150                     setCost(computeCost(currentResiduals));
151                     return current;
152                 }
153             }
154 
155             try {
156                 // solve the linearized least squares problem
157                 RealMatrix mA = new BlockRealMatrix(a);
158                 DecompositionSolver solver = useLU ?
159                         new LUDecomposition(mA).getSolver() :
160                         new QRDecomposition(mA).getSolver();
161                 final double[] dX = solver.solve(new ArrayRealVector(b, false)).toArray();
162                 // update the estimated parameters
163                 for (int i = 0; i < nC; ++i) {
164                     currentPoint[i] += dX[i];
165                 }
166             } catch (SingularMatrixException e) {
167                 throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM);
168             }
169         }
170         // Must never happen.
171         throw new MathInternalError();
172     }
173 
174     /**
175      * @throws MathUnsupportedOperationException if bounds were passed to the
176      * {@link #optimize(OptimizationData[]) optimize} method.
177      */
178     private void checkParameters() {
179         if (getLowerBound() != null ||
180             getUpperBound() != null) {
181             throw new MathUnsupportedOperationException(LocalizedFormats.CONSTRAINT);
182         }
183     }
184 }