SciPy cheat sheet

A scannable SciPy reference: 15 short snippets across 11 topics, each linking back to the lesson it came from.

At a glance

Quick snippets

SciPy arrays and NumPy interop

SciPy is NumPy plus algorithms

import numpy as np
import scipy.linalg as la

a = np.array([[3.0, 1.0], [1.0, 2.0]])
b = np.array([9.0, 8.0])

x = la.solve(a, b)              # SciPy, not numpy.linalg
print(x, x.dtype, x.shape)      # [3. 2.] float64 (2,)

print(la.det(a), la.eigvals(a))
print(np.allclose(a @ x, b))

Sparse matrices: the array you never materialise

import numpy as np
from scipy import sparse

n = 10_000
main = 2.0 * np.ones(n)
off = -1.0 * np.ones(n - 1)

A = sparse.diags([off, main, off], offsets=[-1, 0, 1], format="csc")
print(A.shape, A.nnz)          # (10000, 10000) 29998

y = sparse.linalg.spsolve(A, np.ones(n))
print(y[:3])                   # a dense vector, because the solution is dense

Full lesson: SciPy arrays and NumPy interop →

Optimisation and curve fitting

Minimising an objective

import numpy as np
from scipy.optimize import minimize, Bounds, LinearConstraint

def rosen(x):
    return np.sum(100.0 * (x[1:] - x[:-1]**2.0)**2.0 + (1.0 - x[:-1])**2.0)

res = minimize(rosen, x0=np.array([-1.2, 1.0]), method="BFGS")

print(res.x)              # close to [1.0, 1.0]
print(res.fun)            # near zero
print(res.success, res.nit, res.message)

Minimising an objective

# constrained problems use the modern bounds/constraints objects, not lambdas
res = minimize(rosen, x0=[0.5, 0.5], method="SLSQP",
               bounds=Bounds([0.0, 0.0], [1.5, 1.5]),
               constraints=LinearConstraint([[1.0, 1.0]], 1.0, np.inf))

print(res.x, res.message)

Full lesson: Optimisation and curve fitting →

Installing SciPy and the subpackage map

Installation and version pairing

python -m venv .venv && source .venv/bin/activate

# the usual route: wheels that already bundle the compiled libraries
pip install numpy scipy

# from a fully specified environment
pip install -r requirements.txt
python -c "import scipy, numpy; print(scipy.__version__, numpy.__version__)"

Full lesson: Installing SciPy and the subpackage map →

Linear algebra with scipy.linalg

Factorisations

lu, piv = linalg.lu_factor(A)
x1 = linalg.lu_solve((lu, piv), b)        # cheap for each new b

c, low = linalg.cho_factor(A)             # A must be SPD
x2 = linalg.cho_solve((c, low), b)

U, s, Vt = linalg.svd(M, full_matrices=False)
print("condition number:", s[0] / s[-1])

w, V = linalg.eigh(A)                     # symmetric: real, sorted eigenvalues
print(w)

Full lesson: Linear algebra with scipy.linalg →

Sparse matrices in depth

Solving sparse systems

# direct solver: exact, but can fill in badly
x = spsolve(A_csc, np.ones(3))

# reuse a factorisation across right-hand sides
lu = sparse.linalg.splu(A_csc.tocsc())
x1 = lu.solve(np.ones(3))
x2 = lu.solve(np.arange(3.0))

# iterative solver: little memory, needs a well-behaved matrix
sym = (A_csr @ A_csr.T).tocsr()
x_it, info = cg(sym, np.ones(3), rtol=1e-10, atol=0.0)
print("converged" if info == 0 else "failed")

Sparsity-preserving practice

# WRONG: dense intermediate, memory explodes
bad = np.linalg.inv(A_csr.toarray())

# RIGHT: stay sparse end to end
good = sparse.linalg.spsolve(A_csc, b)

# check before you commit: how much would dense cost?
dense_gb = A_csr.shape[0] * A_csr.shape[1] * 8 / 1e9
print(f"dense would need {dense_gb:.1f} GB, stored nnz={A_csr.nnz}")

Full lesson: Sparse matrices in depth →

Solving differential equations with solve_ivp

Stiff versus non-stiff

