scipy
SciPy is the scientific computing layer that sits on top of NumPy arrays: submodules for optimization, numerical integration, ODE solvers, linear algebra, statistics, signal and image processing, sparse matrices, interpolation and Fourier transforms. Most of the heavy code is compiled Fortran/C/C++ wrapped in Python, so you get battle-tested numerical routines (many from classic libraries like LAPACK) at native speed with a Python API.
The bedrock of scientific Python: NumFOCUS-backed, actively developed since 2001, and the correct default whenever you need a numerical algorithm you did not invent. Just do not install it for array math NumPy already does, and go GPU-native frameworks for deep learning scale.
Use it if
- You need serious numerics: curve fitting, minimization, ODE solving, hypothesis tests or filtering, and want algorithms with decades of vetting instead of hand-rolled math
- You already live in the NumPy world; SciPy functions consume and return plain NumPy arrays
- You need sparse matrices or spatial structures (KD-trees, distance computations) that pure NumPy does not provide
- You want scientific rigor with citations; docstrings reference the actual papers and algorithms
- You only need array math, means and dot products; NumPy alone covers it and SciPy adds a heavyweight compiled dependency to your image for nothing
- You need GPU acceleration or automatic differentiation; SciPy is CPU-only, look at JAX, PyTorch or CuPy instead
- You want machine learning; scikit-learn (which builds on SciPy) has the estimators, cross-validation and metrics you are actually after
- You are in a size-constrained deployment like a slim Lambda; the compiled binaries plus NumPy eat a large share of the budget
Setup reality
pip install scipy is painless on major platforms because binary wheels exist for Linux, macOS and Windows. Pain starts off the happy path: exotic platforms or building from source require Fortran and C compilers plus BLAS/LAPACK, which is a genuinely miserable build. The API is stable but sprawling across 20-ish submodules, and knowing that scipy.optimize.minimize hides a dozen method choices is half the learning curve. Watch deprecation warnings between minor releases; legacy interfaces like interp1d are discouraged in favor of newer ones.
Patterns
Minimize a scalar functionminimize-function
import numpy as np
from scipy.optimize import minimize
def rosen(x):
return sum(100.0 * (x[1:] - x[:-1]**2)**2 + (1 - x[:-1])**2)
res = minimize(rosen, x0=np.array([1.3, 0.7, 0.8]), method='Nelder-Mead')
print(res.x, res.success)Always check res.success; optimizers fail quietly and still return a result object.
Fit a model to noisy datacurve-fit
import numpy as np
from scipy.optimize import curve_fit
def model(x, a, b):
return a * np.exp(-b * x)
xdata = np.linspace(0, 4, 50)
ydata = model(xdata, 2.5, 1.3) + 0.1 * np.random.default_rng(0).normal(size=50)
popt, pcov = curve_fit(model, xdata, ydata, p0=[1.0, 1.0])
perr = np.sqrt(np.diag(pcov)) # 1-sigma parameter errorsFor exponential-ish models a bad initial guess p0 sends the fit to nonsense; never rely on the default all-ones start.
Numerically integrate a functionintegrate-quad
import numpy as np
from scipy.integrate import quad
value, abserr = quad(lambda x: np.exp(-x**2), 0, np.inf)
print(value, abserr) # ~0.8862, error estimatequad handles infinite limits directly; for sharp peaks pass their locations via the points argument or it can miss them.
Solve an ODE system with solve_ivpsolve-ode
import numpy as np
from scipy.integrate import solve_ivp
def lotka(t, z):
x, y = z
return [1.5 * x - x * y, -3 * y + x * y]
sol = solve_ivp(lotka, t_span=(0, 15), y0=[10, 5],
dense_output=True, rtol=1e-8)
z = sol.sol(np.linspace(0, 15, 300)) # smooth trajectorysolve_ivp is the modern interface (odeint is legacy); switch to method='Radau' or 'BDF' when the system is stiff.
Compare two samples with a t-testhypothesis-test
import numpy as np
from scipy import stats
rng = np.random.default_rng(0)
a = rng.normal(0.0, 1.0, 200)
b = rng.normal(0.2, 1.0, 200)
t = stats.ttest_ind(a, b, equal_var=False)
print(t.statistic, t.pvalue)equal_var=False (Welch's test) is the safer default; the classic test assumes equal variances your data rarely has.
Smooth interpolation with CubicSplineinterpolate
import numpy as np
from scipy.interpolate import CubicSpline
x = np.array([0, 1, 2, 3, 4])
y = np.array([0.0, 0.8, 0.9, 0.1, -0.8])
cs = CubicSpline(x, y)
xs = np.linspace(0, 4, 100)
ys = cs(xs)Reach for CubicSpline or make_interp_spline; interp1d is legacy and the docs advise against it for new code.
Low-pass filter a noisy signalfilter-signal
from scipy.signal import butter, filtfilt
fs = 500.0 # sample rate Hz
cutoff = 20 # Hz
b, a = butter(N=4, Wn=cutoff, fs=fs, btype='low')
clean = filtfilt(b, a, noisy_signal)filtfilt runs the filter forward and backward for zero phase shift; plain lfilter delays your signal.
Get the frequency spectrum of a signalfft-spectrum
import numpy as np
from scipy.fft import rfft, rfftfreq
fs = 1000
t = np.arange(0, 1, 1 / fs)
sig = np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 120 * t)
spectrum = np.abs(rfft(sig))
freqs = rfftfreq(len(sig), 1 / fs)Use scipy.fft (not the deprecated scipy.fftpack), and rfft for real signals to skip redundant negative frequencies.
Build and solve a sparse systemsparse-matrix
import numpy as np
from scipy import sparse
from scipy.sparse.linalg import spsolve
A = sparse.diags([1, -2, 1], offsets=[-1, 0, 1], shape=(1000, 1000), format='csr')
b = np.ones(1000)
x = spsolve(A, b)Use COO or LIL format while building, then convert to CSR/CSC for math; inserting into CSR row by row is very slow.
Solve a dense linear system properlylinear-solve
import numpy as np
from scipy import linalg
A = np.array([[3.0, 2.0], [1.0, 4.0]])
b = np.array([5.0, 6.0])
x = linalg.solve(A, b)solve(A, b) is faster and more numerically sound than computing inv(A) @ b; almost never invert a matrix explicitly.
Fast nearest-neighbor lookup with a KD-treenearest-neighbors
import numpy as np
from scipy.spatial import cKDTree
points = np.random.default_rng(0).random((10000, 2))
tree = cKDTree(points)
dist, idx = tree.query([0.5, 0.5], k=5) # 5 nearest pointsBuild the tree once and reuse it; rebuilding per query throws away the entire benefit.
Alternatives
| Package | Registry | Pick it when |
|---|---|---|
| numpy | PyPI | Plain array operations and linear algebra basics are all you need |
| statsmodels | PyPI | You need regression models, econometrics and statistical reporting rather than raw tests |
| scikit-learn | PyPI | The task is machine learning: classifiers, clustering, pipelines and model evaluation |
| sympy | PyPI | You need symbolic math (exact derivatives, algebra) instead of numerical answers |