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.exception.util.LocalizedFormats;
022 import org.apache.commons.math.linear.BlockRealMatrix;
023 import org.apache.commons.math.linear.DecompositionSolver;
024 import org.apache.commons.math.linear.InvalidMatrixException;
025 import org.apache.commons.math.linear.LUDecompositionImpl;
026 import org.apache.commons.math.linear.QRDecompositionImpl;
027 import org.apache.commons.math.linear.RealMatrix;
028 import org.apache.commons.math.optimization.OptimizationException;
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: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 f??vr. 2011) $
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
052 * org.apache.commons.math.optimization.SimpleVectorialValueChecker}
053 * and the maximal number of evaluation is set to
054 * {@link AbstractLeastSquaresOptimizer#DEFAULT_MAX_ITERATIONS}.
055 * @param useLU if true, the normal equations will be solved using LU
056 * decomposition, otherwise they will be solved using QR decomposition
057 */
058 public GaussNewtonOptimizer(final boolean useLU) {
059 this.useLU = useLU;
060 }
061
062 /** {@inheritDoc} */
063 @Override
064 public VectorialPointValuePair doOptimize()
065 throws FunctionEvaluationException, OptimizationException, IllegalArgumentException {
066
067 // iterate until convergence is reached
068 VectorialPointValuePair current = null;
069 for (boolean converged = false; !converged;) {
070
071 incrementIterationsCounter();
072
073 // evaluate the objective function and its jacobian
074 VectorialPointValuePair previous = current;
075 updateResidualsAndCost();
076 updateJacobian();
077 current = new VectorialPointValuePair(point, objective);
078
079 // build the linear problem
080 final double[] b = new double[cols];
081 final double[][] a = new double[cols][cols];
082 for (int i = 0; i < rows; ++i) {
083
084 final double[] grad = jacobian[i];
085 final double weight = residualsWeights[i];
086 final double residual = objective[i] - targetValues[i];
087
088 // compute the normal equation
089 final double wr = weight * residual;
090 for (int j = 0; j < cols; ++j) {
091 b[j] += wr * grad[j];
092 }
093
094 // build the contribution matrix for measurement i
095 for (int k = 0; k < cols; ++k) {
096 double[] ak = a[k];
097 double wgk = weight * grad[k];
098 for (int l = 0; l < cols; ++l) {
099 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 }