Nonlinear Curve Fitting with scipy.optimize.curve fit

Nonlinear curve fitting is one of the most practically important tools in scientific computing and data analysis. While linear regression and polynomial fitting (via numpy.polyfit) handle a broad class of problems, many real-world relationships between variables are fundamentally nonlinear — exponential decay, logistic growth, sinusoidal oscillations, Gaussian peaks, power laws, and countless others cannot be adequately described by a polynomial or a linear combination of predictors. The function scipy.optimize.curve_fit is the standard Python tool for fitting arbitrary nonlinear model functions to observed data, and understanding it deeply means understanding not just the API, but the optimization strategy it uses, how to set it up properly, and how to diagnose and recover from failure.

When to Use Nonlinear Fitting

The essential criterion for reaching for nonlinear fitting is that your model contains parameters that cannot be extracted by linear algebra alone. A linear model, in the statistical sense, is one where the predicted value is a linear combination of known basis functions — for example, y = a + b*x + c*x**2 is linear in the parameters a, b, and c even though it curves as a parabola. You can write such a system as a matrix equation and solve it directly. By contrast, a model like y = A * exp(-k * x) is nonlinear in the parameter k — you cannot isolate k by simple matrix inversion. The same applies to y = A * sin(omega * x + phi), where omega and phi appear inside the transcendental function.

Internally, scipy.optimize.curve_fit solves the nonlinear least-squares problem by an iterative numerical algorithm — by default, the Levenberg-Marquardt method (with fallback to Trust Region Reflective when bounds are specified). The optimizer repeatedly evaluates your model function with trial parameter values, computes the residuals (differences between model predictions and observed data), and adjusts the parameters in a direction that reduces the sum of squared residuals. Because the problem is not convex in general, the algorithm must start somewhere in parameter space, and that starting point is supplied by the user through the p0 argument. This is a crucial difference from linear fitting: the result you get depends on where you start, and poor starting points can lead to a locally optimal but globally wrong solution, or outright failure to converge.

Nonlinear fitting is more flexible than polynomial fitting but demands more from the user. You must encode your physical or mathematical intuition about the relationship into a specific functional form, supply reasonable initial guesses, optionally provide parameter bounds, and interpret the uncertainty estimates produced by the fit. Done carefully, it produces highly interpretable results with direct physical meaning attached to each parameter.

Defining a Model Function for curve_fit

The model function is the heart of the fitting procedure. Its signature must follow a strict convention: the first argument is the independent variable (often written x, but it can be named anything and can even be a 2D array for fitting functions of multiple variables), followed by the parameters one by one as separate positional arguments. The function must return a NumPy array of predicted values with the same shape as the input x.

import numpy as np
from scipy.optimize import curve_fit

# Exponential decay model: y = A * exp(-k * x) + C
def exp_decay(x, A, k, C):
    return A * np.exp(-k * x) + C

Notice the use of np.exp rather than Python's math.exp. NumPy universal functions (ufuncs) operate element-wise on entire arrays, so when x is a NumPy array of 100 data points, np.exp(-k * x) returns an array of 100 values in a single vectorized call. If you used math.exp, it would raise a TypeError when handed an array. Any NumPy ufunc — np.sin, np.cos, np.log, np.sqrt, and so on — is appropriate here.

The parameter names A, k, and C are purely for human readability and code clarity. curve_fit does not parse argument names; it identifies parameters entirely by their positional order. The first extra argument after x is parameter 0 in popt, the second is parameter 1, and so forth. This means the order you write the parameters in the function definition determines the order they appear in the output and in your initial guess list.

A slightly more complex example illustrates a Gaussian peak, common in spectroscopy and signal analysis:

# Gaussian model: y = A * exp(-(x - mu)^2 / (2 * sigma^2))
def gaussian(x, A, mu, sigma):
    return A * np.exp(-(x - mu)**2 / (2 * sigma**2))

