Gauss–Newton Methods

For minimization problems for which the objective function is a sum of squares,

it is often advantageous to use the special structure of the problem. Time and effort can be saved by computing the residual function r(x), and its derivative, the Jacobian J(x). The Gauss–Newton method is an elegant way to do this. Rather than using the complete second-order Hessian matrix for the quadratic model, the Gauss–Newton method uses (Subscript[B, k]) == ((Subscript[J, k])^T)(Subscript[J, k]) in (1) such that the step (Subscript[p, k]) is computed from the formula

where (Subscript[J, k]) == J((Subscript[x, k])), and so on. Note that this is an approximation to the full Hessian, which is  (J^T)J +(∑_(j==1)^m)(Subscript[r, j])(∇^2)(Subscript[r, j]). In the zero residual case, where r==0 is the minimum, or when r varies nearly as a linear function near the minimum point, the approximation to the Hessian is quite good and the quadratic convergence of Newton's method is commonly observed.

Objective functions, which are sums of squares, are quite common, and, in fact, this is the form of the objective function when (RefLink[FindFit, paclet:ref/FindFit]) is used with the default value of the (RefLink[NormFunction, paclet:ref/NormFunction]) option. One way to view the Gauss–Newton method is in terms of least-squares problems. Solving the Gauss–Newton step is the same as solving a linear least-squares problem, so applying a Gauss–Newton method is in effect applying a sequence of linear least-squares fits to a nonlinear function. With this view, it makes sense that this method is particularly appropriate for the sort of nonlinear fitting that (RefLink[FindFit, paclet:ref/FindFit]) does.

This loads a package that contains some utility functions:

<<Optimization`UnconstrainedProblems`

This uses the Unconstrained Problems Package to set up the classic Rosenbrock function, which has a narrow curved valley:

p = GetFindMinimumProblem[Rosenbrock]

When the Wolfram Language encounters a problem that is expressly a sum of squares, such as the Rosenbrock example, or a function that is the dot product of a vector with itself, the Gauss–Newton method will be used automatically.

This shows the steps taken by (RefLink[FindMinimum, paclet:ref/FindMinimum]) with the Gauss–Newton method for Rosenbrock's function using a trust region method for step control:

FindMinimumPlot[p, Method->Automatic]

If you compare this with the same example done with Newton's method, you can see that it was done with fewer steps and evaluations because the Gauss–Newton method is taking advantage of the special structure of the problem. The convergence rate near the minimum is just as good as for Newton's method because the residual is zero at the minimum.

The Levenberg–Marquardt method is a Gauss–Newton method with trust region step control (though it was originally proposed before the general notion of trust regions had been developed). You can request this method specifically by using the (RefLink[FindMinimum, paclet:ref/FindMinimum]) option (RefLink[Method, paclet:ref/Method])->"LevenbergMarquardt" or equivalently (RefLink[Method, paclet:ref/Method])->"GaussNewton".

Sometimes it is awkward to express a function so that it will explicitly be a sum of squares or a dot product of a vector with itself. In these cases, it is possible to use the "Residual" method option to specify the residual directly. Similarly, you can specify the derivative of the residual with the "Jacobian" method option. Note that when the residual is specified through the "Residual" method option, it is not checked for consistency with the first argument of (RefLink[FindMinimum, paclet:ref/FindMinimum]). The values returned will depend on the value given through the option.

This finds the minimum of Rosenbrock's function using the specification of the residual:

FindMinimum[(1/2)((((1-(Subscript[X, 1])))^2)+100 (((-(Subscript[X, 1]^2)+(Subscript[X, 2])))^2)),{{(Subscript[X, 1]),-1.2`},{(Subscript[X, 2]),1.`}}, Method->{"LevenbergMarquardt", "Residual"->{1-(Subscript[X, 1]),10(-(Subscript[X, 1]^2)+(Subscript[X, 2]))}}]

Method options for (RefLink[Method, paclet:ref/Method])->"LevenbergMarquardt".

Another natural way of setting up sums of squares problems in the Wolfram Language is with (RefLink[FindFit, paclet:ref/FindFit]), which computes nonlinear fits to data. A simple example follows.

Here is a model function:

fm[a_, b_, c_, x_] := a If[x > 0, Cos[b x], Exp[c x]]

Here is some data generated by the function with some random perturbations added:

Block[{ϵ = 0.1, a = 1.2, b = 3.4, c = 0.98}, data = Table[{x, fm[a,b,c, x]+ ϵ RandomReal[ {-.5,.5}]},{x, -5,5,.1}]];

This finds a nonlinear least-squares fit to the model function:

fit = FindFit[data, fm[a,b, c, x], {{a,1}, {b,3}, {c, 1}},x]

This shows the fit model with the data:

Show[{ListPlot[data], Plot[fm[a,b,c, x] /.fit, {x,-5,5}, PlotStyle->RGBColor[0,1,0]]}]

In the example, (RefLink[FindFit, paclet:ref/FindFit]) internally constructs a residual function and Jacobian, which are in turn used by the Gauss–Newton method to find the minimum of the sum of squares, or the nonlinear least-squares fit. Of course, (RefLink[FindFit, paclet:ref/FindFit]) can be used with other methods, but because a residual function that evaluates rapidly can be constructed, it is often faster than the other methods.