mrkeyoor.com_
Wed 05 Aug 05:06 UTC
PyPIDataupdated 05 Aug 2026

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.

Verdict

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.

API stability5/5Twenty-plus years of the same core API on a 1.x line; changes go through long deprecation cycles and the release notes call out every removal.
Docs4/5Reference docs are complete with algorithm citations and examples per function, and the tutorial section is decent; discoverability across the huge submodule surface is the weak spot.
Maintenance5/5Pushed daily, two feature releases per year on a predictable cadence, NumFOCUS fiscal sponsorship and a large contributor base.
Ecosystem5/5Foundation of the scientific Python stack: scikit-learn, statsmodels and much of the ecosystem depend on it, and any numerics question about it has a decade of answers.

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
Skip it if

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 errors

For 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 estimate

quad 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 trajectory

solve_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 points

Build the tree once and reuse it; rebuilding per query throws away the entire benefit.

Alternatives

PackageRegistryPick it when
numpyPyPIPlain array operations and linear algebra basics are all you need
statsmodelsPyPIYou need regression models, econometrics and statistical reporting rather than raw tests
scikit-learnPyPIThe task is machine learning: classifiers, clustering, pipelines and model evaluation
sympyPyPIYou need symbolic math (exact derivatives, algebra) instead of numerical answers