SciPy — the visual cheat sheet

import numpy as np
from scipy import optimize, stats, signal        # submodule as namespace
import scipy                                     # then scipy.io.loadmat(...)

scipy.__version__                                # 1.18.x  (Python 3.12-3.14, NumPy >= 2.0)

GOTCHA from scipy import * imports nothing useful — there is no flat namespace. And scipy.optimize is not auto-imported by import scipy in old code paths; import the submodule.

Rule of thumb: if it is a loop over elements, NumPy has it. If it is a numerical method with a name (Brent, BFGS, Radau, Welch, Dijkstra), SciPy has it.

res = optimize.minimize(f, x0)
res.x, res.fun, res.success, res.nit, res.message   # not a tuple

r = stats.ttest_ind(a, b)
r.statistic, r.pvalue, r.confidence_interval()     # named, tab-completable

s = integrate.solve_ivp(rhs, [0, 10], y0, dense_output=True)
s.t, s.y, s.sol(3.7), s.t_events

TIP Unpacking legacy tuples still works for old APIs, but reach for the attribute names — they survive version bumps.

from scipy import linalg

x = linalg.solve(A, b)                     # general square
x = linalg.solve(A, b, assume_a='pos')     # SPD  -> Cholesky, ~2x faster
x = linalg.solve(A, b, assume_a='sym')     # symmetric indefinite
x = linalg.lstsq(A, b)[0]                  # over/under-determined, min-norm
x = linalg.pinv(A) @ b                     # same, but slower - prefer lstsq

GOTCHA linalg.inv(A) @ b is slower and less accurate than solve. Inverses are for textbooks.

lu, piv = linalg.lu_factor(A)              # O(n^3) once
x1 = linalg.lu_solve((lu, piv), b1)        # O(n^2) each
x2 = linalg.lu_solve((lu, piv), b2)

c, low = linalg.cho_factor(A)              # SPD variant
x  = linalg.cho_solve((c, low), b)

TIP Many right-hand sides? Also just stack them: linalg.solve(A, B) with B.shape == (n, k).

GOTCHA eig on a symmetric matrix wastes time and hands you complex dtype noise. Use eigh.

GOTCHA np.exp(A) exponentiates each entry. That is a different, usually wrong, matrix.

Batched (..., n, n) stacks? Use np.linalg. Everything else: scipy.linalg.

from scipy import optimize as opt

res = opt.minimize(f, x0, args=(a, b),
                   method='L-BFGS-B',
                   jac=grad_f,                  # analytic gradient -> big speedup
                   bounds=opt.Bounds([0, 0], [np.inf, 5]),
                   options={'maxiter': 500, 'ftol': 1e-10})

res.x, res.fun, res.success, res.nit, res.message

No jac? SciPy finite-differences it (n+1 extra f-evals per step). Supply the gradient whenever you can, or get it from scipy.differentiate.jacobian.

TIP Unbounded smooth problem? 'BFGS'. Add bounds → 'L-BFGS-B'. Add real constraints → 'trust-constr'.

lin = opt.LinearConstraint([[1, 1], [1, -1]], lb=[-np.inf, 0], ub=[3, np.inf])
non = opt.NonlinearConstraint(lambda x: x[0]**2 + x[1]**2, 0, 1, jac=cjac)

res = opt.minimize(f, x0, method='trust-constr',
                   bounds=opt.Bounds(lb, ub),
                   constraints=[lin, non])

GOTCHA Equality constraint? Set lb == ub. The old {'type':'eq','fun':...} dict form only works for SLSQP/COBYLA.

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

p, cov = opt.curve_fit(model, xdata, ydata,
                       p0=[1, 1, 0],            # ALWAYS give p0
                       sigma=yerr, absolute_sigma=True,
                       bounds=([0, 0, -1], [10, 5, 1]))
perr = np.sqrt(np.diag(cov))                    # 1-sigma parameter errors

# lower level, more control, robust losses:
res = opt.least_squares(residuals, x0, loss='soft_l1', f_scale=0.1)

GOTCHA curve_fit without p0 starts at all-ones and silently converges to garbage on exponentials and power laws.

All take bounds as a list of (lo, hi). Seeding is now rng= (SPEC-7); seed= still works but is legacy.

