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.math4.legacy.fitting.leastsquares;
18
19 import org.apache.commons.math4.legacy.exception.ConvergenceException;
20 import org.apache.commons.math4.legacy.exception.NullArgumentException;
21 import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
22 import org.apache.commons.math4.legacy.fitting.leastsquares.LeastSquaresProblem.Evaluation;
23 import org.apache.commons.math4.legacy.linear.ArrayRealVector;
24 import org.apache.commons.math4.legacy.linear.CholeskyDecomposition;
25 import org.apache.commons.math4.legacy.linear.LUDecomposition;
26 import org.apache.commons.math4.legacy.linear.MatrixUtils;
27 import org.apache.commons.math4.legacy.linear.NonPositiveDefiniteMatrixException;
28 import org.apache.commons.math4.legacy.linear.QRDecomposition;
29 import org.apache.commons.math4.legacy.linear.RealMatrix;
30 import org.apache.commons.math4.legacy.linear.RealVector;
31 import org.apache.commons.math4.legacy.linear.SingularMatrixException;
32 import org.apache.commons.math4.legacy.linear.SingularValueDecomposition;
33 import org.apache.commons.math4.legacy.optim.ConvergenceChecker;
34 import org.apache.commons.math4.legacy.core.IntegerSequence;
35 import org.apache.commons.math4.legacy.core.Pair;
36
37 /**
38 * Gauss-Newton least-squares solver.
39 * <p> This class solve a least-square problem by
40 * solving the normal equations of the linearized problem at each iteration. Either LU
41 * decomposition or Cholesky decomposition can be used to solve the normal equations,
42 * or QR decomposition or SVD decomposition can be used to solve the linear system. LU
43 * decomposition is faster but QR decomposition is more robust for difficult problems,
44 * and SVD can compute a solution for rank-deficient problems.
45 * </p>
46 *
47 * @since 3.3
48 */
49 public class GaussNewtonOptimizer implements LeastSquaresOptimizer {
50
51 /** The decomposition algorithm to use to solve the normal equations. */
52 //TODO move to linear package and expand options?
53 public enum Decomposition {
54 /**
55 * Solve by forming the normal equations (J<sup>T</sup>Jx=J<sup>T</sup>r) and
56 * using the {@link LUDecomposition}.
57 *
58 * <p> Theoretically this method takes mn<sup>2</sup>/2 operations to compute the
59 * normal matrix and n<sup>3</sup>/3 operations (m > n) to solve the system using
60 * the LU decomposition. </p>
61 */
62 LU {
63 @Override
64 protected RealVector solve(final RealMatrix jacobian,
65 final RealVector residuals) {
66 try {
67 final Pair<RealMatrix, RealVector> normalEquation =
68 computeNormalMatrix(jacobian, residuals);
69 final RealMatrix normal = normalEquation.getFirst();
70 final RealVector jTr = normalEquation.getSecond();
71 return new LUDecomposition(normal, SINGULARITY_THRESHOLD)
72 .getSolver()
73 .solve(jTr);
74 } catch (SingularMatrixException e) {
75 throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM, e);
76 }
77 }
78 },
79 /**
80 * Solve the linear least squares problem (Jx=r) using the {@link
81 * QRDecomposition}.
82 *
83 * <p> Theoretically this method takes mn<sup>2</sup> - n<sup>3</sup>/3 operations
84 * (m > n) and has better numerical accuracy than any method that forms the normal
85 * equations. </p>
86 */
87 QR {
88 @Override
89 protected RealVector solve(final RealMatrix jacobian,
90 final RealVector residuals) {
91 try {
92 return new QRDecomposition(jacobian, SINGULARITY_THRESHOLD)
93 .getSolver()
94 .solve(residuals);
95 } catch (SingularMatrixException e) {
96 throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM, e);
97 }
98 }
99 },
100 /**
101 * Solve by forming the normal equations (J<sup>T</sup>Jx=J<sup>T</sup>r) and
102 * using the {@link CholeskyDecomposition}.
103 *
104 * <p> Theoretically this method takes mn<sup>2</sup>/2 operations to compute the
105 * normal matrix and n<sup>3</sup>/6 operations (m > n) to solve the system using
106 * the Cholesky decomposition. </p>
107 */
108 CHOLESKY {
109 @Override
110 protected RealVector solve(final RealMatrix jacobian,
111 final RealVector residuals) {
112 try {
113 final Pair<RealMatrix, RealVector> normalEquation =
114 computeNormalMatrix(jacobian, residuals);
115 final RealMatrix normal = normalEquation.getFirst();
116 final RealVector jTr = normalEquation.getSecond();
117 return new CholeskyDecomposition(
118 normal, SINGULARITY_THRESHOLD, SINGULARITY_THRESHOLD)
119 .getSolver()
120 .solve(jTr);
121 } catch (NonPositiveDefiniteMatrixException e) {
122 throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM, e);
123 }
124 }
125 },
126 /**
127 * Solve the linear least squares problem using the {@link
128 * SingularValueDecomposition}.
129 *
130 * <p> This method is slower, but can provide a solution for rank deficient and
131 * nearly singular systems.
132 */
133 SVD {
134 @Override
135 protected RealVector solve(final RealMatrix jacobian,
136 final RealVector residuals) {
137 return new SingularValueDecomposition(jacobian)
138 .getSolver()
139 .solve(residuals);
140 }
141 };
142
143 /**
144 * Solve the linear least squares problem Jx=r.
145 *
146 * @param jacobian the Jacobian matrix, J. the number of rows >= the number or
147 * columns.
148 * @param residuals the computed residuals, r.
149 * @return the solution x, to the linear least squares problem Jx=r.
150 * @throws ConvergenceException if the matrix properties (e.g. singular) do not
151 * permit a solution.
152 */
153 protected abstract RealVector solve(RealMatrix jacobian,
154 RealVector residuals);
155 }
156
157 /**
158 * The singularity threshold for matrix decompositions. Determines when a {@link
159 * ConvergenceException} is thrown. The current value was the default value for {@link
160 * LUDecomposition}.
161 */
162 private static final double SINGULARITY_THRESHOLD = 1e-11;
163
164 /** Indicator for using LU decomposition. */
165 private final Decomposition decomposition;
166
167 /**
168 * Creates a Gauss Newton optimizer.
169 * <p>
170 * The default for the algorithm is to solve the normal equations using QR
171 * decomposition.
172 */
173 public GaussNewtonOptimizer() {
174 this(Decomposition.QR);
175 }
176
177 /**
178 * Create a Gauss Newton optimizer that uses the given decomposition algorithm to
179 * solve the normal equations.
180 *
181 * @param decomposition the {@link Decomposition} algorithm.
182 */
183 public GaussNewtonOptimizer(final Decomposition decomposition) {
184 this.decomposition = decomposition;
185 }
186
187 /**
188 * Get the matrix decomposition algorithm used to solve the normal equations.
189 *
190 * @return the matrix {@link Decomposition} algoritm.
191 */
192 public Decomposition getDecomposition() {
193 return this.decomposition;
194 }
195
196 /**
197 * Configure the decomposition algorithm.
198 *
199 * @param newDecomposition the {@link Decomposition} algorithm to use.
200 * @return a new instance.
201 */
202 public GaussNewtonOptimizer withDecomposition(final Decomposition newDecomposition) {
203 return new GaussNewtonOptimizer(newDecomposition);
204 }
205
206 /** {@inheritDoc} */
207 @Override
208 public Optimum optimize(final LeastSquaresProblem lsp) {
209 //create local evaluation and iteration counts
210 final IntegerSequence.Incrementor evaluationCounter = lsp.getEvaluationCounter();
211 final IntegerSequence.Incrementor iterationCounter = lsp.getIterationCounter();
212 final ConvergenceChecker<Evaluation> checker
213 = lsp.getConvergenceChecker();
214
215 // Computation will be useless without a checker (see "for-loop").
216 if (checker == null) {
217 throw new NullArgumentException();
218 }
219
220 RealVector currentPoint = lsp.getStart();
221
222 // iterate until convergence is reached
223 Evaluation current = null;
224 while (true) {
225 iterationCounter.increment();
226
227 // evaluate the objective function and its jacobian
228 Evaluation previous = current;
229 // Value of the objective function at "currentPoint".
230 evaluationCounter.increment();
231 current = lsp.evaluate(currentPoint);
232 final RealVector currentResiduals = current.getResiduals();
233 final RealMatrix weightedJacobian = current.getJacobian();
234 currentPoint = current.getPoint();
235
236 // Check convergence.
237 if (previous != null &&
238 checker.converged(iterationCounter.getCount(), previous, current)) {
239 return new OptimumImpl(current,
240 evaluationCounter.getCount(),
241 iterationCounter.getCount());
242 }
243
244 // solve the linearized least squares problem
245 final RealVector dX = this.decomposition.solve(weightedJacobian, currentResiduals);
246 // update the estimated parameters
247 currentPoint = currentPoint.add(dX);
248 }
249 }
250
251 /** {@inheritDoc} */
252 @Override
253 public String toString() {
254 return "GaussNewtonOptimizer{" +
255 "decomposition=" + decomposition +
256 '}';
257 }
258
259 /**
260 * Compute the normal matrix, J<sup>T</sup>J.
261 *
262 * @param jacobian the m by n jacobian matrix, J. Input.
263 * @param residuals the m by 1 residual vector, r. Input.
264 * @return the n by n normal matrix and the n by 1 J<sup>Tr</sup> vector.
265 */
266 private static Pair<RealMatrix, RealVector> computeNormalMatrix(final RealMatrix jacobian,
267 final RealVector residuals) {
268 //since the normal matrix is symmetric, we only need to compute half of it.
269 final int nR = jacobian.getRowDimension();
270 final int nC = jacobian.getColumnDimension();
271 //allocate space for return values
272 final RealMatrix normal = MatrixUtils.createRealMatrix(nC, nC);
273 final RealVector jTr = new ArrayRealVector(nC);
274 //for each measurement
275 for (int i = 0; i < nR; ++i) {
276 //compute JTr for measurement i
277 for (int j = 0; j < nC; j++) {
278 jTr.setEntry(j, jTr.getEntry(j) +
279 residuals.getEntry(i) * jacobian.getEntry(i, j));
280 }
281
282 // add the contribution to the normal matrix for measurement i
283 for (int k = 0; k < nC; ++k) {
284 //only compute the upper triangular part
285 for (int l = k; l < nC; ++l) {
286 normal.setEntry(k, l, normal.getEntry(k, l) +
287 jacobian.getEntry(i, k) * jacobian.getEntry(i, l));
288 }
289 }
290 }
291 //copy the upper triangular part to the lower triangular part.
292 for (int i = 0; i < nC; i++) {
293 for (int j = 0; j < i; j++) {
294 normal.setEntry(i, j, normal.getEntry(j, i));
295 }
296 }
297 return new Pair<>(normal, jTr);
298 }
299 }