SciPy Cheatsheet
Optimization
Use this SciPy reference while you build software engineering projects, review code for technical interview prep, or polish examples for a software engineer resume.
Core Import
from scipy.optimize import ( minimize, minimize_scalar, root, root_scalar, curve_fit, least_squares, linprog, milp, differential_evolution, dual_annealing, shgo, direct, brentq, bisect, newton, check_grad, approx_fprime )
minimize — Unconstrained & Constrained Optimization
from scipy.optimize import minimize import numpy as np # Rosenbrock function def rosen(x): return sum(100*(x[1:] - x[:-1]**2)**2 + (1 - x[:-1])**2) result = minimize(rosen, x0=[0, 0], method='BFGS') print(result.x) # optimal point print(result.fun) # optimal value print(result.success) # bool print(result.message) # status string print(result.nit) # iterations
Method Selection Table
| Method | Needs Gradient | Needs Hessian | Handles Constraints | Best For |
|---|---|---|---|---|
'Nelder-Mead' | No | No | No | Noisy/non-smooth, low-dim |
'BFGS' | Yes (or FD) | No (approx) | No | Smooth unconstrained |
'L-BFGS-B' | Yes (or FD) | No (approx) | Bounds only | Large-scale, bounds |
'CG' | Yes (or FD) | No | No | Smooth, memory-limited |
'Newton-CG' | Yes | Yes (or FD) | No | Smooth, exact Hessian avail |
'trust-ncg' | Yes | Yes | No | Smooth, large-scale |
'SLSQP' | Yes (or FD) | No | Bounds + ineq + eq | General constrained |
'trust-constr' | Yes | Yes (or FD) | Bounds + ineq + eq | Constrained, robust |
'cobyla' | No | No | Ineq only | Derivative-free constrained |
'TNC' | Yes (or FD) | No | Bounds only | Bounded, large-scale |
# With gradient (jac) def rosen_grad(x): g = np.zeros_like(x) g[:-1] = -400*x[:-1]*(x[1:] - x[:-1]**2) - 2*(1 - x[:-1]) g[1:] += 200*(x[1:] - x[:-1]**2) return g result = minimize(rosen, [0, 0], method='BFGS', jac=rosen_grad) # With Hessian from scipy.optimize import rosen_hess result = minimize(rosen, [0, 0], method='Newton-CG', jac=rosen_grad, hess=rosen_hess) # With bounds (L-BFGS-B, SLSQP, TNC, trust-constr) bounds = [(-2, 2), (-2, 2)] result = minimize(rosen, [0, 0], method='L-BFGS-B', bounds=bounds) # With constraints (SLSQP, trust-constr) constraints = [ {'type': 'eq', 'fun': lambda x: x[0] + x[1] - 1}, # x0+x1=1 {'type': 'ineq', 'fun': lambda x: x[0] - 0.5}, # x0>=0.5 ] result = minimize(rosen, [0.5, 0.5], method='SLSQP', bounds=bounds, constraints=constraints) # Options result = minimize(rosen, [0, 0], method='BFGS', options={'maxiter': 1000, 'gtol': 1e-8, 'disp': True})
Bounds Objects (New API)
from scipy.optimize import Bounds b = Bounds(lb=[-2, -2], ub=[2, 2], keep_feasible=True) result = minimize(rosen, [0, 0], bounds=b)
Constraint Objects (New API)
from scipy.optimize import LinearConstraint, NonlinearConstraint # A @ x in [lb, ub] lc = LinearConstraint([[1, 1]], lb=0.5, ub=1.5) # lb <= fun(x) <= ub nlc = NonlinearConstraint(lambda x: x[0]*x[1], lb=0.1, ub=np.inf) result = minimize(rosen, [0.5, 0.5], method='trust-constr', constraints=[lc, nlc])
minimize_scalar — 1-D Minimization
from scipy.optimize import minimize_scalar f = lambda x: (x - 2)**2 + 1 # Unbounded (Brent's method) result = minimize_scalar(f) print(result.x, result.fun) # 2.0, 1.0 # Bounded result = minimize_scalar(f, bounds=(0, 5), method='bounded') # With bracket result = minimize_scalar(f, bracket=(-1, 3), method='brent')
Root Finding
root_scalar — Scalar Equations
from scipy.optimize import root_scalar f = lambda x: x**3 - x - 2 # Bracketing methods (guaranteed convergence) result = root_scalar(f, bracket=[1, 2], method='brentq') result = root_scalar(f, bracket=[1, 2], method='bisect') result = root_scalar(f, bracket=[1, 2], method='ridder') # Derivative methods (faster but need good x0) result = root_scalar(f, x0=1.5, method='newton', fprime=lambda x: 3*x**2 - 1) # Secant (no derivative needed, but needs two points) result = root_scalar(f, x0=1.0, x1=2.0, method='secant') print(result.root, result.converged, result.iterations)
Standalone Bracketing Functions
from scipy.optimize import brentq, bisect, brenth, ridder, toms748 root = brentq(lambda x: x**2 - 2, 1, 2) # most robust root = bisect(lambda x: x**2 - 2, 1, 2) # slowest, most reliable root = toms748(lambda x: x**2 - 2, 1, 2) # modern, efficient # All accept: xtol, rtol, maxiter, full_output root, info = brentq(lambda x: x**2 - 2, 1, 2, full_output=True)
root — Systems of Equations
from scipy.optimize import root def equations(x): return [x[0]**2 + x[1]**2 - 4, x[0] - x[1]] result = root(equations, x0=[1, 1]) print(result.x) # [sqrt(2), sqrt(2)] # With Jacobian def jac(x): return [[2*x[0], 2*x[1]], [1, -1 ]] result = root(equations, [1, 1], jac=jac, method='hybr')
| Root Method | Notes |
|---|---|
'hybr' | Default; MINPACK hybrid (Powell's method) |
'lm' | Levenberg-Marquardt, overdetermined ok |
'broyden1' | Broyden's first method |
'broyden2' | Broyden's second method |
'krylov' | Large sparse systems |
'df-sane' | Derivative-free |
Curve Fitting
from scipy.optimize import curve_fit # Define model 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.randn(50) # Fit popt, pcov = curve_fit(model, xdata, ydata) print(popt) # best-fit parameters [a, b, c] perr = np.sqrt(np.diag(pcov)) # 1-sigma uncertainties # With bounds and initial guess popt, pcov = curve_fit(model, xdata, ydata, p0=[2, 1, 0], bounds=([0, 0, -np.inf], [10, 10, np.inf]), method='trf', # 'trf', 'dogbox', 'lm' maxfev=5000) # With sigma (weights) popt, pcov = curve_fit(model, xdata, ydata, sigma=0.2*np.ones_like(ydata), absolute_sigma=True)
Least Squares
from scipy.optimize import least_squares def residuals(params, x, y): a, b, c = params return a * np.exp(-b * x) + c - y result = least_squares(residuals, x0=[1, 1, 0], args=(xdata, ydata), method='trf', # 'trf', 'dogbox', 'lm' bounds=([0, 0, -np.inf], [10, 10, np.inf]), loss='soft_l1', # robust loss f_scale=0.1) print(result.x, result.cost, result.success)
loss option | Notes |
|---|---|
'linear' | Standard least squares (default) |
'soft_l1' | Smooth approximation of L1 |
'huber' | Huber loss (robust to outliers) |
'cauchy' | Very robust, heavy tails |
'arctan' | Bounded loss |
Linear Programming
from scipy.optimize import linprog # Minimize c @ x subject to A_ub @ x <= b_ub, A_eq @ x == b_eq, bounds c = [-1, -2] # maximize x0 + 2*x1 → minimize negation A_ub = [[1, 1], [-1, 2]] b_ub = [10, 8] bounds = [(0, None), (0, None)] # x0, x1 >= 0 result = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs') print(result.x, result.fun)
Mixed-Integer LP
from scipy.optimize import milp, LinearConstraint, Bounds c = np.array([-1.0, -2.0]) constraints = LinearConstraint([[1, 1], [-1, 2]], lb=-np.inf, ub=[10, 8]) integrality = np.array([1, 1]) # 1=integer, 0=continuous bounds = Bounds(lb=[0, 0], ub=[np.inf, np.inf]) result = milp(c, constraints=constraints, integrality=integrality, bounds=bounds)
Global Optimization
from scipy.optimize import differential_evolution, dual_annealing, shgo, direct bounds = [(-5, 5), (-5, 5)] # Differential Evolution (stochastic, parallel-ready) result = differential_evolution(rosen, bounds, seed=42, workers=1, popsize=15, mutation=(0.5, 1), recombination=0.7, maxiter=1000, tol=1e-7, polish=True) # Dual Annealing (combines SA + local search) result = dual_annealing(rosen, bounds, seed=42, maxiter=1000) # SHGO (simplicial homology, finds ALL local minima) result = shgo(rosen, bounds, n=100, iters=5) print(result.xl) # all local minima found # DIRECT (Dividing Rectangles, deterministic) result = direct(rosen, bounds, maxiter=1000)
Gradient Checking
from scipy.optimize import check_grad, approx_fprime f = lambda x: x[0]**2 + x[1]**2 g = lambda x: np.array([2*x[0], 2*x[1]]) err = check_grad(f, g, [1.0, 1.0]) print(f"Gradient error: {err:.2e}") # should be near machine epsilon # Finite-difference gradient fd_grad = approx_fprime([1.0, 1.0], f, epsilon=1e-7)
Common Gotchas
Always use
floatarrays. Integer inputs cause silent failures or unexpected downcasting.
curve_fituses Levenberg-Marquardt ('lm') by default only for unconstrained problems. Passingboundsautomatically switches the method to'trf'; explicitly combiningmethod='lm'withboundsraisesValueError('lm'cannot handle bounds).
minimizeresult is always returned even on failure. Always checkresult.successbefore usingresult.x.
Constraint type
'ineq'meansfun(x) >= 0(not<= 0). A constraintx >= 0.5is{'type': 'ineq', 'fun': lambda x: x[0] - 0.5}.