r = opt.root_scalar(f, bracket=[0, 2], method='brentq')      # 1-D, sign change known
r = opt.root_scalar(f, x0=1.0, fprime=df, method='newton')   # 1-D, derivative known
r.root, r.converged, r.iterations

sol = opt.root(F, x0, jac=J, method='hybr')                  # system F(x) = 0
sol.x, sol.success

from scipy.optimize import elementwise                       # vectorized, 1.15+
elementwise.find_root(lambda x: x**2 - c, (lo, hi))          # solves a whole array at once

TIP brentq on a bracket is guaranteed to converge. newton without a bracket can fly off to infinity — bracket first with elementwise.bracket_root.

res = opt.linprog(c, A_ub=A, b_ub=b, A_eq=Ae, b_eq=be,
                  bounds=(0, None), method='highs')     # HiGHS: fast, default

res = opt.milp(c=c, constraints=opt.LinearConstraint(A, lb, ub),
               integrality=[1, 1, 0])                   # 1 = integer, 0 = continuous

row, col = opt.linear_sum_assignment(cost)              # Hungarian, O(n^3)
cost[row, col].sum()

linprog minimizes c @ x. To maximize, negate c.

from scipy import integrate as ig

val, err = ig.quad(f, 0, 1, args=(k,))          # adaptive, QUADPACK
val, err = ig.quad(f, -np.inf, np.inf)          # infinite limits are fine
val, err = ig.quad(f, 0, 1, points=[0.5])       # declare known singularities

ig.dblquad(f, xa, xb, ya, yb)                   # note: f(y, x) - y FIRST
ig.tplquad(f, ...) / ig.nquad(f, ranges)        # 3-D / n-D

ig.cubature(f, [0,0], [1,1])                    # 1.15+: vectorized n-D, preferred
ig.tanhsinh(f, a, b)                            # 1.15+: endpoint singularities
ig.nsum(lambda k: 1/k**2, 1, np.inf)            # 1.15+: infinite series

GOTCHA dblquad's integrand signature is f(y, x) — inner variable first. It is the most-reported SciPy footgun.

GOTCHA trapz, cumtrapz, simps were removed in 1.14. Use the long names.

def rhs(t, y, k):                                # signature: (t, y, *args)
    return [y[1], -k * y[0]]

sol = ig.solve_ivp(rhs, t_span=(0, 20), y0=[1, 0], args=(2.0,),
                   method='RK45',                # LSODA if unsure, Radau/BDF if stiff
                   t_eval=np.linspace(0, 20, 400),
                   dense_output=True,
                   rtol=1e-8, atol=1e-10)

sol.t, sol.y          # y has shape (n_states, n_times)
sol.sol(3.7)          # continuous interpolant, any t

GOTCHA Defaults are rtol=1e-3, atol=1e-6 — loose. Tighten them before you trust a trajectory.

def hits_ground(t, y):  return y[0]              # root of this = event
hits_ground.terminal  = True                     # stop integrating
hits_ground.direction = -1                       # only downward crossings

sol = ig.solve_ivp(rhs, (0, 100), y0, events=hits_ground)
sol.t_events[0], sol.y_events[0]

# two-point BVP
sol = ig.solve_bvp(fun, bc, x_mesh, y_guess)     # sol.sol(x) -> solution

odeint still works but is the legacy LSODA wrapper with a (y, t) argument order. New code: solve_ivp.

from scipy import differentiate as df

df.derivative(np.sin, 1.0)          # .df, .error - adaptive-step, high order
df.jacobian(F, x)                   # vector-in, vector-out  -> .df is (m, n)
df.hessian(f, x)                    # vector-in, scalar-out  -> .df is (n, n)

TIP Feed jacobian(...).df straight into minimize(jac=...) when you have no analytic gradient but want better-than-default finite differences.

from scipy import interpolate as ip

np.interp(xnew, x, y)                       # linear, no object, fastest
spl = ip.make_interp_spline(x, y, k=3)      # cubic B-spline - the recommended default
cs  = ip.CubicSpline(x, y, bc_type='natural')
pch = ip.PchipInterpolator(x, y)            # monotone data -> no overshoot
ak  = ip.Akima1DInterpolator(x, y)          # tames wiggles near outliers

