Nonlinear least squares

After the solution of square linear systems, we generalized to the case of having more constraints to satisfy than available variables. Our next step is to do the same for nonlinear equations, thus filling out this table:

As in the linear case, we consider only overdetermined problems, where $m>n$. Minimizing a positive quantity is equivalent to minimizing its square, so we could also define the result as minimizing $\phi(\mathbf{x})=\mathbf{f}(\mathbf{x})^T\mathbf{f}(\mathbf{x})$.

4.7.1Gauss–Newton method¶

You should not be surprised to learn that we can formulate an algorithm by substituting a linear model function for $\mathbf{f}$. At a current estimate $\mathbf{x}_k$ we define

where $\mathbf{A}_k$ is the exact $m\times n$ Jacobian matrix, $\mathbf{J}(\mathbf{x}_k)$, or an approximation of it as described in Quasi-Newton methods.

In the square case, we solved $\mathbf{q}=\boldsymbol{0}$ to define the new value for $\mathbf{x}$, leading to the condition $\mathbf{A}_k\mathbf{s}_k=-\mathbf{f}_k$, where $\mathbf{s}_k=\mathbf{x}_{k+1}-\mathbf{x}_k$. Now, with $m>n$, we cannot expect to solve $\mathbf{q}=\boldsymbol{0}$, so instead we define $\mathbf{x}_{k+1}$ as the value that minimizes $\| \mathbf{q} \|_2$.

In brief, Gauss–Newton solves a series of linear least-squares problems in order to solve a nonlinear least-squares problem.

Surprisingly, Function 4.5.2 and Function 4.6.3, which were introduced for the case of $m=n$ nonlinear equations, work without modification as the Gauss–Newton method for the overdetermined case! The reason is that the backslash operator applies equally well to the linear system and linear least-squares problems, and nothing else in those functions was written with explicit reference to $n$.

4.7.2Convergence¶

In the multidimensional Newton method for a nonlinear system, we expect quadratic convergence to a solution in the typical case. For the Gauss–Newton method, the picture is more complicated.

As always in least-squares problems, the residual $\mathbf{f}(\mathbf{x})$ will not necessarily be zero when $\|\mathbf{f}\|$ is minimized. Suppose that the minimum value of $\|\mathbf{f}\|$ is $R>0$. In general, we might observe quadratic-like convergence until the iterate $\|\mathbf{x}_k\|$ is within distance $R$ of a true minimizer, and linear convergence thereafter. When $R$ is not sufficiently small, the convergence can be quite slow.

4.7.3Nonlinear data fitting¶

In Fitting functions to data we saw how to fit functions to data values, provided that the set of candidate fitting functions depends linearly on the undetermined coefficients. We now have a tool to generalize that process to fitting functions that depend nonlinearly on unknown parameters.

Suppose that $(t_i,y_i)$ for $i=1,\ldots,m$ are given points. We wish to model the data by a function $g(t,\mathbf{x})$ that depends on unknown parameters $x_1,\ldots,x_n$ in an arbitrary way. A standard approach is to minimize the discrepancy between the model and the observations, in a least-squares sense. Define

We call $\mathbf{f}$ a misfit function. By minimizing $\bigl\| \mathbf{f}(\mathbf{c}) \bigr\|^2$, we get the best possible fit to the data. If an explicit Jacobian matrix is desired for the minimization, we can compute

The form of $g$ is up to the modeler. There may be compelling theoretical choices, or you may just be looking for enough algebraic power to express the data well. Naturally, in the special case where the dependence on $\mathbf{c}$ is linear, i.e.,

then the misfit function is also linear in $\mathbf{c}$ and the fitting problem reduces to linear least squares.

4.7.4Exercises¶

1. ✍ Define $\mathbf{f}(x)=[ x-8, \; x^2-4 ]$.

   (a) Write out the linear model of $\mathbf{f}$ at $x=2$.

   (b) Find the estimate produced by one step of the Gauss–Newton method, starting at $x=2$.