def robertson(t, y):
    k1, k2, k3 = 0.04, 3e7, 1e4
    return [-k1 * y[0] + k3 * y[1] * y[2],
             k1 * y[0] - k2 * y[1] ** 2 - k3 * y[1] * y[2],
             k2 * y[1] ** 2]

fast = solve_ivp(robertson, (0, 1e5), [1, 0, 0], method="RK45")
slow = solve_ivp(robertson, (0, 1e5), [1, 0, 0], method="BDF")
print(fast.nfev, slow.nfev)     # the explicit solver takes orders of magnitude more function calls

Full lesson: Solving differential equations with solve_ivp →

Root finding and special functions

Using them honestly

# always check convergence, do not just take the number
res = optimize.root_scalar(f, bracket=[2, 3], method="brentq")
print(res.converged, res.iterations, res.root, f(res.root))

# scan for brackets when you have no idea where the roots are
xs = np.linspace(-5, 5, 200)
fs = f(xs)
candidates = [(a, b) for a, b, fa, fb in zip(xs[:-1], xs[1:], fs[:-1], fs[1:]) if fa * fb < 0]
print(len(candidates))

# polynomials: get every root at once
print(np.roots([1, 0, -2, -5]))

Full lesson: Root finding and special functions →

Image processing with scipy.ndimage

Filtering and morphology

import numpy as np
from scipy import ndimage as ndi

img = np.random.default_rng(0).random((200, 200))
blurred = ndi.gaussian_filter(img, sigma=2.0)
edges = ndi.sobel(blurred, axis=0)

mask = blurred > 0.6
opened = ndi.binary_opening(mask, structure=np.ones((3, 3)), iterations=1)
filled = ndi.binary_fill_holes(opened)
dilated = ndi.binary_dilation(filled, iterations=2)

Transforms and interpolation

rotated = ndi.rotate(img, angle=30, reshape=True, order=1)
zoomed = ndi.zoom(img, zoom=2.0, order=3)
shifted = ndi.shift(img, shift=(5, -3), mode="nearest")

matrix = np.array([[1, 0.2, 0], [0, 1, 0], [0, 0, 1]])
sheared = ndi.affine_transform(img, matrix, offset=0.0, order=1)

import matplotlib.pyplot as plt
plt.imshow(sheared, cmap="gray", origin="upper")
plt.axis("off")

Full lesson: Image processing with scipy.ndimage →

Spatial algorithms: KDTree, distances and Delaunay

Distance matrices

a = rng.random((1000, 4))
b = rng.random((800, 4))

compact = distance.pdist(a, metric="euclidean")      # 1D, n*(n-1)/2 entries
square = distance.squareform(compact)                # (1000, 1000) dense
cross = distance.cdist(a, b, metric="cosine")        # (1000, 800)

print(compact.shape, square.shape, cross.shape)

Full lesson: Spatial algorithms: KDTree, distances and Delaunay →

FFT and spectral analysis with scipy.fft

Leakage and windowing

window = fft.get_window("hann", sig.size)
windowed = sig * window

f, t_spec, Sxx = fft.spectrogram(sig, fs, nperseg=256, noverlap=128, window="hann")
print(f.shape, t_spec.shape, Sxx.shape)

n = sig.size
nf = np.arange(0, n // 2 + 1)
freq = nf * fs / n
plain = np.abs(np.fft.rfft(sig)) / n
tapered = np.abs(np.fft.rfft(windowed)) / (window.sum())
print(plain[:5].round(4), tapered[:5].round(4))

Full lesson: FFT and spectral analysis with scipy.fft →

Choosing SciPy vs specialised libraries

The awkward boundary cases

# a curve fit with bounds and uncertainty: SciPy's natural job
from scipy.optimize import curve_fit

def decay(t, a, tau):
    return a * np.exp(-t / tau)

popt, pcov = curve_fit(decay, t, y, p0=[1.0, 1.0], bounds=([0, 0], [np.inf, np.inf]))
perr = np.sqrt(np.diag(pcov))
print(popt, perr)

Full lesson: Choosing SciPy vs specialised libraries →

FAQ

Related cheat sheets

Python 3 NumPy pandas Matplotlib Jupyter Notebook Flask

Last refreshed 2026-09-27.