SciPy least_squares: Methods, Loss Functions and Results

scipy.optimize.least_squares solves nonlinear least-squares problems: you give it a function that returns a vector of residuals and a starting guess, and it adjusts the parameters to minimise the sum of squared residuals. Unlike curve_fit, it lets you choose a robust loss that down-weights outliers, pick the solver (trf, dogbox or lm), and read back why it stopped. This guide covers “scipy least squares” usage with a runnable outlier-robust fit, a method decision table, and a decoder for the result object. It follows the SciPy v1.18.0 reference and tutorial.

What scipy least_squares is used for in research

The SciPy reference defines the problem as minimising F(x) = 0.5 * sum(ρ(f_i(x)²)) subject to lower and upper bounds on x, where f is the vector of residuals and ρ is a loss function. With the default loss='linear', ρ(z) = z, so this is ordinary nonlinear least squares. Typical research uses are fitting a decay, growth or dose-response model to measurements, calibrating a simulation against data, and any problem where you can write “model minus observation” as a function of the parameters.

The SciPy tutorial writes the data-fitting case as residuals f_i(x) = w_i(φ(t_i; x) – y_i), where φ is the model, y_i the observation and w_i an optional observation weight. You return that vector from your function; least_squares never sees the data except through it.

If you only need parameters and their uncertainties for a single curve, start with the CASRAI guide to SciPy curve_fit. Use least_squares when you need robust loss, a choice of method, sparse Jacobians, or the full diagnostics. For the wider module, see scipy.optimize.

Install or access

The official SciPy install page (scipy.org/install, checked 2026-10-29) gives python -m pip install scipy and says to install Python first and work inside a virtual environment. Then import the function:

import numpy as np
from scipy.optimize import least_squares

The signature in SciPy v1.18.0 is least_squares(fun, x0, jac='2-point', bounds=(-inf, inf), method='trf', ftol=1e-08, xtol=1e-08, gtol=1e-08, x_scale=None, loss='linear', f_scale=1.0, diff_step=None, tr_solver=None, tr_options=None, jac_sparsity=None, max_nfev=None, verbose=0, args=(), kwargs=None, callback=None, workers=None). The callback and workers arguments were added in 1.16.0, and the default of x_scale changed from 1 to None in 1.16.0, so older installs behave differently.

Step-by-step: a robust fit with outliers

The task: fit y = a + b·exp(c·t) to 30 points that include three outliers, and compare losses. The reference documents this same model and the f_scale=0.1 setting; the data below are deterministic (no random numbers), so you can reproduce every printed value.

1. Write the residual function. fun(x, *args, **kwargs) must return the residual vector, with parameters first. Here: def residuals(x, t, y): return (x[0] + x[1] * np.exp(x[2] * t)) - y.
2. Choose a starting guess x0 (an array of shape (n,) or a float).
3. Pass your data through args=(t, y) and pick a loss and f_scale.
4. Call least_squares and read result.x, result.cost and result.status.

import numpy as np
from scipy.optimize import least_squares

def model(x, t):
    return x[0] + x[1] * np.exp(x[2] * t)

def residuals(x, t, y):
    return model(x, t) - y

t = np.linspace(0, 10, 30)
true = np.array([0.5, 2.0, -0.4])
y = model(true, t) + np.array([0.02, -0.03, 0.01, 0.0, -0.02, 0.03] * 5)
y[[4, 15, 22]] += [3.0, -2.5, 4.0]      # three outliers

x0 = np.array([1.0, 1.0, 0.0])
for loss in ["linear", "soft_l1", "huber", "cauchy"]:
    r = least_squares(residuals, x0, loss=loss, f_scale=0.1, args=(t, y))
    print(loss, np.round(r.x, 3), "cost=%.4f" % r.cost, "status", r.status, "nfev", r.nfev)

Output from this script (computed for this guide with SciPy 1.18.0 and NumPy 2.5.1; the true parameters are 0.5, 2.0, -0.4):

loss	fitted a, b, c	cost	status	nfev
linear	0.645, 2.227, -0.446	15.0831	2	21
soft_l1	0.507, 2.012, -0.402	0.9219	2	12
huber	0.506, 2.011, -0.402	0.9365	2	12
cauchy	0.504, 2.000, -0.401	0.1089	2	12

The linear fit is pulled toward the three outliers; every robust loss recovers parameters close to the true values. The costs are not comparable across losses, because each loss defines cost differently (see the next section). The reference suggests trying soft_l1 or huber first and notes that cauchy may cause optimisation difficulties; it behaved here, but treat that as one dataset.

Bounds

Pass bounds=(lower, upper) as two arrays or a Bounds instance; use np.inf for unbounded sides. In the same example, bounds=([0, 0, -1], [1, 1.5, 0]) with soft_l1 stops with b at its upper limit of 1.5 (fitted a, b, c = 0.547, 1.5, -0.343, cost 1.0696 with the data above), and result.active_mask reports it as [0 1 0].

Supplying the Jacobian

The default jac='2-point' estimates the Jacobian by finite differences. '3-point' is more accurate at about twice the cost, and 'cs' needs a function that accepts complex input. The SciPy tutorial states it is “highly recommended to compute this matrix analytically and pass it to least_squares” because finite differences can be very inaccurate in hard cases. Pass a callable as jac for that.

