Nonlinear least-squares regression - Executable Engineering: Interactive Numerical Methods in Python

There is one more important unconstrained minimization that we should discuss: non-linear least-squares, which is the extension of the linear least-square regression we did for curve fitting.

Now, we have several measurements in pairs $(x_i, y_i)$ and we wish to fit a function of $x$ which depends in some complex, nonlinear fashion on a set of parameters $\beta$ . E.g.: $f(x,\beta) = \sin(\beta_1 x) e^{\beta_2 x^3}$

For each data point we have, we can define a residual, $r_i = f(x_i,\beta) - y_i$ . Summing the square of these residuals, we get $F = \mathbf{r} \cdot \mathbf{r}$ . Our goal is to choose the set of parameters $\beta$ that minimizes $F$ :

\min_\beta F = \min_\beta \sum_i r_i^2 = \min_\beta \sum_i [f(x_i,\beta) - y_i]^2

(1)

Note: $F$ is necessarily symmetric and positive semidefinite, therefore any extremum is a minimum.

import numpy as np
from scipy.optimize import curve_fit
import plotly.graph_objects as go
from IPython.display import Math, display

try:
    import executable_engineering as exe
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

Gauss-Newton regression¶

Let’s consider Newton’s optimization:

\frac{\partial^2 F}{\partial \beta^2} \Delta \beta = - \frac{\partial F}{\partial \beta}

(2)

Since $F$ is the sum of squared $r_i$ we can expand:

\frac{\partial F}{\partial \beta} = 2 r \frac{\partial r}{\partial \beta} = 2 r J

(3)

where $J = \frac{\partial r}{\partial \beta}$ is the Jacobian of the residuals with respect to the parameters. Since we typically have more data points (residuals) than parameters, this is rectangular.

The Hessian is expanded:

\begin{aligned} \frac{\partial^2 F}{\partial \beta^2} &= 2 \sum_i \frac{\partial r_i}{\partial \beta_j} \frac{\partial r_i}{\partial \beta_k} + r \frac{\partial^2 r_i}{\partial \beta_j \partial \beta_k} \\ &= 2\left[ J^T J + r \frac{\partial^2 r_i}{\partial \beta_j \partial \beta_k} \right] \end{aligned}

(4)

The second derivative in the Hessian is troublesome for computation. Luckily, near the root $r\rightarrow 0$ and we can disregard the term entirely:

\frac{\partial^2 F}{\partial \beta^2} \approx 2 J^T J

(5)

Our minimization now becomes:

\begin{aligned} J^T J \Delta \beta &= -J^T r \\ \Delta \beta &= -[J^T J]^{-1} J^T r \\ &= -J^\dagger r \end{aligned}

(6)

where $J^\dagger$ is the pseudoinverse of the Jacobian!

This is an interesting result since it suggests we could write the original problem as solving $J \Delta \beta = -r$ which is the Newton-Raphson method for root finding, not minimization... Remember however, that this is the pseudoinverse, not the true inverse. Indeed, if we had exactly as many data points as parameters, $J$ would be square, and we would be able to find the exact solution which is the root! The fact that this is analogously true for least squares regression comes from 1) approximating the Hessian as the product of Jacobians, and 2) properties of the pseudo-inverse in performing linear least-squares.

The Levenberg-Marquardt algorithm¶

The problem with Gauss-Newton is that the Hessian approximation is only good near the root, and the method fails unless the initial guess is good.

Levenberg-Marquardt blends Gauss-Newton with Gradient descent:

[J^T J + \lambda I ] \Delta \beta = -J^T r

(7)

When $\lambda$ is small, the method reduces to Gauss-Newton. When it is large, the $J^T J$ term can be disregarded and the method is simply gradient descent. Determination of $\lambda$ is a matter of heuristics and implementation.

The Levenberg-Marquardt algorithm is the go-to for common least-squares solvers including SciPy, NumPy, and tools like Excel.

Example: Nonlinear curve fit¶

We fit a damped sinusoidal model to noisy observations.

The original (data-generating) function is:

f_{\text{orig}}(x) = 2.5 \sin(1.3 x) e^{-0.5 x}

(8)

We fit the parameterized model function:

f(x; a, b, c) = a \sin(b x) e^{-c x}

(9)

Nonlinear least-squares determines the optimal parameters $\hat{a}, \hat{b}, \hat{c}$ for the fit function:

f_{\text{fit}}(x) = \hat{a} \sin(\hat{b} x) e^{-\hat{c} x}

(10)

# Define the function to fit
def func(x, a, b, c):
    return a * np.sin(b * x) * np.exp(-c * x)

# Generate sample data with reproducible noise
np.random.seed(42)
x_data = np.linspace(0, 5, 50)
y_data = func(x_data, 2.5, 1.3, 0.5) + np.random.normal(0, 0.2, 50)

