Python scipy optimize minimize and curve fitting

When you need to minimize a scalar function or fit a model to data, Python's scipy.optimize module provides two primary entry points: minimize for general optimization and curve_fit for nonlinear least squares. This article covers python scipy optimize minimize and curve fitting, explaining how to use both, when each is appropriate, and how to handle common issues.

Using minimize for General Optimization

scipy.optimize.minimize solves problems of the form:

[ \min_x f(x) ]

where (f) is a scalar-valued function of one or more variables. The function takes an initial guess x0 and returns an OptimizeResult containing the optimal values, the objective value, convergence status, and other diagnostics.

Here is a minimal example that minimizes the Rosenbrock function, a common test problem:

import numpy as np
from scipy.optimize import minimize

def rosen(x):
    return sum(100.0 * (x[1:] - x[:-1]**2)**2 + (1 - x[:-1])**2)

x0 = np.array([-1.2, 1.0])
result = minimize(rosen, x0, method='BFGS')
print(result.x)
print(result.fun)

The method parameter selects the optimization algorithm. Common choices include 'BFGS' for unconstrained problems, 'L-BFGS-B' for bounded problems, 'SLSQP' for constrained problems, and 'trust-constr' for problems with both bounds and constraints. The default method depends on whether bounds or constraints are present.

When the objective function is expensive to evaluate, you can supply a gradient via the jac parameter. For example:

def rosen_grad(x):
    grad = np.zeros_like(x)
    grad[0] = -400 * x[0] * (x[1] - x[0]**2) - 2 * (1 - x[0])
    grad[1] = 200 * (x[1] - x[0]**2)
    return grad

result = minimize(rosen, x0, jac=rosen_grad, method='BFGS')

Providing an analytic gradient often improves convergence speed and reliability.

Fitting a Model to Data with curve_fit

scipy.optimize.curve_fit is a specialized tool for nonlinear least squares fitting. Given a model function f(x, *params) and observed data (xdata, ydata), it finds the parameter values that minimize the sum of squared residuals:

[ \min_{\theta} \sum_i (y_i - f(x_i, \theta))^2 ]

The basic usage is straightforward:

from scipy.optimize import curve_fit

def model(x, a, b, c):
    return a * np.exp(-b * x) + c

xdata = np.linspace(0, 4, 50)
ydata = model(xdata, 2.5, 1.3, 0.5) + 0.2 * np.random.normal(size=50)

popt, pcov = curve_fit(model, xdata, ydata, p0=[1, 1, 1])
print(popt)

The function returns two values: popt, the optimal parameter values, and pcov, the estimated covariance of the parameters. The square root of the diagonal of pcov gives the standard errors of the fitted parameters.

You can constrain parameters using the bounds argument. It accepts a tuple of two lists: the lower and upper bounds for each parameter. Use -np.inf and np.inf for unbounded parameters:

popt, pcov = curve_fit(model, xdata, ydata, p0=[1, 1, 1],
                       bounds=([0, 0, -np.inf], [np.inf, np.inf, np.inf]))

Key Differences Between minimize and curve_fit

curve_fit is essentially a wrapper around scipy.optimize.least_squares that handles the residual computation and parameter unpacking. If your problem is purely about fitting a model to data, curve_fit is usually more convenient. If you need to optimize an arbitrary objective function, minimize is the appropriate tool.

Handling Bounds and Constraints

For minimize, bounds are passed as a sequence of (min, max) pairs for each variable. For example:

result = minimize(rosen, x0, method='L-BFGS-B',
                  bounds=[(-2, 2), (-2, 2)])

General constraints are defined using LinearConstraint or NonlinearConstraint objects and passed to the constraints parameter. Here is an example with a nonlinear inequality constraint:

from scipy.optimize import NonlinearConstraint

def constraint(x):
    return x[0]**2 + x[1]**2

nlc = NonlinearConstraint(constraint, 0, 1)
result = minimize(rosen, x0, method='trust-constr', constraints=[nlc])

curve_fit only supports box bounds, not general constraints. If you need to enforce a relationship between parameters, you must reparameterize the model or use minimize with a custom residual function.

Choosing the Right Method

The performance and reliability of minimize depend heavily on the chosen method. For unconstrained smooth problems, BFGS is a good default. If the problem has many variables and you need to conserve memory, L-BFGS-B is often faster. For constrained problems, SLSQP is a solid choice, but trust-constr is more robust for nonlinear constraints.

curve_fit uses the Levenberg-Marquardt algorithm by default, which is efficient for small and medium-sized problems. You can change the method via the method parameter, but the default is usually adequate. If your model is highly nonlinear or the initial guess is poor, consider using method='trf' (Trust Region Reflective) which is more robust to outliers.

Numerical Considerations and Performance

Scaling matters. If parameters differ by orders of magnitude, the optimizer may struggle. Consider rescaling variables or providing a custom gradient. For curve_fit, the sigma parameter can be used to weight residuals, which is important when data points have different uncertainties.

Tolerances control when the optimizer stops. minimize accepts tol for the gradient norm, and curve_fit accepts ftol, xtol, and gtol via the maxfev and epsfcn parameters. Tightening tolerances too much can lead to excessive iterations without meaningful improvement.

When the objective function is expensive, consider supplying a Jacobian (for minimize) or a residual Jacobian (for curve_fit via the jac parameter). This can reduce the number of function evaluations significantly.

Common Pitfalls and How to Avoid Them

A poor initial guess is the most common cause of failure. For curve_fit, start with parameter values that are plausible given the data. Plotting the model with the initial guess against the data can quickly reveal bad starting points.

Local minima are a risk in non-convex problems. Run the optimization from multiple starting points and compare the resulting objective values. For curve_fit, you can use p0 as a list of multiple guesses and loop over them.

Singular covariance matrices indicate that the model is overparameterized or that parameters are not identifiable from the data. Reduce the number of parameters or add regularization.

Finally, remember that curve_fit assumes the residuals are independent and normally distributed. If that assumption is violated, the covariance estimates may be misleading.