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 agents.anac.y2019.harddealer.math3.fitting.leastsquares;
|
---|
18 |
|
---|
19 | import agents.anac.y2019.harddealer.math3.exception.ConvergenceException;
|
---|
20 | import agents.anac.y2019.harddealer.math3.exception.NullArgumentException;
|
---|
21 | import agents.anac.y2019.harddealer.math3.exception.util.LocalizedFormats;
|
---|
22 | import agents.anac.y2019.harddealer.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation;
|
---|
23 | import agents.anac.y2019.harddealer.math3.linear.ArrayRealVector;
|
---|
24 | import agents.anac.y2019.harddealer.math3.linear.CholeskyDecomposition;
|
---|
25 | import agents.anac.y2019.harddealer.math3.linear.LUDecomposition;
|
---|
26 | import agents.anac.y2019.harddealer.math3.linear.MatrixUtils;
|
---|
27 | import agents.anac.y2019.harddealer.math3.linear.NonPositiveDefiniteMatrixException;
|
---|
28 | import agents.anac.y2019.harddealer.math3.linear.QRDecomposition;
|
---|
29 | import agents.anac.y2019.harddealer.math3.linear.RealMatrix;
|
---|
30 | import agents.anac.y2019.harddealer.math3.linear.RealVector;
|
---|
31 | import agents.anac.y2019.harddealer.math3.linear.SingularMatrixException;
|
---|
32 | import agents.anac.y2019.harddealer.math3.linear.SingularValueDecomposition;
|
---|
33 | import agents.anac.y2019.harddealer.math3.optim.ConvergenceChecker;
|
---|
34 | import agents.anac.y2019.harddealer.math3.util.Incrementor;
|
---|
35 | import agents.anac.y2019.harddealer.math3.util.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 | public Optimum optimize(final LeastSquaresProblem lsp) {
|
---|
208 | //create local evaluation and iteration counts
|
---|
209 | final Incrementor evaluationCounter = lsp.getEvaluationCounter();
|
---|
210 | final Incrementor iterationCounter = lsp.getIterationCounter();
|
---|
211 | final ConvergenceChecker<Evaluation> checker
|
---|
212 | = lsp.getConvergenceChecker();
|
---|
213 |
|
---|
214 | // Computation will be useless without a checker (see "for-loop").
|
---|
215 | if (checker == null) {
|
---|
216 | throw new NullArgumentException();
|
---|
217 | }
|
---|
218 |
|
---|
219 | RealVector currentPoint = lsp.getStart();
|
---|
220 |
|
---|
221 | // iterate until convergence is reached
|
---|
222 | Evaluation current = null;
|
---|
223 | while (true) {
|
---|
224 | iterationCounter.incrementCount();
|
---|
225 |
|
---|
226 | // evaluate the objective function and its jacobian
|
---|
227 | Evaluation previous = current;
|
---|
228 | // Value of the objective function at "currentPoint".
|
---|
229 | evaluationCounter.incrementCount();
|
---|
230 | current = lsp.evaluate(currentPoint);
|
---|
231 | final RealVector currentResiduals = current.getResiduals();
|
---|
232 | final RealMatrix weightedJacobian = current.getJacobian();
|
---|
233 | currentPoint = current.getPoint();
|
---|
234 |
|
---|
235 | // Check convergence.
|
---|
236 | if (previous != null &&
|
---|
237 | checker.converged(iterationCounter.getCount(), previous, current)) {
|
---|
238 | return new OptimumImpl(current,
|
---|
239 | evaluationCounter.getCount(),
|
---|
240 | iterationCounter.getCount());
|
---|
241 | }
|
---|
242 |
|
---|
243 | // solve the linearized least squares problem
|
---|
244 | final RealVector dX = this.decomposition.solve(weightedJacobian, currentResiduals);
|
---|
245 | // update the estimated parameters
|
---|
246 | currentPoint = currentPoint.add(dX);
|
---|
247 | }
|
---|
248 | }
|
---|
249 |
|
---|
250 | /** {@inheritDoc} */
|
---|
251 | @Override
|
---|
252 | public String toString() {
|
---|
253 | return "GaussNewtonOptimizer{" +
|
---|
254 | "decomposition=" + decomposition +
|
---|
255 | '}';
|
---|
256 | }
|
---|
257 |
|
---|
258 | /**
|
---|
259 | * Compute the normal matrix, J<sup>T</sup>J.
|
---|
260 | *
|
---|
261 | * @param jacobian the m by n jacobian matrix, J. Input.
|
---|
262 | * @param residuals the m by 1 residual vector, r. Input.
|
---|
263 | * @return the n by n normal matrix and the n by 1 J<sup>Tr vector.
|
---|
264 | */
|
---|
265 | private static Pair<RealMatrix, RealVector> computeNormalMatrix(final RealMatrix jacobian,
|
---|
266 | final RealVector residuals) {
|
---|
267 | //since the normal matrix is symmetric, we only need to compute half of it.
|
---|
268 | final int nR = jacobian.getRowDimension();
|
---|
269 | final int nC = jacobian.getColumnDimension();
|
---|
270 | //allocate space for return values
|
---|
271 | final RealMatrix normal = MatrixUtils.createRealMatrix(nC, nC);
|
---|
272 | final RealVector jTr = new ArrayRealVector(nC);
|
---|
273 | //for each measurement
|
---|
274 | for (int i = 0; i < nR; ++i) {
|
---|
275 | //compute JTr for measurement i
|
---|
276 | for (int j = 0; j < nC; j++) {
|
---|
277 | jTr.setEntry(j, jTr.getEntry(j) +
|
---|
278 | residuals.getEntry(i) * jacobian.getEntry(i, j));
|
---|
279 | }
|
---|
280 |
|
---|
281 | // add the the contribution to the normal matrix for measurement i
|
---|
282 | for (int k = 0; k < nC; ++k) {
|
---|
283 | //only compute the upper triangular part
|
---|
284 | for (int l = k; l < nC; ++l) {
|
---|
285 | normal.setEntry(k, l, normal.getEntry(k, l) +
|
---|
286 | jacobian.getEntry(i, k) * jacobian.getEntry(i, l));
|
---|
287 | }
|
---|
288 | }
|
---|
289 | }
|
---|
290 | //copy the upper triangular part to the lower triangular part.
|
---|
291 | for (int i = 0; i < nC; i++) {
|
---|
292 | for (int j = 0; j < i; j++) {
|
---|
293 | normal.setEntry(i, j, normal.getEntry(j, i));
|
---|
294 | }
|
---|
295 | }
|
---|
296 | return new Pair<RealMatrix, RealVector>(normal, jTr);
|
---|
297 | }
|
---|
298 |
|
---|
299 | }
|
---|