spl(xnew); spl.derivative()(xnew); spl.antiderivative()(xnew)
cs.roots(); cs.integrate(a, b)

GOTCHA interp1d is legacy — the docs recommend np.interp (linear) or make_interp_spline (spline) in new code. It still works; it just won't improve.

# data ON A REGULAR GRID: pass the AXES, not a meshgrid
rgi = ip.RegularGridInterpolator((xs, ys, zs), values, method='cubic',
                                 bounds_error=False, fill_value=None)   # None -> extrapolate
rgi(np.array([[0.1, 2.3, 4.5]]))

# SCATTERED points
ip.griddata(points, values, xi, method='cubic')      # one-shot convenience
lin  = ip.LinearNDInterpolator(points, values)       # reusable object
ct   = ip.CloughTocher2DInterpolator(points, values) # C1-smooth, 2-D
rbf  = ip.RBFInterpolator(points, values, kernel='thin_plate_spline',
                          smoothing=0.0, neighbors=64)

GOTCHA interp2d was removed. Regular grid → RegularGridInterpolator; scattered → LinearNDInterpolator / RBFInterpolator.

ip.make_smoothing_spline(x, y)              # GCV picks lambda for you  <- start here
ip.make_splrep(x, y, s=0.5)                 # explicit smoothing factor s
ip.UnivariateSpline(x, y, s=len(y))         # OO wrapper, s ~ n*sigma^2
ip.AAA(x, y)                                # 1.15+: barycentric rational, poles & all

s=0 means interpolate exactly. Increasing s trades fidelity for smoothness. splrep/splev are the older procedural FITPACK interface.

from scipy import stats

rv = stats.norm(loc=mu, scale=sd)        # FROZEN: params bound once
rv.pdf(x)    rv.cdf(x)    rv.ppf(0.975)  # density | P(X<=x) | quantile
rv.sf(x)     rv.isf(0.05)                # 1-cdf (accurate in the far tail!) | its inverse
rv.rvs(size=1000, random_state=rng)      # samples
rv.mean(), rv.var(), rv.std(), rv.interval(0.95)
rv.stats(moments='mvsk')                 # mean, var, skew, kurtosis

stats.gamma.fit(data)                    # MLE -> (a, loc, scale)
stats.gamma.fit(data, floc=0)            # pin location: almost always what you want

GOTCHA fit() estimates loc too, which almost never makes sense for gamma/lognorm/weibull. Pass floc=0.

TIP stats.fit(dist, data, bounds=...) is the modern MLE front-end and returns a result object with .plot().

X = stats.Normal(mu=0.0, sigma=1.0)       # a real object, not a frozen tuple
X.pdf(0.0); X.mean(); X.sample(shape=5)

Y = stats.make_distribution(stats.gamma)(a=2.0)     # upgrade any old distribution
Z = stats.Mixture([stats.Normal(), stats.Normal(mu=3)], weights=[0.4, 0.6])

stats.abs(X); stats.exp(X); stats.log(X)  # folded normal, lognormal, ... by transformation
3 * X + 1                                 # arithmetic on random variables

Faster, more accurate, and composable. The classic rv_continuous API is not going away — but new numerical work is landing here.

stats.ttest_ind(a, b, equal_var=False)      # Welch - the safe default for 2 groups
stats.ttest_rel(before, after)              # paired
stats.ttest_1samp(x, popmean=0)

stats.mannwhitneyu(a, b, alternative='greater')   # rank-based counterparts
stats.wilcoxon(before, after)
stats.kruskal(g1, g2, g3)

stats.f_oneway(g1, g2, g3)                  # one-way ANOVA
stats.tukey_hsd(g1, g2, g3)                 # post-hoc, all pairs
stats.chi2_contingency(table)               # -> chi2, p, dof, expected
stats.fisher_exact(table_2x2)
stats.shapiro(x); stats.anderson(x); stats.ks_2samp(a, b)
stats.levene(a, b)                          # variance equality, robust