2. ✍ (Continuation of Exercise 1.) The Gauss–Newton method replaces $\mathbf{f}(\mathbf{x})$ by a linear model and minimizes the norm of its residual. An alternative is to replace $\| \mathbf{f}(\mathbf{x}) \|_2^2$ by a scalar quadratic model $q(\mathbf{x})$ and minimize that.

   (a) Using $\mathbf{f}(x) = [ x-8, \; x^2-4 ]$, let $q(x)$ be defined by the first three terms in the Taylor series for $\| \mathbf{f}(x) \|_2^2$ at $x=2$.

   (b) Find the unique $x$ that minimizes $q(x)$. Is the result the same as the estimate produced by Gauss–Newton?

3. ⌨ A famous result by Kermack and McKendrick in 1927 Kermack et al. (1927) suggests that in epidemics that kill only a small fraction of a susceptible population, the death rate as a function of time is well modeled by

   $$w'(t) = A \operatorname{sech}^2[B(t-C)]$$

   (4.7.6)

   for constant values of the parameters $A,B,C$. Since the maximum of sech is $\operatorname{sech}(0)=1$, $A$ is the maximum death rate and $C$ is the time of peak deaths. You will use this model to fit the deaths per week from plague recorded in Mumbai during 1906:

   5, 10, 17, 22, 30, 50, 51, 90, 120, 180, 292, 395, 445, 775, 780,
   700, 698, 880, 925, 800, 578, 400, 350, 202, 105, 65, 55, 40, 30, 20

   (a) Use Function 4.6.3 to find the best least-squares fit to the data using the $\operatorname{sech}^2$ model. Make a plot of the model fit superimposed on the data.

   (b) Repeat part (a) using only the first 15 data values.

4. ⌨ (Variation on Exercise 4.5.6.) Suppose the points $(x_i,y_i)$ for $i=1,\ldots,m$ are given, and the goal is to find the circle with center $(a,b)$ and radius $r$ that best fits the points. Define

   $$f_i(a,b,r) = (a-x_i)^2 + (b-y_i)^2 - r^2, \qquad i=1,\ldots,m.$$

   (4.7.7)

   Then we can define the best circle as the one that minimizes $\|\mathbf{f}\|$.

   Define data points as follows:

   m = 30; t = 2π*rand(m);
   x = @. -2 + 5*cos(t); y = @. 1 + 5*sin(t);
   x += 0.2*randn(m); y += 0.2*randn(m);

   Use Function 4.6.3 to find the best-fit circle, and make a plot of the circle superimposed on the points.

5. ⌨ The position of the upper lid during an eye blink can be measured from high-speed video Wu et al. (2014), and it may be possible to classify blinks based in part on fits to the lid position Brosch et al. (2017). The lid position functions proposed to fit blinks is a product of a monomial or polynomial multiplying a decaying exponential Berke & Mueller (1998). In this problem, you will generate representative data, add a small amount of noise to it, and then perform nonlinear least-squares fits to the data.

   (a) Consider the function $y(\mathbf{a}) = a_1 t^2 \exp \left( -a_2 t^{a_3} \right)$, using the vector of coefficients $\mathbf{a} = [a_1,a_2,a_3]$, and create synthetic eyelid position data as follows:

   N = 20;                            # number of time values
   t = (1:N)/N;                       # equally spaced to t=1
   a = [10, 10, 2];                   # baseline values
   y = @. a(1)*t^2*exp(-a(2)*t^a(3)); # ideal data
   ym = copy(y);                      # vector for data
   ir = 1:N-1;                        # range to add noise
   noise = 0.03;                      # amplitude of noise
   ym[ir] += noise*rand(N-1);         # add noise

   (b) Using the data (t,ym), find the nonlinear least-squares fit using Function 4.6.3.

   (c) Plot the fits using np = 100 points over t=(1:np)/np together with symbols for the N measured data points ym.

   (d) Increase the noise to 5% and 10%. You may have to increase the number of measured points N and/or the maximum number of iterations. How close are the coefficients? Plot the data and the resulting fit for each case.

6. ⌨ Repeat the previous problem using the fitting function $y(\mathbf{a}) = (a_1+a_2 t + a_3 t^2) t^2 \exp \left( -a_4 t^{a_5} \right)$, using the vector of coefficients $\mathbf{a} = [a_1,\ldots,a_5]$. (This was the choice used in Brosch et al Brosch et al. (2017).) Use a = [20, -10, -8, 7, 2] to create the data and as an initial guess for the coefficients for the fit to the noisy data.