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.math.optimization.general;
19
20 import org.apache.commons.math.FunctionEvaluationException;
21 import org.apache.commons.math.linear.BlockRealMatrix;
22 import org.apache.commons.math.linear.DecompositionSolver;
23 import org.apache.commons.math.linear.InvalidMatrixException;
24 import org.apache.commons.math.linear.LUDecompositionImpl;
25 import org.apache.commons.math.linear.QRDecompositionImpl;
26 import org.apache.commons.math.linear.RealMatrix;
27 import org.apache.commons.math.optimization.OptimizationException;
28 import org.apache.commons.math.optimization.SimpleVectorialValueChecker;
29 import org.apache.commons.math.optimization.VectorialPointValuePair;
30
31 /**
32 * Gauss-Newton least-squares solver.
33 * <p>
34 * This class solve a least-square problem by solving the normal equations
35 * of the linearized problem at each iteration. Either LU decomposition or
36 * QR decomposition can be used to solve the normal equations. LU decomposition
37 * is faster but QR decomposition is more robust for difficult problems.
38 * </p>
39 *
40 * @version $Revision: 795978 $ $Date: 2009-07-20 15:57:08 -0400 (Mon, 20 Jul 2009) $
41 * @since 2.0
42 *
43 */
44
45 public class GaussNewtonOptimizer extends AbstractLeastSquaresOptimizer {
46
47 /** Indicator for using LU decomposition. */
48 private final boolean useLU;
49
50 /** Simple constructor with default settings.
51 * <p>The convergence check is set to a {@link SimpleVectorialValueChecker}
52 * and the maximal number of evaluation is set to
53 * {@link AbstractLeastSquaresOptimizer#DEFAULT_MAX_ITERATIONS}.
54 * @param useLU if true, the normal equations will be solved using LU
55 * decomposition, otherwise they will be solved using QR decomposition
56 */
57 public GaussNewtonOptimizer(final boolean useLU) {
58 this.useLU = useLU;
59 }
60
61 /** {@inheritDoc} */
62 @Override
63 public VectorialPointValuePair doOptimize()
64 throws FunctionEvaluationException, OptimizationException, IllegalArgumentException {
65
66 // iterate until convergence is reached
67 VectorialPointValuePair current = null;
68 for (boolean converged = false; !converged;) {
69
70 incrementIterationsCounter();
71
72 // evaluate the objective function and its jacobian
73 VectorialPointValuePair previous = current;
74 updateResidualsAndCost();
75 updateJacobian();
76 current = new VectorialPointValuePair(point, objective);
77
78 // build the linear problem
79 final double[] b = new double[cols];
80 final double[][] a = new double[cols][cols];
81 for (int i = 0; i < rows; ++i) {
82
83 final double[] grad = jacobian[i];
84 final double weight = weights[i];
85 final double residual = objective[i] - target[i];
86
87 // compute the normal equation
88 final double wr = weight * residual;
89 for (int j = 0; j < cols; ++j) {
90 b[j] += wr * grad[j];
91 }
92
93 // build the contribution matrix for measurement i
94 for (int k = 0; k < cols; ++k) {
95 double[] ak = a[k];
96 double wgk = weight * grad[k];
97 for (int l = 0; l < cols; ++l) {
98 ak[l] += wgk * grad[l];
99 }
100 }
101
102 }
103
104 try {
105
106 // solve the linearized least squares problem
107 RealMatrix mA = new BlockRealMatrix(a);
108 DecompositionSolver solver = useLU ?
109 new LUDecompositionImpl(mA).getSolver() :
110 new QRDecompositionImpl(mA).getSolver();
111 final double[] dX = solver.solve(b);
112
113 // update the estimated parameters
114 for (int i = 0; i < cols; ++i) {
115 point[i] += dX[i];
116 }
117
118 } catch(InvalidMatrixException e) {
119 throw new OptimizationException("unable to solve: singular problem");
120 }
121
122 // check convergence
123 if (previous != null) {
124 converged = checker.converged(getIterations(), previous, current);
125 }
126
127 }
128
129 // we have converged
130 return current;
131
132 }
133
134 }