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