Nonlinear least squares - Fundamentals of Numerical Computation

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:

	linear	nonlinear
square	$\mathbf{A}\mathbf{x}=\mathbf{b}$	$\mathbf{f}(\mathbf{x})=\boldsymbol{0}$
overdetermined	$\min\, \bigl\|\mathbf{A}\mathbf{x} - \mathbf{b}\bigr\|_2$	$\min\, \bigl\|\mathbf{f}(\mathbf{x}) \bigr\|_2$

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

\mathbf{q}(\mathbf{x}) = \mathbf{f}(\mathbf{x}_k) + \mathbf{A}_k(\mathbf{x}-\mathbf{x}_k),

(4.7.1)

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.

Example 4.7.1 (Convergence of nonlinear least squares)

Julia

MATLAB

Python

Example 4.7.1

We will observe the convergence of Function 4.6.3 for different levels of the minimum least-squares residual. We start with a function mapping from $\real^2$ into $\real^3$ , and a point that will be near the optimum.

g(x) = [sin(x[1] + x[2]), cos(x[1] - x[2]), exp(x[1] - x[2])]
p = [1, 1];

The function $\mathbf{g}(\mathbf{x}) - \mathbf{g}(\mathbf{p})$ obviously has a zero residual at $\mathbf{p}$ . We’ll make different perturbations of that function in order to create nonzero residuals.

using Printf
plt = plot(xlabel="iteration", yaxis=(:log10, "error"),
    title="Convergence of Gauss–Newton")
for R in [1e-3, 1e-2, 1e-1]
    # Define the perturbed function.
    f(x) = g(x) - g(p) + R * normalize([-1, 1, -1])
    x = FNC.levenberg(f, [0, 0])
    r = x[end]
    err = [norm(x - r) for x in x[1:end-1]]
    normres = norm(f(r))
    plot!(err, label=@sprintf("R=%.2g", normres))
end
plt

In the least perturbed case, where the minimized residual is less than 10^-3, the convergence is plausibly quadratic. At the next level up, the convergence starts similarly but suddenly stagnates for a long time. In the most perturbed case, the quadratic phase is nearly gone and the overall shape looks linear.

Example 4.7.1

g = @(x) [sin(x(1) + x(2)); cos(x(1) - x(2)); exp(x(1) - x(2))];
p = [1; 1];

clf
labels = [];
for R = [1e-3, 1e-2, 1e-1]
    % Define the perturbed function.
    f = @(x) g(x) - g(p) + R * [-1; 1; -1] / sqrt(3)
    x = levenberg(f, [0; 0]);
    r = x(:, end);
    err = abs(x(1, 1:end-1) - r(1));
    normres = norm(f(r));
    semilogy(err), hold on
    labels = [labels; sprintf("R=%.2g", normres)];
end
xlabel("iteration"), ylabel("error")
legend(labels);

Example 4.7.1

g = lambda x: array([sin(x[0] + x[1]), cos(x[0] - x[1]), exp(x[0] - x[1])])
p = array([1, 1])

for R in [1e-3, 1e-2, 1e-1]:
    # Define the perturbed function.
    f = lambda x: g(x) - g(p) + R * array([-1, 1, -1]) / sqrt(3)
    x = FNC.levenberg(f, [0, 0])
    r = x[-1]
    err = [norm(x[j] - r) for j in range(len(x) - 1)]
    normres = norm(f(r))
    semilogy(err, label=f"R={normres:.2g}")
title("Convergence of Gauss–Newton")
xlabel("iteration"), ylabel("error")
legend();

/Users/driscoll/Documents/GitHub/fnc/python/fncbook/fncbook/chapter04.py:167: UserWarning: Iteration did not find a root.
  warnings.warn("Iteration did not find a root.")

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

\mathbf{f}(\mathbf{x}) = \left[\, g(t_i,\mathbf{x})-y_i \, \right]_{\,i=1,\ldots,m}.

(4.7.2)

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

\mathbf{f}{\,}'(\mathbf{x}) = \left[ \frac{\partial}{\partial x_j} g(t_i,\mathbf{x}) \right]_{\,i=1,\ldots,m;\,j=1,\ldots,n.}

(4.7.3)

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.,

g(t,\mathbf{c}) = c_1 g_1(t) + c_2 g_2(t) + \cdots + c_m g_m(t),

(4.7.4)

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

Example 4.7.2 (Nonlinear data fitting)

Inhibited enzyme reactions often follow what are known as Michaelis–Menten kinetics, in which a reaction rate $w$ follows a law of the form

w(s) = \frac{V s}{K_m + s},

(4.7.5)

where $s$ is the concentration of a substrate. The real values $V$ and $K_m$ are parameters that are free to fit to data. For this example, we cook up some artificial data with $V=2$ and $K_m=1/2$ .

Julia

MATLAB

Python

Example 4.7.2

