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.exception.util.LocalizedFormats; 22 import org.apache.commons.math.linear.BlockRealMatrix; 23 import org.apache.commons.math.linear.DecompositionSolver; 24 import org.apache.commons.math.linear.InvalidMatrixException; 25 import org.apache.commons.math.linear.LUDecompositionImpl; 26 import org.apache.commons.math.linear.QRDecompositionImpl; 27 import org.apache.commons.math.linear.RealMatrix; 28 import org.apache.commons.math.optimization.OptimizationException; 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: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 fvr. 2011) $ 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 52 * org.apache.commons.math.optimization.SimpleVectorialValueChecker} 53 * and the maximal number of evaluation is set to 54 * {@link AbstractLeastSquaresOptimizer#DEFAULT_MAX_ITERATIONS}. 55 * @param useLU if true, the normal equations will be solved using LU 56 * decomposition, otherwise they will be solved using QR decomposition 57 */ 58 public GaussNewtonOptimizer(final boolean useLU) { 59 this.useLU = useLU; 60 } 61 62 /** {@inheritDoc} */ 63 @Override 64 public VectorialPointValuePair doOptimize() 65 throws FunctionEvaluationException, OptimizationException, IllegalArgumentException { 66 67 // iterate until convergence is reached 68 VectorialPointValuePair current = null; 69 for (boolean converged = false; !converged;) { 70 71 incrementIterationsCounter(); 72 73 // evaluate the objective function and its jacobian 74 VectorialPointValuePair previous = current; 75 updateResidualsAndCost(); 76 updateJacobian(); 77 current = new VectorialPointValuePair(point, objective); 78 79 // build the linear problem 80 final double[] b = new double[cols]; 81 final double[][] a = new double[cols][cols]; 82 for (int i = 0; i < rows; ++i) { 83 84 final double[] grad = jacobian[i]; 85 final double weight = residualsWeights[i]; 86 final double residual = objective[i] - targetValues[i]; 87 88 // compute the normal equation 89 final double wr = weight * residual; 90 for (int j = 0; j < cols; ++j) { 91 b[j] += wr * grad[j]; 92 } 93 94 // build the contribution matrix for measurement i 95 for (int k = 0; k < cols; ++k) { 96 double[] ak = a[k]; 97 double wgk = weight * grad[k]; 98 for (int l = 0; l < cols; ++l) { 99 ak[l] += wgk * grad[l]; 100 } 101 } 102 103 } 104 105 try { 106 107 // solve the linearized least squares problem 108 RealMatrix mA = new BlockRealMatrix(a); 109 DecompositionSolver solver = useLU ? 110 new LUDecompositionImpl(mA).getSolver() : 111 new QRDecompositionImpl(mA).getSolver(); 112 final double[] dX = solver.solve(b); 113 114 // update the estimated parameters 115 for (int i = 0; i < cols; ++i) { 116 point[i] += dX[i]; 117 } 118 119 } catch(InvalidMatrixException e) { 120 throw new OptimizationException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM); 121 } 122 123 // check convergence 124 if (previous != null) { 125 converged = checker.converged(getIterations(), previous, current); 126 } 127 128 } 129 130 // we have converged 131 return current; 132 133 } 134 135 } 136
