scipy 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)

• from scipy import stats
  Idiomatic. Submodules are lazily loaded, so this is cheap.
• import scipy # for scipy.io
  io collides with the stdlib module — docs explicitly prefer this form.
• from scipy.stats import norm
  Fine. Going one level deeper is only safe for documented public subpackages.

Gotchafrom 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.

• NumPy
  The ndarray + elementwise ops. No algorithms.
• SciPy
  Thin, well-tested Python over Fortran/C: LAPACK, ARPACK, QUADPACK, ODEPACK, FITPACK, SuperLU, HiGHS, Qhull.
• Above it
  scikit-learn, statsmodels, scikit-image all sit on these primitives.

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

TipUnpacking 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

• assume_a=
  'gen' | 'sym' | 'her' | 'pos' — picks the LAPACK driver.
• linalg.solve(A, b, assume_a='pos')
  The single highest-leverage keyword in the module.

Gotchalinalg.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)

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

• P, L, U = linalg.lu(A)
  Permuted LU. lu_factor is the one you solve with.
• Q, R = linalg.qr(A, mode='economic')
  Least squares, orthonormal bases, Gram–Schmidt done right.
• U, s, Vh = linalg.svd(A, full_matrices=False)
  The workhorse: rank, pinv, PCA, conditioning. s is a 1-D vector.
• L = linalg.cholesky(A, lower=True)
  SPD only — raises LinAlgError otherwise (a cheap SPD test).
• w, v = linalg.eigh(A)
  Symmetric/Hermitian — real eigenvalues, ascending, orthonormal v.
• w, v = linalg.eig(A)
  General — complex eigenvalues, no ordering guarantee.
• linalg.eigh(A, subset_by_index=[0, 4])
  Just the 5 smallest. Do not compute all n and slice.
• linalg.qz(A, B) / linalg.schur(A)
  Generalized / Schur forms.

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

• linalg.expm(A)
  True matrix exponential — solves ỹ = Ay. Not np.exp(A).
• linalg.logm / sqrtm / funm
  Inverse, square root, arbitrary f(A) via Schur–Parlett.
• linalg.fractional_matrix_power(A, 0.5)
  A**p for real p.
• linalg.norm(A, 2) / linalg.det / linalg.matrix_rank
  2-norm = largest singular value.

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

• linalg.solve_triangular(L, b, lower=True)
  O(n²). Free win after a factorization.
• linalg.solve_banded((l, u), ab, b)
  Tridiagonal ODE/PDE systems. ab is the diagonal-ordered form.
• linalg.solve_toeplitz(c_r, b)
  Convolution / AR systems, O(n²).
• linalg.block_diag(A, B, C)
  Assemble block-diagonal.
• linalg.circulant / toeplitz / hankel / hilbert
  Constructors.
• linalg.null_space(A) / orth(A)
  SVD-based bases.
• linalg.solve_sylvester / solve_lyapunov / solve_continuous_are
  Control theory (Riccati, Lyapunov).

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.

TipUnbounded 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])

GotchaEquality 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)

• loss='linear'
  Ordinary least squares (default).
• loss='soft_l1' / 'huber' / 'cauchy'
  Down-weights outliers. f_scale sets where 'outlier' begins.
• res.jac
  Use it for covariance: inv(J.T @ J) * s^2.

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

• opt.differential_evolution(f, bounds)
  Robust default for multimodal. Try workers=-1, polish=True.
• opt.shgo(f, bounds)
  Simplicial homology; finds all local minima, good with constraints.
• opt.dual_annealing(f, bounds)
  Simulated annealing + local polish.
• opt.basinhopping(f, x0)
  Monte-Carlo hops + local minimize. Good for rugged landscapes.
• opt.direct(f, bounds)
  Deterministic, Lipschitz. Reproducible.
• opt.brute(f, ranges)
  Grid search + polish. Sanity check only.

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

Tipbrentq 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.

• opt.minimize_scalar(f, bounds=(0,5), method='bounded')
  Bounded scalar — use this, not minimize with a 1-vector.
• opt.minimize_scalar(f, bracket=(a,b,c))
  Brent, unbounded.
• opt.curve_fit / opt.nnls(A, b)
  Non-negative least squares.

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

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

• ig.trapezoid(y, x)
  Robust default for arbitrary spacing.
• ig.simpson(y, x)
  Higher order; wants an odd number of points.
• ig.cumulative_trapezoid(y, x, initial=0)
  Running integral — e.g. velocity → position.
• ig.cumulative_simpson(y, x=x)
  Higher-order cumulative.

Gotchatrapz, 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

• method=
  RK45 (default) · DOP853 · Radau / BDF (stiff) · LSODA (auto)
• jac= / jac_sparsity=
  Huge speedups for stiff or large systems.
• first_step=, max_step=
  Only reach for these when the solver stalls.

GotchaDefaults 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)

TipFeed 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)

• bc_type=
  'not-a-knot' (default) · 'natural' · 'clamped' · 'periodic'
• extrapolate=False
  Return NaN outside the data range instead of inventing values.

Gotchainterp1d 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)