m = 25;
s = range(0.05, 6, length=m)
ŵ = @. 2 * s / (0.5 + s)                      # exactly on the curve
w = @. ŵ + 0.15 * cos(2 * exp(s / 16) * s);     # smooth noise added

scatter(s, w, label="noisy data",
    xlabel="s", ylabel="v", leg=:bottomright)
plot!(s, ŵ, l=:dash, color=:black, label="perfect data")

The idea is to pretend that we know nothing of the origins of this data and use nonlinear least squares to find the parameters in the theoretical model function $v(s)$ . In (4.7.2), the $s$ variable plays the role of $t$ , and $v$ plays the role of $g$ .

function misfit(x)
    V, Km = x   # rename components for clarity
    return @. V * s / (Km + s) - w
end

misfit (generic function with 1 method)

In the Jacobian the derivatives are with respect to the parameters in $\mathbf{x}$ .

function misfitjac(x)
    V, Km = x   # rename components for clarity
    J = zeros(m, 2)
    J[:, 1] = @. s / (Km + s)              # dw/dV
    J[:, 2] = @. -V * s / (Km + s)^2         # dw/d_Km
    return J
end

misfitjac (generic function with 1 method)

x₁ = [1, 0.75]
x = FNC.newtonsys(misfit, misfitjac, x₁)

@show V, Km = x[end]  # final values

(V, Km) = x[end] = [1.968652598378229, 0.4693037307416775]

2-element Vector{Float64}:
 1.968652598378229
 0.4693037307416775

The final values are reasonably close to the values $V=2$ , $K_m=0.5$ that we used to generate the noise-free data. Graphically, the model looks close to the original data.

model(s) = V * s / (Km + s)
plot!(model, 0, 6, label="nonlinear fit")

For this particular model, we also have the option of linearizing the fit process. Rewrite the model as

\frac{1}{w} = \frac{\alpha}{s} + \beta = \alpha \cdot s^{-1} + \beta

for the new fitting parameters $\alpha=K_m/V$ and $\beta=1/V$ . This corresponds to the misfit function whose entries are

f_i([\alpha,\beta]) = \left(\alpha \cdot \frac{1}{s_i} + \beta\right) - \frac{1}{w_i}

(1)

for $i=1,\ldots,m$ . Although this misfit is nonlinear in $s$ and $w$ , it’s linear in the unknown parameters α and β. This lets us pose and solve it as a linear least-squares problem.

A = [s .^ (-1) s .^ 0]
u = 1 ./ w
α, β = A \ u

2-element Vector{Float64}:
 0.12476333709901538
 0.5713959100431231

The two fits are different because they do not optimize the same quantities.

linmodel(x) = 1 / (β + α / x)
plot!(linmodel, 0, 6, label="linearized fit")

The truly nonlinear fit is clearly better in this case. It optimizes a residual for the original measured quantity rather than a transformed one we picked for algorithmic convenience.

Example 4.7.2

m = 25; V = 2; Km = 0.5;
s = linspace(0.05, 6, m)';
w = V * s ./ (Km + s);                      % exactly on the curve
w = w + 0.15 * cos(2 * exp(s / 16) .* s);   % noise added
clf, fplot(@(s) V * s ./ (Km + s), [0, 6], '--')
hold on, scatter(s, w)
xlabel('concentration'), ylabel('reaction rate')    
labels = ["ideal", "noisy data"];    
legend(labels, 'location', 'east');

f47_misfit.m

function [f, J] = misfit(c, s, w)
    V = c(1);   Km = c(2);
    f = V * s ./ (Km + s) - w;
    J(:,1) = s ./ (Km + s);            % d/d(V)
    J(:,2) = -V * s ./ (Km + s).^2;    % d/d(Km)
end

The misfit function above has to know the parameters x that are being optimized as well as the data s and w that remain fixed. We use a closure to pass the data values along.

f = @(x) f47_misfit(x, s, w);

Now we have a function that accepts a single 2-vector input and returns a 25-vector output. We can pass this function to levenberg to find the best-fit parameters.

x1 = [1; 0.75];
x = newtonsys(f, x1);
V = x(1, end),  Km = x(2, end)     % final values
model = @(s) V * s ./ (Km + s);    % best-fit model

The final values are reasonably close to the values $V=2$ , $K_m=0.5$ that we used to generate the noise-free data. Graphically, the model looks close to the original data:

final_misfit_norm = norm(model(s) - w) 
hold on, fplot(model, [0, 6])
title('Michaelis-Menten fitting')    
labels = [labels, "nonlinear fit"];    
legend(labels, 'location', 'east');

For this particular model, we also have the option of linearizing the fit process. Rewrite the model as

\frac{1}{w} = \frac{\alpha}{s} + \beta = \alpha \cdot s^{-1} + \beta

for the new fitting parameters $\alpha=K_m/V$ and $\beta=1/V$ . This corresponds to the misfit function whose entries are

