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