001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018 package org.apache.commons.math.optimization.general; 019 020 import org.apache.commons.math.FunctionEvaluationException; 021 import org.apache.commons.math.linear.BlockRealMatrix; 022 import org.apache.commons.math.linear.DecompositionSolver; 023 import org.apache.commons.math.linear.InvalidMatrixException; 024 import org.apache.commons.math.linear.LUDecompositionImpl; 025 import org.apache.commons.math.linear.QRDecompositionImpl; 026 import org.apache.commons.math.linear.RealMatrix; 027 import org.apache.commons.math.optimization.OptimizationException; 028 import org.apache.commons.math.optimization.SimpleVectorialValueChecker; 029 import org.apache.commons.math.optimization.VectorialPointValuePair; 030 031 /** 032 * Gauss-Newton least-squares solver. 033 * <p> 034 * This class solve a least-square problem by solving the normal equations 035 * of the linearized problem at each iteration. Either LU decomposition or 036 * QR decomposition can be used to solve the normal equations. LU decomposition 037 * is faster but QR decomposition is more robust for difficult problems. 038 * </p> 039 * 040 * @version $Revision: 795978 $ $Date: 2009-07-20 15:57:08 -0400 (Mon, 20 Jul 2009) $ 041 * @since 2.0 042 * 043 */ 044 045 public class GaussNewtonOptimizer extends AbstractLeastSquaresOptimizer { 046 047 /** Indicator for using LU decomposition. */ 048 private final boolean useLU; 049 050 /** Simple constructor with default settings. 051 * <p>The convergence check is set to a {@link SimpleVectorialValueChecker} 052 * and the maximal number of evaluation is set to 053 * {@link AbstractLeastSquaresOptimizer#DEFAULT_MAX_ITERATIONS}. 054 * @param useLU if true, the normal equations will be solved using LU 055 * decomposition, otherwise they will be solved using QR decomposition 056 */ 057 public GaussNewtonOptimizer(final boolean useLU) { 058 this.useLU = useLU; 059 } 060 061 /** {@inheritDoc} */ 062 @Override 063 public VectorialPointValuePair doOptimize() 064 throws FunctionEvaluationException, OptimizationException, IllegalArgumentException { 065 066 // iterate until convergence is reached 067 VectorialPointValuePair current = null; 068 for (boolean converged = false; !converged;) { 069 070 incrementIterationsCounter(); 071 072 // evaluate the objective function and its jacobian 073 VectorialPointValuePair previous = current; 074 updateResidualsAndCost(); 075 updateJacobian(); 076 current = new VectorialPointValuePair(point, objective); 077 078 // build the linear problem 079 final double[] b = new double[cols]; 080 final double[][] a = new double[cols][cols]; 081 for (int i = 0; i < rows; ++i) { 082 083 final double[] grad = jacobian[i]; 084 final double weight = weights[i]; 085 final double residual = objective[i] - target[i]; 086 087 // compute the normal equation 088 final double wr = weight * residual; 089 for (int j = 0; j < cols; ++j) { 090 b[j] += wr * grad[j]; 091 } 092 093 // build the contribution matrix for measurement i 094 for (int k = 0; k < cols; ++k) { 095 double[] ak = a[k]; 096 double wgk = weight * grad[k]; 097 for (int l = 0; l < cols; ++l) { 098 ak[l] += wgk * grad[l]; 099 } 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 }