# Perform the curve fit
popt, _ = curve_fit(func, x_data, y_data)

# Print and display the fitted parameters and function
print("Fitted parameters [a, b, c]:", np.round(popt, 4))
display(Math(rf"f_{{\text{{fit}}}}(x) = {popt[0]:.3f} \sin({popt[1]:.3f} x) e^{{-{popt[2]:.3f} x}}"))

# Generate the fitted curve and true curve
y_fit = func(x_data, *popt)
y_true = func(x_data, 2.5, 1.3, 0.5)

# Plot the data, original function, and fitted curve
fig = go.Figure()
fig.add_trace(go.Scatter(x=x_data, y=y_true, mode='lines', name='Original: 2.5 sin(1.3x) e^{-0.5x}', line=dict(dash='dash', color='green')))
fig.add_trace(go.Scatter(x=x_data, y=y_data, mode='markers', name='Noisy data', marker=dict(color='blue')))
fig.add_trace(go.Scatter(x=x_data, y=y_fit, mode='lines', name='Fit', line=dict(color='red')))
fig.update_layout(
    xaxis_title='x',
    yaxis_title='y',
    template='plotly_white',
    width=700,
    height=450
)
fig.show()

Example: Nonlinear 2D curve fit¶

We extend nonlinear regression to higher dimensions by fitting a 2D surface.

The original (data-generating) function is a standard 2D Gaussian centered at the origin:

f_{\text{orig}}(x, y) = \exp\left( -\frac{x^2 + y^2}{2} \right)

(11)

We fit an elliptical 2D Gaussian model with a vertical baseline offset:

f(x, y; a, b, c, d, e, f) = a \exp\left( -\left[ \frac{(x - b)^2}{2 c^2} + \frac{(y - d)^2}{2 e^2} \right] \right) + f

(12)

Nonlinear least-squares estimates $\hat{\boldsymbol{\beta}} = [\hat{a}, \hat{b}, \hat{c}, \hat{d}, \hat{e}, \hat{f}]^T$ to obtain the fit function:

f_{\text{fit}}(x, y) = \hat{a} \exp\left( -\left[ \frac{(x - \hat{b})^2}{2 \hat{c}^2} + \frac{(y - \hat{d})^2}{2 \hat{e}^2} \right] \right) + \hat{f}

(13)

# Define the 2D Gaussian model function with offset
def func(X, a, b, c, d, e, f):
    x, y = X
    return a * np.exp(-((x - b)**2 / (2 * c**2) + (y - d)**2 / (2 * e**2))) + f

# Generate sample data with reproducible noise
np.random.seed(42)
x_data = np.linspace(-5, 5, 50)
y_data = np.linspace(-5, 5, 50)
X, Y = np.meshgrid(x_data, y_data)
Z_data = func((X, Y), 1, 0, 1, 0, 1, 0) + np.random.normal(0, 0.1, (50, 50))

# Flatten coordinates and grid data for curve_fit
x_data_flat = X.flatten()
y_data_flat = Y.flatten()
z_data_flat = Z_data.flatten()

# Perform curve fit with initial guess
initial_guess = [1, 0, 1, 0, 1, 0]
popt, _ = curve_fit(func, (x_data_flat, y_data_flat), z_data_flat, p0=initial_guess)

print("Fitted parameters [a, b, c, d, e, f]:", np.round(popt, 4))
x_shift = f"+ {-popt[1]:.3f}" if popt[1] < 0 else f"- {popt[1]:.3f}"
y_shift = f"+ {-popt[3]:.3f}" if popt[3] < 0 else f"- {popt[3]:.3f}"
offset_str = f"- {-popt[5]:.3f}" if popt[5] < 0 else f"+ {popt[5]:.3f}"
display(Math(rf"f_{{\text{{fit}}}}(x, y) = {popt[0]:.3f} \exp\left( -\left[ \frac{{(x {x_shift})^2}}{{2({popt[2]:.3f})^2}} + \frac{{(y {y_shift})^2}}{{2({popt[4]:.3f})^2}} \right] \right) {offset_str}"))

# Generate the fitted surface
Z_fit = func((X, Y), *popt)

# Plot the sample data and fitted surface
fig = go.Figure()
fig.add_trace(go.Scatter3d(
    x=x_data_flat[::4], y=y_data_flat[::4], z=z_data_flat[::4],
    mode='markers', marker=dict(size=2, color='blue', opacity=0.5),
    name='Noisy data'
))
fig.add_trace(go.Surface(
    x=x_data, y=y_data, z=Z_fit,
    colorscale='viridis', opacity=0.7,
    name='Fit surface'
))
fig.update_layout(
    scene=dict(xaxis_title='x', yaxis_title='y', zaxis_title='z'),
    width=750, height=550, margin=dict(l=20, r=20, t=30, b=20)
)
fig.show()