least_squares

ndoptimize.least_squares(model, x, data, variance, guess, lower=-inf, upper=inf, prior_mean=0, prior_sigma=inf, correlation=None, max_iterations=1000, tolerance=1e-08)

Fit a model to many independent sets of data at once, using the Levenberg-Marquardt method.

Every fit has the same model, but its own independent variable, data, noise, initial guess, bounds and priors. The fits run in parallel, each with its own damping and convergence test.

Each fit minimizes

\[\chi^2 = (y - f)^T C^{-1} (y - f) + \sum_k \left( \frac{p_k - \bar{p}_k}{\sigma_k} \right)^2,\]

where \(y\) is the data, \(f\) is the model, \(p_k\) are the parameters, \(\bar{p}_k\) and \(\sigma_k\) are the means and standard deviations of their Gaussian priors, and \(C\) is the covariance of the noise,

\[C_{ij} = \sigma_i \sigma_j \rho_{|i - j|},\]

with \(\sigma_i^2\) the variance of data point \(i\) and \(\rho_k\) the correlation of the noise of two points \(k\) samples apart. Uncorrelated noise, \(\rho = (1,)\), is ordinary weighted least squares. Correlated noise, from resampling for instance, makes this generalized least squares, which is computed by whitening the residual and the Jacobian with the banded Cholesky factor of \(C\).

Parameters:

• model (Model) – The model to fit, made with ndoptimize.model().

• x (ArrayLike) – The independent variable, an array of shape (..., num_points, num_inputs), where num_inputs is the number of names in ndoptimize.Model.inputs.

• data (ArrayLike) – The data to fit, an array of shape (..., num_points). Points which are NaN are left out of their fit.

• variance (ArrayLike) – The variance of the noise of each data point, an array of shape (..., num_points). Points with a variance which is NaN or not positive are left out of their fit.

• guess (ArrayLike) – The initial guess for the parameters, an array of shape (..., num_parameters), where num_parameters is the number of names in ndoptimize.Model.parameters. It is clipped to the bounds.

• lower (ArrayLike) – The lower bound of each parameter, an array of shape (..., num_parameters).

• upper (ArrayLike) – The upper bound of each parameter, an array of shape (..., num_parameters), no less than lower.

• prior_mean (ArrayLike) – The mean of the Gaussian prior of each parameter, an array of shape (..., num_parameters). It is ignored where prior_sigma is infinite.

• prior_sigma (ArrayLike) – The standard deviation of the Gaussian prior of each parameter, a positive array of shape (..., num_parameters). An infinite standard deviation, the default, means no prior.

• correlation (None | ArrayLike) – The correlation of the noise of two data points by their separation along the last axis of data, in samples, starting with 1 for no separation and zero beyond its length. The separation counts the points left out of a fit. It must describe a positive-definite noise covariance. If None (the default), the noise is uncorrelated.

• max_iterations (int) – The maximum number of accepted steps for each fit. A step which does not decrease \(\chi^2\) is retried with more damping, which evaluates the model again without counting as a step.

• tolerance (float) – The convergence tolerance on \(g^T A^{-1} g\), the decrease in \(\chi^2\) that a Gauss-Newton step predicts, where \(g = J^T C^{-1} r\) is minus half the gradient of \(\chi^2\) and \(A = J^T C^{-1} J\) is the normal matrix, both including the priors. This is half the square of the Newton decrement of \(\chi^2\). The default, \(10^{-8}\), puts each parameter within about \(10^{-4}\) of its uncertainty of the minimum.

Returns:

• The best-fit parameters, their covariance, chi-squared, and the

• convergence state of each fit.

• Fits with no usable data points are NaN.

Return type:

LeastSquaresResult

Notes

The batch axes of the arguments, all but their last one (two for x), are broadcast against each other.

Each step solves \((A + \lambda \, \text{diag} A) \delta = g\), Marquardt’s damping, and the damping parameter \(\lambda\) is updated by the gain-ratio rule of Madsen, Nielsen, and Tingleff (2004), “Methods for Non-Linear Least Squares Problems”, which converges in curved valleys of \(\chi^2\) (where the amplitude and width of a weak line trade off) in about a third of the steps of the classic fixed factors.

Steps are clipped to the bounds, and a parameter at a bound which the gradient pushes outward is held fixed for the step and for the convergence test, as is a parameter which has no effect on the model.

A fit whose \(\chi^2\) is NaN at the initial guess, where the model is NaN for instance, takes no steps and does not converge.

Examples

Fit a Gaussian plus a straight line to a thousand noisy spectral lines.

import numpy as np
import matplotlib.pyplot as plt
import ndoptimize

rng = np.random.default_rng(0)

# the true parameters: amplitude, center, width, slope and intercept
num = 1000
truth = np.stack([
    rng.uniform(2, 10, num),
    rng.normal(0, 5, num),
    rng.uniform(10, 20, num),
    np.zeros(num),
    rng.uniform(1, 2, num),
], axis=-1)

# the noisy spectra
velocity = np.linspace(-100, 100, 41)
x = velocity[:, np.newaxis]
f = np.empty((num, velocity.size))
jacobian = np.empty((velocity.size, 5))
for i in range(num):
    ndoptimize.models.gaussian_linear(x, truth[i], f[i], jacobian)
variance = 0.5
data = f + rng.normal(0, np.sqrt(variance), f.shape)

# fit every spectrum
result = ndoptimize.least_squares(
    model=ndoptimize.models.gaussian_linear,
    x=x,
    data=data,
    variance=variance,
    guess=[5, 0, 15, 0, 1],
)

# compare the fitted centers to the truth
sigma = np.sqrt(result.covariance[..., 1, 1])
error = (result.parameters[..., 1] - truth[..., 1]) / sigma

fig, ax = plt.subplots(constrained_layout=True)
ax.hist(error, bins=40, density=True, label="fits");
z = np.linspace(-4, 4, 101)
ax.plot(z, np.exp(-z**2 / 2) / np.sqrt(2 * np.pi), label="unit normal");
ax.set_xlabel("(fitted center - true center) / uncertainty");
ax.legend();