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.VectorialPointValuePair;
029
030 /**
031 * Gauss-Newton least-squares solver.
032 * <p>
033 * This class solve a least-square problem by solving the normal equations
034 * of the linearized problem at each iteration. Either LU decomposition or
035 * QR decomposition can be used to solve the normal equations. LU decomposition
036 * is faster but QR decomposition is more robust for difficult problems.
037 * </p>
038 *
039 * @version $Revision: 811685 $ $Date: 2009-09-05 13:36:48 -0400 (Sat, 05 Sep 2009) $
040 * @since 2.0
041 *
042 */
043
044 public class GaussNewtonOptimizer extends AbstractLeastSquaresOptimizer {
045
046 /** Indicator for using LU decomposition. */
047 private final boolean useLU;
048
049 /** Simple constructor with default settings.
050 * <p>The convergence check is set to a {@link
051 * org.apache.commons.math.optimization.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 = residualsWeights[i];
085 final double residual = objective[i] - targetValues[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 }