Gotchainterp2d 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

• loc / scale
  Universal shift & stretch. expon(scale=1/lam) — SciPy uses scale, not rate.
• sf / isf
  Use instead of 1 - cdf: survives underflow past ~1e-16.
• shape params
  Come first: stats.gamma(a, loc=, scale=).

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

Tipstats.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.describe(x)
  nobs, minmax, mean, variance, skewness, kurtosis in one shot.
• stats.zscore(x, axis=0, ddof=1)
  Standardize.
• stats.iqr / trim_mean(x, 0.1) / sem
  Robust spread, trimmed mean, standard error.
• stats.gmean / hmean / pmean
  Geometric / harmonic / power means.
• stats.skew / kurtosis(x, fisher=True)
  Fisher → excess kurtosis (normal = 0).
• stats.mode(x, keepdims=False)
  Modal value + count.
• stats.rankdata(x, method='average')
  Ranks — the basis of every nonparametric test.
• stats.gaussian_kde(x)(grid)
  KDE. bw_method='silverman' or a float.
• stats.ecdf(x).cdf.evaluate(t)
  Empirical CDF with confidence bands.

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

• alternative=
  'two-sided' (default) · 'less' · 'greater'
• res.confidence_interval()
  Available on t-tests and several others — report it, not just p.

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

• stats.pearsonr(x, y)
  Linear. Assumes normality for its p-value and CI.
• stats.spearmanr(x, y)
  Monotone, rank-based. Handles nonlinearity and outliers.
• stats.kendalltau(x, y)
  Concordance. Better for small n / many ties.
• stats.linregress(x, y)
  Simple OLS: .slope .intercept .rvalue .pvalue .stderr.
• stats.theilslopes / siegelslopes
  Robust regression slopes — outlier-proof.
• stats.pointbiserialr / chisquare
  Binary–continuous / goodness of fit.

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

TipWhen 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

GotchaSobol 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

• workers=-1
  Parallelize over the leading axes across all cores. Free speedup.
• fft.next_fast_len(n)
  Pad to a 5-smooth length → can be 10× faster than a prime n.
• norm=
  'backward' (default, 1/n on inverse) · 'ortho' · 'forward'

Gotchascipy.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!

• fs=
  Pass it and give Wn in Hz. Otherwise Wn is normalized to Nyquist (0–1).
• output='sos'
  'ba' (transfer-function) coefficients go numerically unstable above ~order 8.
• IIR vs FIR
  IIR: cheap, nonlinear phase. FIR: expensive, exactly linear phase.

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

• sg.sosfiltfilt(sos, x)
  Forward + backward → zero phase, doubled order. Offline analysis.
• sg.sosfilt(sos, x, zi=zi)
  Causal, streaming, keeps state. Real-time.
• sg.filtfilt(b, a, x) / sg.lfilter(b, a, x)
  Same pair, transfer-function form.
• sg.savgol_filter(x, 51, 3)
  Polynomial smoothing — preserves peak height/width.
• sg.medfilt(x, 5) / sg.wiener(x)
  Impulse-noise removal / adaptive denoise.
• sg.detrend(x, type='linear')
  Kill the trend before you FFT.
• sg.hilbert(x) / sg.envelope(x)
  Analytic signal → instantaneous amplitude & phase.

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

• sg.fftconvolve(a, v, mode='same')
  O(n log n). Fastest for long kernels.
• sg.oaconvolve(a, v)
  Overlap-add — long signal, short kernel.
• sg.convolve(a, v, method='auto')
  Picks direct vs FFT for you.
• sg.choose_conv_method(a, v)
  Tells you what 'auto' would choose.
• sg.correlate(a, b, mode='full')
  Cross-correlation.
• sg.correlation_lags(len(a), len(b), 'full')
  The lag axis. Do not hand-roll it.
• sg.convolve2d / correlate2d
  2-D (images).

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

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

• sg.resample_poly(x, up, down)
  Rational rate change, polyphase. The default choice.
• sg.decimate(x, q)
  Downsample with an anti-alias filter built in.
• sg.resample(x, num)
  FFT-based — assumes the signal is periodic. Edge artifacts otherwise.
• sys = sg.lti(num, den) / sg.TransferFunction / ZerosPolesGain / StateSpace
  LTI models.
• sg.bode(sys) / sg.step(sys) / sg.impulse(sys) / sg.lsim(sys, u, t)
  Responses.
• sg.tf2zpk / zpk2sos / tf2ss / bilinear
  Representation conversions.

GotchaNever 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

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

GotchaMixing 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)

Gotchainfo != 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)

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

• distance.pdist(X, 'euclidean')
  Condensed — all pairs within X, as a flat vector of length n(n-1)/2.
• distance.cdist(XA, XB, 'cosine')
  Full m×n matrix between two sets.
• distance.squareform(d)
  Condensed ↔ square. Round-trips both ways.
• metrics
  euclidean, cityblock, chebyshev, minkowski, cosine, correlation, mahalanobis, jaccard, hamming, jensenshannon, canberra

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

• hull = ConvexHull(pts)
  .vertices, .simplices, .volume, .area.