f_i([\alpha,\beta]) = \left(\alpha \cdot \frac{1}{s_i} + \beta\right) - \frac{1}{w_i}

(1)

for $i=1,\ldots,m$ . Although this misfit is nonlinear in $s$ and $w$ , it’s linear in the unknown parameters α and β. This lets us pose and solve it as a linear least-squares problem.

u = 1 ./ w;
A = [s.^(-1), s.^0];  
z = A \ u;
alpha = z(1);  beta = z(2);

The two fits are different because they do not optimize the same quantities.

linmodel = @(s) 1 ./ (beta + alpha ./ s);
final_misfit_linearized = norm(linmodel(s) - w)
fplot(linmodel, [0, 6])
labels = [labels, "linearized fit"];    
legend(labels, 'location', 'east');

The truly nonlinear fit is clearly better in this case. It optimizes a residual for the original measured quantity rather than a transformed one we picked for algorithmic convenience.

Example 4.7.2

m = 25
V, Km = 2, 0.5
s = linspace(0.05, 6, m)
model = lambda x: V * x / (Km + x)
w = model(s) + 0.15 * cos(2 * exp(s / 16) * s)    # noise added

fig, ax = subplots()
ax.scatter(s, w, label="data")
ax.plot(s, model(s), 'k--', label="unperturbed model")
xlabel("s"), ylabel("w")
legend();

def misfit(c):
    V, Km = c  # rename components for clarity
    f = V * s / (Km + s) - w
    return f

In the Jacobian the derivatives are with respect to the parameters in $\mathbf{x}$ .

def misfitjac(x):
    V, Km = x   # rename components for clarity
    J = zeros([m, 2])
    J[:, 0] = s / (Km + s)          # d/d(V)
    J[:, 1] = -V * s / (Km + s)**2  # d/d(Km)
    return J

x1 = [1, 0.75]
x = FNC.newtonsys(misfit, misfitjac, x1)
V, Km = x[-1]  # final values
print(f"estimates are V = {V:.3f}, Km = {Km:.3f}")

estimates are V = 1.969, Km = 0.469

The final values are reasonably close to the values $V=2$ , $K_m=0.5$ that we used to generate the noise-free data. Graphically, the model looks close to the original data.

# since V and Km have been updated, model() is too
ax.plot(s, model(s), label="nonlinear fit")

For this particular model, we also have the option of linearizing the fit process. Rewrite the model as

\frac{1}{w} = \frac{\alpha}{s} + \beta = \alpha \cdot s^{-1} + \beta

for the new fitting parameters $\alpha=K_m/V$ and $\beta=1/V$ . This corresponds to the misfit function whose entries are

f_i([\alpha,\beta]) = \left(\alpha \cdot \frac{1}{s_i} + \beta\right) - \frac{1}{w_i}

(1)

for $i=1,\ldots,m$ . Although this misfit is nonlinear in $s$ and $w$ , it’s linear in the unknown parameters α and β. This lets us pose and solve it as a linear least-squares problem.

from numpy.linalg import lstsq
A = array( [[1 / s[i], 1.0] for i in range(len(s))] )
z = lstsq(A, 1 / w, rcond=None)[0]
alpha, beta = z
print("alpha:", alpha, "beta:", beta)

alpha: 0.12476333709901544 beta: 0.571395910043123

The two fits are different; they do not optimize the same quantities.

linmodel = lambda x: 1 / (beta + alpha / x)
ax.plot(s, linmodel(s), label="linear fit")
ax.legend()
fig

The truly nonlinear fit is clearly better in this case. It optimizes a residual for the original measured quantity rather than a transformed one we picked for algorithmic convenience.

4.7.4Exercises¶

Exercise 4.7.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.

Exercise 4.7.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.

Exercise 4.7.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.

References¶

Kermack, W. O., McKendrick, A. G., & Walker, G. T. (1927). A Contribution to the Mathematical Theory of Epidemics. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 115(772), 700–721. 10.1098/rspa.1927.0118
Wu, Z., Begley, C. G., Situ, P., Simpson, T., & Liu, H. (2014). The Effects of Mild Ocular Surface Stimulation and Concentration on Spontaneous Blink Parameters. Current Eye Research, 39(1), 9–20. 10.3109/02713683.2013.822896
Brosch, J. K., Wu, Z., Begley, C. G., Driscoll, T. A., & Braun, R. J. (2017). Blink Characterization Using Curve Fitting and Clustering Algorithms. Journal for Modeling in Ophthalmology, 1(3), 60–81. 10.35119/maio.v1i3.38
Berke, A., & Mueller, S. (1998). The Kinetics of Lid Motion and Its Effects on the Tear Film. In D. A. Sullivan, D. A. Dartt, & M. A. Meneray (Eds.), Lacrimal Gland, Tear Film, and Dry Eye Syndromes 2 (Vol. 438, pp. 417–424). Springer US. 10.1007/978-1-4615-5359-5_58