Here A is the peak amplitude, mu is the center position, and sigma is the standard deviation (width). All three are meaningful physical quantities that the fit will estimate from data. Writing such model functions is an exercise in translating equations directly into code, ensuring every operation is array-compatible.

Calling curve_fit and Interpreting Output

curve_fit returns a tuple of exactly two objects: popt and pcov.

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# Simulate noisy exponential decay data
np.random.seed(42)
x_data = np.linspace(0, 4, 50)
y_true = 3.0 * np.exp(-1.5 * x_data) + 0.5
y_data = y_true + np.random.normal(0, 0.2, size=x_data.size)

def exp_decay(x, A, k, C):
    return A * np.exp(-k * x) + C

p0 = [2.5, 1.0, 0.3]  # initial guesses for A, k, C
popt, pcov = curve_fit(exp_decay, x_data, y_data, p0=p0)

print("Optimal parameters:", popt)
print("Covariance matrix:\n", pcov)

popt is a 1D NumPy array containing the parameter values at the least-squares optimum — one value per parameter in the order they appear in your model function. In this example, popt[0] is the fitted A, popt[1] is the fitted k, and popt[2] is the fitted C. To generate the fitted curve and overlay it on the data, simply call your model function with popt unpacked:

x_fit = np.linspace(0, 4, 300)
y_fit = exp_decay(x_fit, *popt)  # unpacks popt as A, k, C

plt.scatter(x_data, y_data, label='Data', color='steelblue', s=20)
plt.plot(x_fit, y_fit, label='Fitted curve', color='tomato')
plt.legend()
plt.show()

pcov is a 2D NumPy array — the estimated covariance matrix of the parameters. Its shape is (n_params, n_params). The diagonal entries pcov[i, i] give the estimated variance of parameter i. Taking the square root of the diagonal gives the one-sigma (one standard deviation) uncertainty on each parameter:

perr = np.sqrt(np.diag(pcov))
print("Parameter uncertainties (1-sigma):")
for name, val, err in zip(['A', 'k', 'C'], popt, perr):
    print(f"  {name} = {val:.4f} ± {err:.4f}")

These uncertainties reflect how precisely the data constrain each parameter under the assumption that the model is correct and that the residuals are approximately normally distributed. The off-diagonal elements of pcov encode correlations between parameters — a large off-diagonal value means that the data cannot independently pin down both parameters; if one shifts, the other shifts to compensate. High parameter correlation is often a sign that the model is overparameterized or that the data do not span a sufficient range to disentangle the effects of different parameters.

It is worth emphasizing that if curve_fit cannot estimate the covariance (for example, because the Jacobian is singular near the solution), entries of pcov will be inf. This is a warning that the fit result may be unreliable even if popt looks plausible.

Setting Initial Guesses and Bounds

The p0 argument is a list or array of initial guesses for each parameter, in the same order as the function signature. If p0 is not provided, curve_fit defaults to all ones, which is often a poor choice. Because nonlinear optimization is sensitive to starting conditions, choosing good initial guesses is genuinely important — it is not a technicality to be glossed over.

The recommended approach is to use domain knowledge and visual inspection together. Plot the raw data first and ask: what is the rough amplitude of the signal? Where does it peak or level off? What timescale or frequency seems to govern the behavior? Translate those visual impressions directly into numerical guesses.

# Suppose the data looks like it starts near 3, decays fast, and flattens near 0.5
# A rough guess: A ~ 3, k ~ 1.5, C ~ 0.5
p0 = [3.0, 1.5, 0.5]
popt, pcov = curve_fit(exp_decay, x_data, y_data, p0=p0)

A powerful diagnostic before calling curve_fit is to overlay your model evaluated at your guessed parameters on the raw data. If the guessed curve is qualitatively in the right neighborhood — roughly the right shape, magnitude, and location — the optimizer has a good chance of converging to the correct minimum. If the guessed curve looks completely wrong, refine your guesses until it looks reasonable.