Choosing the method and the loss

Method decision table

Situation	Use	What the SciPy reference says
Default, large sparse problems, bounds	method='trf'	Default; “particularly suitable for large sparse problems with bounds”.
Small bounded problems	method='dogbox'	Suited to small bounded problems; “not recommended for problems with rank-deficient Jacobian”.
Unconstrained, at least as many residuals as parameters, plain squared error	method='lm'	MINPACK; “should be your first choice for unconstrained problems”. No bounds, no sparse Jacobians, only loss='linear', and “doesn’t work when m < n”.

Here m is the number of residuals and n the number of parameters. callback works only with trf and dogbox. For large sparse problems pass jac_sparsity, which forces tr_solver='lsmr'; the reference example uses n = 100000.

Loss functions and f_scale

loss	ρ(z)	Note from the reference
'linear' (default)	z	Standard least squares; f_scale has no effect.
'soft_l1'	2((1+z)^0.5 – 1)	The reference suggests trying it first.
'huber'	z if z ≤ 1, else 2z^0.5 – 1	The reference suggests trying it first.
'cauchy'	ln(1 + z)	The reference warns it may cause optimisation difficulties.
'arctan'	arctan(z)	Listed as a built-in option.

A callable loss is also accepted. f_scale (default 1.0) rescales the loss: the scaled loss is C² · ρ(f²/C²) with C = f_scale. The reference example uses 0.1; choose a value that matches the residual size you consider normal for your data.

Output and how to read it

least_squares returns an OptimizeResult. The fields to check first:

Field	Meaning
x	Solution parameters.
cost	0.5 * sum(ρ(f_i²)) at the solution; depends on the chosen loss.
fun	Residual vector at the solution. Inspect it for outliers; see how to calculate residuals.
jac	Modified Jacobian at the solution; JᵀJ approximates the Gauss-Newton Hessian.
optimality	First-order optimality; for unconstrained problems, the uniform norm of the gradient.
active_mask	Per parameter: 0 inactive, -1 lower bound active, 1 upper bound active.
nfev, njev	Function and Jacobian evaluation counts.
success	True when status is greater than 0.

Status codes

status	Meaning
-2	The callback raised StopIteration.
-1	Improper input to MINPACK (lm).
0	Maximum number of function evaluations exceeded.
1	gtol termination condition satisfied.
2	ftol termination condition satisfied.
3	xtol termination condition satisfied.
4	Both ftol and xtol satisfied.

The tolerances ftol, xtol and gtol each default to 1e-08, and passing None disables a criterion unless method='lm'. success=True means a tolerance was met, not that the model is right: check the residuals. least_squares does not return parameter uncertainties; if you need them from a single fit, curve_fit returns a covariance matrix directly.

Troubleshooting

status 0 and success False

The default max_nfev is 100 * n (evaluations used for finite-difference Jacobians are not counted). Raise max_nfev, improve x0, or supply an analytic jac. With max_nfev=2 on the example above the result is status 0, “The maximum number of function evaluations is exceeded.” and success=False.

ValueError with method=’lm’

These three messages are raised by SciPy, each matching a documented limit of lm:

• Method 'lm' doesn't support bounds. Use trf or dogbox.
• method='lm' supports only 'linear' loss function. Use trf for robust loss.
• Method 'lm' doesn't work when the number of residuals is less than the number of variables. Add data or fix parameters.

Poor or unstable convergence

Finite-difference Jacobians can be inaccurate: try jac='3-point', set diff_step, or pass an analytic Jacobian. If parameters differ by orders of magnitude, set x_scale (an array, or 'jac'). With dogbox, a rank-deficient Jacobian is a documented weak spot. Residual functions that return complex values must be split into real and imaginary parts, as in the reference example.

Parameters pinned at a bound

A nonzero active_mask entry tells you the optimum is on the constraint. Widen the bound or reconsider the model before trusting the value. For the norm of a residual vector, see numpy.linalg.norm.

Free alternatives

Option	Use it when
scipy.optimize.curve_fit	You have x and y data and want fitted parameters with a covariance matrix from one call.
scipy.optimize.minimize	Your objective is a single scalar rather than a residual vector.
statsmodels OLS	The model is linear in its parameters and you want regression tables.

FAQ

What is scipy least squares?

It usually means scipy.optimize.least_squares, the SciPy function for nonlinear least-squares problems with optional bounds and robust loss functions. For a straight-line fit, a linear tool such as statsmodels OLS is simpler.

What is the difference between least_squares and curve_fit?

curve_fit takes a model and x, y data and returns parameters plus a covariance matrix. least_squares takes a residual function, exposes the loss, method, Jacobian and stopping options, and returns the full OptimizeResult. See the curve_fit guide.

Which loss should I use for outliers?

The reference advises trying soft_l1 or huber first and setting f_scale to suit your data, as in the reference example.

Can I use bounds with method=’lm’?

No. lm does not support bounds; use trf (the default) or dogbox.

Sources

• SciPy v1.18.0 reference, scipy.optimize.least_squares, checked 2026-10-29.
• SciPy tutorial, Optimization (least-squares section), checked 2026-10-29.
• SciPy, Install SciPy, checked 2026-10-29.
• All printed numbers above were produced by running the shown code with SciPy 1.18.0 and NumPy 2.5.1.