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 */ 017package org.apache.commons.math3.optim.nonlinear.vector.jacobian; 018 019import org.apache.commons.math3.exception.ConvergenceException; 020import org.apache.commons.math3.exception.NullArgumentException; 021import org.apache.commons.math3.exception.MathInternalError; 022import org.apache.commons.math3.exception.MathUnsupportedOperationException; 023import org.apache.commons.math3.exception.util.LocalizedFormats; 024import org.apache.commons.math3.linear.ArrayRealVector; 025import org.apache.commons.math3.linear.BlockRealMatrix; 026import org.apache.commons.math3.linear.DecompositionSolver; 027import org.apache.commons.math3.linear.LUDecomposition; 028import org.apache.commons.math3.linear.QRDecomposition; 029import org.apache.commons.math3.linear.RealMatrix; 030import org.apache.commons.math3.linear.SingularMatrixException; 031import org.apache.commons.math3.optim.ConvergenceChecker; 032import org.apache.commons.math3.optim.PointVectorValuePair; 033 034/** 035 * Gauss-Newton least-squares solver. 036 * <br/> 037 * Constraints are not supported: the call to 038 * {@link #optimize(OptimizationData[]) optimize} will throw 039 * {@link MathUnsupportedOperationException} if bounds are passed to it. 040 * 041 * <p> 042 * This class solve a least-square problem by solving the normal equations 043 * of the linearized problem at each iteration. Either LU decomposition or 044 * QR decomposition can be used to solve the normal equations. LU decomposition 045 * is faster but QR decomposition is more robust for difficult problems. 046 * </p> 047 * 048 * @since 2.0 049 * @deprecated All classes and interfaces in this package are deprecated. 050 * The optimizers that were provided here were moved to the 051 * {@link org.apache.commons.math3.fitting.leastsquares} package 052 * (cf. MATH-1008). 053 */ 054@Deprecated 055public class GaussNewtonOptimizer extends AbstractLeastSquaresOptimizer { 056 /** Indicator for using LU decomposition. */ 057 private final boolean useLU; 058 059 /** 060 * Simple constructor with default settings. 061 * The normal equations will be solved using LU decomposition. 062 * 063 * @param checker Convergence checker. 064 */ 065 public GaussNewtonOptimizer(ConvergenceChecker<PointVectorValuePair> checker) { 066 this(true, checker); 067 } 068 069 /** 070 * @param useLU If {@code true}, the normal equations will be solved 071 * using LU decomposition, otherwise they will be solved using QR 072 * decomposition. 073 * @param checker Convergence checker. 074 */ 075 public GaussNewtonOptimizer(final boolean useLU, 076 ConvergenceChecker<PointVectorValuePair> checker) { 077 super(checker); 078 this.useLU = useLU; 079 } 080 081 /** {@inheritDoc} */ 082 @Override 083 public PointVectorValuePair doOptimize() { 084 checkParameters(); 085 086 final ConvergenceChecker<PointVectorValuePair> checker 087 = getConvergenceChecker(); 088 089 // Computation will be useless without a checker (see "for-loop"). 090 if (checker == null) { 091 throw new NullArgumentException(); 092 } 093 094 final double[] targetValues = getTarget(); 095 final int nR = targetValues.length; // Number of observed data. 096 097 final RealMatrix weightMatrix = getWeight(); 098 // Diagonal of the weight matrix. 099 final double[] residualsWeights = new double[nR]; 100 for (int i = 0; i < nR; i++) { 101 residualsWeights[i] = weightMatrix.getEntry(i, i); 102 } 103 104 final double[] currentPoint = getStartPoint(); 105 final int nC = currentPoint.length; 106 107 // iterate until convergence is reached 108 PointVectorValuePair current = null; 109 for (boolean converged = false; !converged;) { 110 incrementIterationCount(); 111 112 // evaluate the objective function and its jacobian 113 PointVectorValuePair previous = current; 114 // Value of the objective function at "currentPoint". 115 final double[] currentObjective = computeObjectiveValue(currentPoint); 116 final double[] currentResiduals = computeResiduals(currentObjective); 117 final RealMatrix weightedJacobian = computeWeightedJacobian(currentPoint); 118 current = new PointVectorValuePair(currentPoint, currentObjective); 119 120 // build the linear problem 121 final double[] b = new double[nC]; 122 final double[][] a = new double[nC][nC]; 123 for (int i = 0; i < nR; ++i) { 124 125 final double[] grad = weightedJacobian.getRow(i); 126 final double weight = residualsWeights[i]; 127 final double residual = currentResiduals[i]; 128 129 // compute the normal equation 130 final double wr = weight * residual; 131 for (int j = 0; j < nC; ++j) { 132 b[j] += wr * grad[j]; 133 } 134 135 // build the contribution matrix for measurement i 136 for (int k = 0; k < nC; ++k) { 137 double[] ak = a[k]; 138 double wgk = weight * grad[k]; 139 for (int l = 0; l < nC; ++l) { 140 ak[l] += wgk * grad[l]; 141 } 142 } 143 } 144 145 // Check convergence. 146 if (previous != null) { 147 converged = checker.converged(getIterations(), previous, current); 148 if (converged) { 149 setCost(computeCost(currentResiduals)); 150 return current; 151 } 152 } 153 154 try { 155 // solve the linearized least squares problem 156 RealMatrix mA = new BlockRealMatrix(a); 157 DecompositionSolver solver = useLU ? 158 new LUDecomposition(mA).getSolver() : 159 new QRDecomposition(mA).getSolver(); 160 final double[] dX = solver.solve(new ArrayRealVector(b, false)).toArray(); 161 // update the estimated parameters 162 for (int i = 0; i < nC; ++i) { 163 currentPoint[i] += dX[i]; 164 } 165 } catch (SingularMatrixException e) { 166 throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM); 167 } 168 } 169 // Must never happen. 170 throw new MathInternalError(); 171 } 172 173 /** 174 * @throws MathUnsupportedOperationException if bounds were passed to the 175 * {@link #optimize(OptimizationData[]) optimize} method. 176 */ 177 private void checkParameters() { 178 if (getLowerBound() != null || 179 getUpperBound() != null) { 180 throw new MathUnsupportedOperationException(LocalizedFormats.CONSTRAINT); 181 } 182 } 183}