• tri = Delaunay(pts)
  .simplices, .find_simplex(p) → point location, -1 if outside.
• vor = Voronoi(pts)
  .vertices, .regions, .ridge_vertices. -1 = a ridge running to infinity.
• SphericalVoronoi(pts, radius, center)
  On a sphere (geodesy, molecular surfaces).
• procrustes(A, B)
  Optimal similarity alignment + disparity.
• distance_matrix(A, B, p=2)
  Dense, no tree.

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

GotchaQuaternion 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)

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

• nd.gaussian_filter(img, sigma=2)
  Separable, the standard blur. sigma can be per-axis.
• nd.uniform_filter / median_filter(img, size=3)
  Box mean / median (salt-and-pepper noise).
• nd.sobel(img, axis=0) / prewitt / laplace
  Edges.
• nd.gaussian_laplace / gaussian_gradient_magnitude
  Blob and edge strength.
• nd.convolve(img, kernel) / correlate
  Custom kernels.
• nd.generic_filter(img, np.std, size=3)
  Any Python callable on a sliding window. Slow but general.
• mode=
  'reflect' 'constant' 'nearest' 'mirror' 'wrap' — boundary handling.

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

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

• nd.zoom(img, 2.0, order=3)
  Resample. order = spline order 0–5.
• nd.rotate(img, 30, reshape=False)
  Rotate about the centre.
• nd.shift(img, (dy, dx))
  Sub-pixel translation.
• nd.affine_transform(img, M, offset=o)
  General affine. M maps output → input.
• nd.map_coordinates(img, coords, order=3)
  Arbitrary warping / interpolated sampling.

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

• special.gamma / gammaln / beta / betaln
  Log versions avoid overflow — use them.
• special.erf / erfc / erfinv
  Error function family.
• special.expit / logit / softmax / log_softmax
  Sigmoid and friends, numerically stable.
• special.logsumexp(a, b=w, axis=0)
  The stable log-sum-exp. Never write it yourself.
• special.comb(n, k, exact=True) / perm / factorial
  Combinatorics.
• special.jv / yn / iv / kn / airy / legendre_p
  Bessel, Airy, orthogonal polynomials.
• special.xlogy(x, y) / rel_entr / kl_div
  0·log0 = 0 done right — entropy and KL.

Tiplogsumexp, 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.io.loadmat(f) / savemat(f, dict)
  MATLAB .mat ≤ v7.2. For v7.3 use h5py.
• scipy.io.wavfile.read/write(f)
  WAV → (rate, ndarray).
• scipy.io.mmread / mmwrite
  Matrix Market — the interchange format for sparse matrices.
• scipy.io.arff.loadarff(f)
  Weka ARFF.
• scipy.io.hb_read / hb_write
  Harwell–Boeing sparse.
• scipy.datasets.face() / ascent() / electrocardiogram()
  Sample data (downloads on first use via pooch).

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

Gotchaseed= 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.

• workers=-1
  fft.*, KDTree.query, differential_evolution, cdist-adjacent code.
• linalg.solve(..., assume_a='pos')
  Cholesky instead of LU.
• lu_factor / cho_factor / splu / factorized
  Amortize the O(n³) across many right-hand sides.
• output='sos'
  Stable and fast for high-order filters.
• fft.next_fast_len(n)
  Pad away a prime-length FFT.
• jac= / jac_sparsity=
  Turns an intractable stiff ODE or optimization into a fast one.
• LinearOperator
  Never materialize A — supply matvec and let cg/eigsh run.
• Build sparse in COO, compute in CSR
  Avoids O(nnz) reallocation per insert.

TipBefore 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).

• dblquad(f, ...)
  The integrand is f(y, x) — inner variable first.
• solve_ivp(rhs, ...) vs odeint(rhs, ...)
  rhs(t, y) vs rhs(y, t). Opposite.
• Rotation.from_quat
  (x, y, z, w) — scalar last.
• ndimage.affine_transform(M)
  M maps output → input (inverse map).
• stats.expon(scale=1/lam)
  SciPy parameterizes by scale, never rate.
• stats.gamma(a, loc, scale)
  Shape params come before loc/scale.
• signal Wn without fs=
  Wn is then normalized to Nyquist (0–1), not Hz.

• ttest_ind default
  equal_var=True. Pass equal_var=False (Welch).
• curve_fit without p0
  Starts at all-ones; converges to nonsense on exponentials.
• dist.fit(data)
  Fits loc too. Pass floc=0 for gamma/lognorm/weibull.
• cg / gmres return value
  info != 0 = did not converge. It does not raise.
• 1 - cdf(x) in the tail
  Underflows to 0. Use sf(x).
• np.exp(A) for a matrix
  Elementwise. You wanted linalg.expm(A).
• A * B on a sparse matrix
  Matrix product for spmatrix, elementwise for sparray.
• x[::4] to downsample
  Aliases. Use signal.decimate.
• solve_ivp default tolerances
  rtol=1e-3 is loose. Tighten before trusting.
• KDTree missing neighbours
  d = inf, i = tree.n — not an error.