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

MethodNeeds GradientNeeds HessianHandles ConstraintsBest For
'Nelder-Mead'NoNoNoNoisy/non-smooth, low-dim
'BFGS'Yes (or FD)No (approx)NoSmooth unconstrained
'L-BFGS-B'Yes (or FD)No (approx)Bounds onlyLarge-scale, bounds
'CG'Yes (or FD)NoNoSmooth, memory-limited
'Newton-CG'YesYes (or FD)NoSmooth, exact Hessian avail
'trust-ncg'YesYesNoSmooth, large-scale
'SLSQP'Yes (or FD)NoBounds + ineq + eqGeneral constrained
'trust-constr'YesYes (or FD)Bounds + ineq + eqConstrained, robust
'cobyla'NoNoIneq onlyDerivative-free constrained
'TNC'Yes (or FD)NoBounds onlyBounded, 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 MethodNotes
'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 optionNotes
'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 float arrays. Integer inputs cause silent failures or unexpected downcasting.

curve_fit uses Levenberg-Marquardt ('lm') by default only for unconstrained problems. Passing bounds automatically switches the method to 'trf'; explicitly combining method='lm' with bounds raises ValueError ('lm' cannot handle bounds).

minimize result is always returned even on failure. Always check result.success before using result.x.

Constraint type 'ineq' means fun(x) >= 0 (not <= 0). A constraint x >= 0.5 is {'type': 'ineq', 'fun': lambda x: x[0] - 0.5}.