# Check initial guess visually before fitting
y_guess = exp_decay(x_data, *p0)
plt.scatter(x_data, y_data, label='Data')
plt.plot(x_data, y_guess, '--', label='Initial guess')
plt.legend()
plt.show()

The bounds argument allows you to constrain parameters to physically meaningful ranges. It is specified as a tuple of two sequences: (lower_bounds, upper_bounds), each of length equal to the number of parameters. Use -np.inf and np.inf to leave a parameter unconstrained on one side.

# Constrain A > 0, k > 0, C can be anything
lower = [0, 0, -np.inf]
upper = [np.inf, np.inf, np.inf]
popt, pcov = curve_fit(exp_decay, x_data, y_data, p0=p0, bounds=(lower, upper))

Bounds serve two purposes: they prevent the optimizer from wandering into regions where the model is undefined or physically nonsensical (e.g., a negative decay rate or a negative amplitude when the data are always positive), and they can help guide convergence when the loss landscape has multiple minima. When bounds are active, curve_fit automatically switches from the Levenberg-Marquardt algorithm to the Trust Region Reflective algorithm, which natively handles box constraints. This is done transparently — you do not need to specify the method yourself unless you want to.

Handling Convergence Failures

The most common failure mode of curve_fit is a RuntimeError with a message such as "Optimal parameters not found: Number of calls to function has reached maxfev". This means the iterative optimizer exhausted its budget of function evaluations without satisfying the convergence criterion. It does not necessarily mean a solution does not exist; it usually means the algorithm got lost on its way to the solution.

The first thing to do when this happens is not to immediately increase maxfev. That is a last resort. The first response should be to diagnose the starting point:

# Step 1: plot your data
plt.scatter(x_data, y_data)
plt.show()

# Step 2: manually evaluate the model at your initial guess and plot it
y_guess = exp_decay(x_data, *p0)
plt.scatter(x_data, y_data, label='Data')
plt.plot(x_data, y_guess, '--', label='Guess')
plt.legend()
plt.show()

If the guessed curve looks wildly different from the data, refine your p0 until it qualitatively matches. Even rough agreement — within an order of magnitude on each parameter — is usually enough for the algorithm to find its way. Also consider whether the scale of parameters is very different from one to another (e.g., one parameter is around 1000 and another is around 0.001), because large scale disparities can make the optimizer's gradient estimates unreliable. In such cases, rescaling your data or your model parameters can help.

If you have genuinely tried good initial guesses and the optimizer still struggles, you can increase the maximum number of allowed function evaluations:

popt, pcov = curve_fit(exp_decay, x_data, y_data, p0=p0, maxfev=10000)

The default value of maxfev is 200 * (number_of_parameters + 1), which is sometimes insufficient for poorly conditioned problems or models with many parameters. Raising it gives the optimizer more attempts to converge, but if it still fails after a very large maxfev, the problem almost certainly lies elsewhere.

A repeated or persistent convergence failure should prompt a deeper reconsideration of the model itself. Ask whether the chosen functional form is truly appropriate for the data. A model that is structurally wrong will never fit, no matter how sophisticated the optimizer or how good the initial guess. It is also worth checking whether the data contain outliers that are pulling the fit toward implausible parameter values — robust fitting methods or outlier removal may be necessary in such cases. Finally, parameter identifiability should be considered: some model forms have redundant parameters that cannot be uniquely determined from the data (for example, if two parameters always appear together as a product, only their product is identifiable, not the individual values). In such cases, the model should be reparameterized to remove the redundancy before fitting.

In summary, effective use of scipy.optimize.curve_fit is a workflow that combines mathematical modeling (writing an appropriate model function), physical intuition (choosing good initial guesses and meaningful bounds), numerical execution (calling the function and unpacking popt and pcov), and statistical interpretation (reading off parameter uncertainties and assessing the quality of the fit through residual analysis). Each step informs the others, and skill in nonlinear curve fitting grows with practice across a variety of real data sets and model types.