GOTCHA ttest_ind defaults to equal_var=True (Student's). Real data rarely obliges; pass equal_var=False.

Multiple regression is out of scope — that is statsmodels or np.linalg.lstsq.

res = stats.bootstrap((x,), np.median, confidence_level=0.95,
                      n_resamples=9999, method='BCa', rng=0)
res.confidence_interval.low, res.standard_error

def diff(a, b, axis):  return np.mean(a, axis=axis) - np.mean(b, axis=axis)
res = stats.permutation_test((a, b), diff, n_resamples=9999,
                             alternative='two-sided', rng=0)
res.pvalue                                   # exact-ish, distribution-free

stats.monte_carlo_test(x, stats.norm.rvs, stats.skew)   # calibrate any statistic
stats.false_discovery_control(pvals, method='bh')       # Benjamini-Hochberg

TIP When the test you want does not exist, permutation_test almost always does. Vectorize your statistic over axis and it stays fast.

from scipy.stats import qmc

sob = qmc.Sobol(d=3, scramble=True, rng=0)
u = sob.random_base2(m=10)                   # 1024 points - powers of 2 only!
x = qmc.scale(u, l_bounds=[0, -1, 5], u_bounds=[1, 1, 10])

qmc.LatinHypercube(d=3, rng=0).random(64)
qmc.Halton(d=3, rng=0).random(64)
qmc.discrepancy(u)                           # lower = better coverage

GOTCHA Sobol loses its balance properties unless n is a power of 2 — use random_base2, and it will warn you if you don't.

from scipy import fft

X = fft.rfft(x)                       # REAL input -> half spectrum. Use this.
f = fft.rfftfreq(len(x), d=1/fs)      # matching frequency axis, in Hz
amp = 2 * np.abs(X) / len(x)          # single-sided amplitude
phase = np.angle(X)
x = fft.irfft(X, n=len(x))

fft.fft / ifft / fft2 / fftn          # complex, n-D
fft.fftshift(fft.fftfreq(n, 1/fs))    # centre the axis for plotting
fft.dct(x, type=2, norm='ortho')      # cosine (JPEG, MFCC) / dst / dctn
fft.fht(a, dln, mu)                   # fast Hankel transform

GOTCHA scipy.fftpack is the legacy interface. scipy.fft supersedes it — different norm semantics, better performance.

from scipy import signal as sg

sos = sg.butter(4, 30, btype='lowpass', fs=1000, output='sos')   # ALWAYS output='sos'
sos = sg.butter(4, [20, 40], btype='bandpass', fs=1000, output='sos')
sos = sg.cheby1(4, rp=1, Wn=30, fs=fs, output='sos')             # ripple in passband
sos = sg.ellip(4, 1, 40, 30, fs=fs, output='sos')                # steepest, most ripple
sos = sg.bessel(4, 30, fs=fs, output='sos')                      # flattest group delay

N, Wn = sg.buttord(wp=30, ws=50, gpass=3, gstop=40, fs=fs)       # let it pick the order
taps  = sg.firwin(101, cutoff=30, fs=fs, window='hamming')       # FIR: linear phase
w, h  = sg.sosfreqz(sos, fs=fs)                                  # verify the response!

GOTCHA A 12th-order 'ba' Butterworth will produce garbage. 'sos' costs you nothing.

GOTCHA filtfilt is non-causal — it uses the future. Never in real time, never for latency claims.

f, Pxx = sg.welch(x, fs, nperseg=1024, noverlap=512, window='hann')  # PSD, averaged
f, Pxx = sg.periodogram(x, fs)                                       # raw, noisy

SFT = sg.ShortTimeFFT(sg.windows.hann(256, sym=True), hop=64, fs=fs, scale_to='psd')
S = SFT.stft(x)                     # 1.12+ replacement for spectrogram/stft
t = SFT.t(len(x)); f = SFT.f

f, Cxy = sg.coherence(x, y, fs)     # 0..1 per frequency
f, Pxy = sg.csd(x, y, fs)
sg.lombscargle(t, y, freqs)         # UNEVENLY sampled data

welch trades frequency resolution for variance. Bigger nperseg → sharper peaks, noisier estimate.

peaks, props = sg.find_peaks(x,
                             height=0.5,        # absolute threshold
                             prominence=0.2,    # how much it stands out  <- the good one
                             distance=20,       # min samples between peaks
                             width=(3, 50))
props['prominences'], props['widths'], props['peak_heights']

sg.peak_widths(x, peaks, rel_height=0.5)        # FWHM

TIP prominence is nearly always the right knob. height alone drowns in noise riding on a trend — detrend first.

GOTCHA Never plain-slice to downsample (x[::4]) — you alias. Use decimate.

Convert with .tocsr() / .tocsc() / .toarray(). A.nnz, A.shape, A.data, A.indices, A.indptr.

from scipy import sparse

A = sparse.coo_array((vals, (rows, cols)), shape=(m, n)).tocsr()
A = sparse.csr_array(dense)
sparse.eye_array(n, k=0)
sparse.diags_array([lo, mid, hi], offsets=[-1, 0, 1])       # tridiagonal
sparse.block_array([[A, B], [None, C]])                     # block assembly
sparse.hstack([A, B]); sparse.vstack([A, B]); sparse.kron(A, B)
sparse.random_array((m, n), density=0.01, rng=0)

sparse.issparse(A) and A.format == 'csr'                    # the modern type check

GOTCHA Use *_array, not *_matrix. The sparse matrix API is on its way out, and it silently changes what * and ** mean.

GOTCHA Mixing the two propagates the left operand's type. Pick sparray and stay there. np.func(A) on a sparse array is also wrong — NumPy treats it as an opaque object.

from scipy.sparse import linalg as spla

x = spla.spsolve(A.tocsc(), b)              # direct (SuperLU). CSC, please.
lu = spla.splu(A.tocsc()); x = lu.solve(b)  # factor once, many rhs
solve = spla.factorized(A.tocsc())          # same, functional style

x, info = spla.cg(A, b, rtol=1e-8, M=M)     # iterative: SPD
x, info = spla.gmres(A, b, M=M)             # iterative: general
M = spla.LinearOperator((n, n), matvec=lambda v: ilu.solve(v))   # preconditioner
ilu = spla.spilu(A.tocsc(), drop_tol=1e-4)

w, v = spla.eigsh(A, k=6, which='SM')       # few eigenpairs (ARPACK), symmetric
u, s, vt = spla.svds(A, k=10)               # truncated SVD
spla.norm(A); spla.expm(A); spla.matrix_power(A, 3)

GOTCHA info != 0 from cg/gmres means it did not converge. Check it. Without a preconditioner, iterative solvers on ill-conditioned systems just stop.

from scipy.sparse import csgraph

dist, pred = csgraph.dijkstra(G, indices=src, return_predecessors=True)
csgraph.shortest_path(G, method='D', directed=True)      # 'D' Dijkstra 'BF' Bellman-Ford 'FW'
n, labels = csgraph.connected_components(G, directed=False)
csgraph.minimum_spanning_tree(G)
csgraph.breadth_first_order(G, i_start)
csgraph.laplacian(G, normed=True)                        # spectral clustering input
csgraph.maximum_bipartite_matching(G)

The adjacency matrix is the graph: entry (i,j) = edge weight, 0 = no edge. Careful with genuine zero-weight edges — use csgraph_from_dense(..., null_value=...).

from scipy.spatial import KDTree

tree = KDTree(points)                       # build: O(n log n)
d, i = tree.query(xq, k=5, workers=-1)      # k nearest      (workers=-1 -> all cores)
idx   = tree.query_ball_point(xq, r=0.3)    # all within radius r
pairs = tree.query_pairs(r=0.1)             # all close pairs in the set
M = tree.sparse_distance_matrix(tree, 0.3)  # sparse, only pairs under the cutoff
d, i = tree.query(xq, k=1, distance_upper_bound=0.5)

GOTCHA Missing neighbours come back as d = inf and i = tree.n (an out-of-range index), not an exception. Mask on np.isinf(d).

GOTCHA pdist is O(n²) in memory. Above ~10k points, use a KDTree with a radius instead of materializing the matrix.

from scipy.spatial.transform import Rotation as R, Slerp

r = R.from_euler('zyx', [90, 0, 45], degrees=True)
r = R.from_quat([x, y, z, w])            # SCALAR-LAST by default
r = R.from_matrix(M); R.from_rotvec(v)

r.as_quat(); r.as_matrix(); r.as_euler('xyz', degrees=True)
r.apply(vectors)                         # rotate points
(r2 * r1).apply(v)                       # compose: r1 THEN r2
r.inv(); r.magnitude()

key = R.from_quat([[0,0,0,1], [0,0,0.707,0.707]])
Slerp([0, 1], key)([0.25, 0.5, 0.75])    # smooth interpolation between orientations
R.random(10, rng=0); R.mean()            # uniform random / average rotation

GOTCHA Quaternion order is (x, y, z, w) — scalar last. ROS, Eigen and most textbooks put w first. Pass scalar_first=True if you need that.

from scipy.cluster import hierarchy as hc, vq

Z = hc.linkage(X, method='ward', metric='euclidean')   # 'single' 'complete' 'average' 'ward'
labels = hc.fcluster(Z, t=3, criterion='maxclust')     # or criterion='distance'
hc.dendrogram(Z, truncate_mode='lastp', p=30)          # needs matplotlib
c, _ = hc.cophenet(Z, distance.pdist(X))               # how faithful is the tree? (want > 0.75)

whitened = vq.whiten(X)                                # scale each feature to unit variance
centroids, labels = vq.kmeans2(whitened, k=3, minit='++', seed=0)

GOTCHA ward is only defined for Euclidean distance. And k-means without whiten lets whichever feature has the biggest units dominate.

from scipy import ndimage as nd

mask = img > threshold
mask = nd.binary_opening(mask, structure=np.ones((3, 3)))    # erode then dilate: despeckle
mask = nd.binary_closing(mask)                               # fill pinholes
mask = nd.binary_fill_holes(mask)

lbl, n = nd.label(mask)                       # connected components
slices = nd.find_objects(lbl)                 # bounding box per label -> img[slices[0]]
sizes  = nd.sum_labels(mask, lbl, range(1, n+1))
com    = nd.center_of_mass(img, lbl, range(1, n+1))
nd.distance_transform_edt(mask)               # distance to nearest background pixel

TIP label + find_objects + sum_labels is a complete particle-analysis pipeline in three lines.

GOTCHA affine_transform's matrix is the inverse map — output coords to input coords. Half of all image warps come out mirrored because of this.

TIP logsumexp, expit, xlogy and gammaln are why scipy.special shows up inside so much ML code.

from scipy import constants as c

c.c, c.g, c.h, c.k, c.N_A, c.R, c.e        # SI values
c.pi, c.golden

c.value('Planck constant')                 # CODATA lookup
c.unit('Planck constant'); c.precision(...)
c.find('electron')                         # search the CODATA table

c.mile, c.hour, c.psi, c.bar, c.eV         # conversion factors to SI
c.convert_temperature(100, 'Celsius', 'Kelvin')

scipy.misc is gone. Its datasets moved to scipy.datasets.

from scipy import odr

model = odr.Model(lambda B, x: B[0] * x + B[1])
data  = odr.RealData(x, y, sx=x_err, sy=y_err)      # errors in BOTH variables
out   = odr.ODR(data, model, beta0=[1., 0.]).run()
out.beta, out.sd_beta, out.res_var

Use ODR instead of curve_fit when your x values have uncertainty too. curve_fit assumes x is exact.

rng = np.random.default_rng(42)              # the modern Generator

stats.bootstrap((x,), np.mean, rng=rng)
stats.qmc.Sobol(d=2, rng=rng)
optimize.differential_evolution(f, bounds, rng=rng)
sparse.random_array((m, n), density=0.01, rng=rng)

stats.norm.rvs(size=10, random_state=rng)    # rv_continuous still says random_state

GOTCHA seed= and random_state= still work but are legacy, and passing an int to rng= vs seed= can produce a different stream. Pass a Generator, not an int, and the ambiguity disappears.

export SCIPY_ARRAY_API=1        # opt in (env var, before import)

# then many scipy functions accept and RETURN the same array type:
#   PyTorch (incl. CUDA)  |  JAX (incl. jit)  |  CuPy  |  Dask  |  array-api-strict

Coverage is broad in stats, fft, signal, special, cluster, and growing each release. 1.18 added lazy-array and JAX-JIT support to a large batch of stats functions. NumPy stays the default and always works.

TIP Before optimizing: check you are not calling a Python-level function inside quad/minimize a million times. Vectorize the integrand instead (cubature, tanhsinh, and elementwise.* are built for this).