SciPy: Advanced Scientific Computing in Python

SciPy is the scientific computing library that extends NumPy with proven numerical algorithms: optimization, curve fitting, statistics, differential equation integration, interpolation, signal processing, and sparse matrices. Where NumPy provides the data structure (the ndarray) and elementary operations, SciPy provides high-level numerical methods, most of them backed by reference Fortran and C libraries (LAPACK, MINPACK, FFTPACK). This guide walks through the most commonly used submodules with complete, runnable examples.

SciPy vs. NumPy

The boundary is easy to remember:

• NumPy: the ndarray array, vectorized arithmetic, basic linear algebra, random number generation.
• SciPy: the algorithms that operate on those arrays, with robust implementations drawn from decades of scientific computing.

SciPy depends on NumPy and is imported submodule by submodule:

import numpy as np
from scipy import optimize, stats, integrate, interpolate, signal, sparse

# Check the installed version
import scipy
print(scipy.__version__)     # 1.14.1 (for example)

Some overlap exists (scipy.linalg extends numpy.linalg with additional decompositions; scipy.fft is a better replacement for numpy.fft), and in those cases the SciPy version is generally more complete.

scipy.optimize: Minimization and Curve Fitting

Minimizing a Function with minimize

minimize finds a local minimum of a scalar function of several variables.

from scipy.optimize import minimize

# Rosenbrock function, a classic benchmark
def rosenbrock(p):
    x, y = p
    return (1 - x) ** 2 + 100 * (y - x ** 2) ** 2

resultat = minimize(rosenbrock, x0=[0.0, 0.0], method="Nelder-Mead")

print(resultat.success)      # True
print(resultat.x)            # [0.99998736 0.99997369] -> minimum at (1, 1)
print(resultat.fun)          # 5.9e-10 (value at the minimum, nearly zero)
print(resultat.nit)          # number of iterations

The choice of method depends on the problem: "Nelder-Mead" does not require a gradient, "BFGS" (the default without constraints) estimates it by finite differences, "L-BFGS-B" accepts bounds, and "SLSQP" handles general constraints.

# Bounded minimization: x in [0, 2], y in [0, 2]
res = minimize(rosenbrock, x0=[0.5, 0.5], method="L-BFGS-B",
               bounds=[(0, 2), (0, 2)])
print(res.x)                 # [1. 1.]

To find the root of an equation rather than a minimum, use optimize.root_scalar or optimize.brentq.

Fitting a Model to Data with curve_fit

curve_fit estimates a model’s parameters using nonlinear least squares. Example: fitting an exponential decay to noisy measurements.

from scipy.optimize import curve_fit

# Model: the first variable is always x, the following ones are the parameters
def modele(x, a, tau, c):
    return a * np.exp(-x / tau) + c

# Noisy synthetic data
rng = np.random.default_rng(42)
x_donnees = np.linspace(0, 10, 50)
y_vrai = modele(x_donnees, a=3.0, tau=2.5, c=0.5)
y_donnees = y_vrai + rng.normal(0, 0.15, size=x_donnees.size)

# Fit (p0: initial parameter guess)
params, covariance = curve_fit(modele, x_donnees, y_donnees, p0=[1, 1, 0])

a_est, tau_est, c_est = params
incertitudes = np.sqrt(np.diag(covariance))

print(f"a   = {a_est:.3f} ± {incertitudes[0]:.3f}")   # a   = 3.011 ± 0.052
print(f"tau = {tau_est:.3f} ± {incertitudes[1]:.3f}") # tau = 2.487 ± 0.098
print(f"c   = {c_est:.3f} ± {incertitudes[2]:.3f}")   # c   = 0.503 ± 0.031

The diagonal of the covariance matrix gives the variance of each estimated parameter: its square root provides an uncertainty of one standard deviation. Always provide a reasonable p0 for nonlinear models, or the fit may converge to a poor optimum.

scipy.stats: Distributions and Statistical Tests

Working with Distributions

Every distribution exposes the same interface: pdf (density), cdf (cumulative distribution function), ppf (quantiles, the inverse of the cdf), and rvs (sampling).

from scipy import stats

# Normal distribution with mean 100 and standard deviation 15
loi = stats.norm(loc=100, scale=15)

print(loi.pdf(100))          # 0.02659... (density at the peak)
print(loi.cdf(115))          # 0.8413... (P(X <= 115))
print(loi.ppf(0.975))        # 129.399... (97.5% quantile)
print(loi.rvs(size=3, random_state=42))   # [107.45  84.44 109.72]

# Other distributions of the same mold: stats.t, stats.chi2, stats.binom,
# stats.poisson, stats.expon, stats.uniform...

t-Test: Comparing Two Means

rng = np.random.default_rng(1)
groupe_a = rng.normal(loc=52, scale=8, size=40)   # e.g., scores with method A
groupe_b = rng.normal(loc=48, scale=8, size=40)   # e.g., scores with method B

t_stat, p_valeur = stats.ttest_ind(groupe_a, groupe_b)
print(f"t = {t_stat:.3f}, p = {p_valeur:.4f}")
# t = 2.279, p = 0.0255 (for example)

if p_valeur < 0.05:
    print("Différence significative au seuil de 5 %")
# Paired variant (before/after measurements on the same subjects): stats.ttest_rel
# Unequal variances: ttest_ind(..., equal_var=False) (Welch's test)

Chi-Square Test: Independence of Two Categorical Variables

# Contingency table: rows = training received, columns = outcome
observations = np.array([[45, 15],    # with training: 45 passes, 15 failures
                         [30, 30]])   # without training: 30 passes, 30 failures

chi2, p, ddl, attendus = stats.chi2_contingency(observations)
print(f"chi2 = {chi2:.3f}, p = {p:.4f}, ddl = {ddl}")
# chi2 = 6.806, p = 0.0091, ddl = 1
print(attendus)              # expected counts under the independence hypothesis
# [[37.5 22.5]
#  [37.5 22.5]]

Correlations

x = rng.normal(size=100)
y = 0.7 * x + rng.normal(scale=0.5, size=100)

r, p = stats.pearsonr(x, y)          # linear correlation
print(f"Pearson  r = {r:.3f}, p = {p:.2e}")   # r = 0.812, p = 1.5e-24

rho, p = stats.spearmanr(x, y)       # rank correlation (monotonic relationships)
print(f"Spearman rho = {rho:.3f}")   # rho = 0.795

Use Spearman when the relationship is monotonic but nonlinear, or when outliers are present.

scipy.integrate: Integrals and Differential Equations

Numerical Integration with quad

from scipy.integrate import quad

# Integral of exp(-x²) from 0 to infinity (equals sqrt(pi)/2)
valeur, erreur = quad(lambda x: np.exp(-x ** 2), 0, np.inf)
print(valeur)                        # 0.8862269254527579
print(np.sqrt(np.pi) / 2)            # 0.8862269254527579
print(erreur)                        # 7.1e-09 (error estimate)

Solving an ODE with solve_ivp

Example: the damped harmonic oscillator, described by the second-order equation x” + 2ζω x’ + ω² x = 0, rewritten as a first-order system.

from scipy.integrate import solve_ivp

omega = 2.0    # natural angular frequency
zeta = 0.15    # damping ratio

def oscillateur(t, y):
    x, v = y                              # position, velocity
    return [v, -2 * zeta * omega * v - omega ** 2 * x]

solution = solve_ivp(
    oscillateur,
    t_span=(0, 20),                       # integration interval
    y0=[1.0, 0.0],                        # position 1, zero velocity at t=0
    t_eval=np.linspace(0, 20, 400),       # points at which to evaluate the solution
    method="RK45",                        # adaptive Runge-Kutta (default)
)

print(solution.success)                   # True
print(solution.y.shape)                   # (2, 400) -> position and velocity
print(solution.y[0, :3])                  # [1.  0.995 0.980] (damped position)

# Plot position over time
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(solution.t, solution.y[0])
ax.set_xlabel("Temps (s)")
ax.set_ylabel("Position")
ax.set_title("Oscillateur harmonique amorti")
plt.show()

For stiff systems (very different time constants), pass method="Radau" or method="BDF".

scipy.interpolate: Data Interpolation

Interpolation builds a continuous function that passes through discrete measurement points.

from scipy.interpolate import CubicSpline, interp1d

# Sparse measurement points
x_mesures = np.array([0, 1, 2, 3, 4, 5], dtype=float)
y_mesures = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])

# Cubic spline: smooth, differentiable, recommended
spline = CubicSpline(x_mesures, y_mesures)

x_fin = np.linspace(0, 5, 100)
y_interpole = spline(x_fin)

print(spline(2.5))                   # 0.5546... (interpolated value between 2 and 3)
print(spline(2.5, 1))                # first derivative at the same point

# Simple linear interpolation
lineaire = interp1d(x_mesures, y_mesures, kind="linear")
print(lineaire(2.5))                 # 0.5 (average of the neighboring points)

Warning: interpolation is only reliable within the domain of the measurements. To extrapolate or to smooth noisy data, prefer fitting a model (curve_fit) or a smoothing spline (scipy.interpolate.make_smoothing_spline).

scipy.signal and scipy.fft: Signal Processing

Spectral Analysis with the FFT

from scipy.fft import fft, fftfreq

# Sampled signal: sum of two sinusoids (5 Hz and 50 Hz) plus noise
fe = 500                              # sampling frequency (Hz)
duree = 2.0
t = np.linspace(0, duree, int(fe * duree), endpoint=False)
rng = np.random.default_rng(7)
signal_brut = (np.sin(2 * np.pi * 5 * t)
               + 0.5 * np.sin(2 * np.pi * 50 * t)
               + 0.3 * rng.normal(size=t.size))

# Fourier transform
spectre = fft(signal_brut)
frequences = fftfreq(t.size, d=1 / fe)

# Amplitudes at positive frequencies
moitie = t.size // 2
amplitudes = 2 / t.size * np.abs(spectre[:moitie])

pics = frequences[:moitie][amplitudes > 0.3]
print(pics)                           # [ 5. 50.] -> both frequencies recovered

Filtering with a Butterworth Filter

from scipy.signal import butter, filtfilt

# 4th-order low-pass filter, cutoff frequency 20 Hz
b, a = butter(N=4, Wn=20, btype="low", fs=fe)

# filtfilt applies the filter in both directions: no phase shift
signal_filtre = filtfilt(b, a, signal_brut)

# The 50 Hz component and the high-frequency noise are removed,
# the 5 Hz sinusoid is preserved.

scipy.signal also provides peak detection (find_peaks), spectrograms (ShortTimeFFT), convolution, and the design of many filter types (Chebyshev, elliptic, FIR).

from scipy.signal import find_peaks

pics_indices, proprietes = find_peaks(signal_filtre, height=0.5, distance=fe // 10)
print(len(pics_indices))              # 10 -> one peak per period of the 5 Hz sinusoid

scipy.sparse: Sparse Matrices

When a matrix contains an overwhelming majority of zeros (graphs, systems of equations derived from meshes, term-document matrices), storing it densely wastes memory and computation.

from scipy import sparse

# Sparse matrix in COO format (construction) then CSR (computation)
lignes = np.array([0, 1, 2, 2])
colonnes = np.array([1, 2, 0, 2])
valeurs = np.array([3.0, 4.0, 5.0, 6.0])

m = sparse.coo_array((valeurs, (lignes, colonnes)), shape=(3, 3)).tocsr()
print(m.toarray())
# [[0. 3. 0.]
#  [0. 0. 4.]
#  [5. 0. 6.]]

print(m.nnz)                          # 4 stored nonzero elements (out of 9)

# Matrix-vector product: only nonzero elements are traversed
v = np.array([1.0, 2.0, 3.0])
print(m @ v)                          # [ 6. 12. 23.]

# Solving a sparse linear system
from scipy.sparse.linalg import spsolve
b = np.array([6.0, 12.0, 23.0])
x = spsolve(m.tocsc(), b)
print(x)                              # [1. 2. 3.]

Orders of magnitude: a dense 100,000 × 100,000 matrix in float64 would take up 80 GB; the same matrix with 0.01% nonzero elements fits in a few megabytes in CSR format.

Submodule Summary Table

Submodule	Domain	Key functions
scipy.optimize	Optimization, fitting, roots	minimize, curve_fit, root_scalar, linprog
scipy.stats	Distributions, tests, correlations	norm, ttest_ind, chi2_contingency, pearsonr
scipy.integrate	Integrals, differential equations	quad, solve_ivp
scipy.interpolate	Interpolation, splines	CubicSpline, interp1d, make_smoothing_spline
scipy.signal	Filtering, peaks, spectrograms	butter, filtfilt, find_peaks
scipy.fft	Fourier transforms	fft, ifft, fftfreq, rfft
scipy.sparse	Sparse matrices	csr_array, coo_array, sparse.linalg.spsolve
scipy.linalg	Extended linear algebra	lu, svd, expm, solve_banded
scipy.spatial	Geometry, distances	KDTree, distance.cdist, ConvexHull
scipy.ndimage	n-D image processing	gaussian_filter, label, rotate

Best Practices

• Import submodules explicitly (from scipy import optimize): import scipy alone does not load the submodules.
• Always check resultat.success after minimize or solve_ivp: an algorithm that has not converged still returns a result.
• Provide a sensible initial guess p0 to curve_fit: nonlinear least squares converges to a local optimum that depends on the starting point.
• Interpret p-values with caution: a p-value below 0.05 measures neither the size of the effect nor its practical significance; also report means and confidence intervals.
• Choose the method suited to the problem: Radau or BDF for stiff ODEs, L-BFGS-B whenever there are bounds, Welch’s test when variances differ.
• Don’t extrapolate with interpolation: outside the domain of the measurements, a cubic spline diverges quickly.
• Switch to sparse matrices as soon as the density of nonzero elements drops below a few percent on large matrices.

Conclusion

SciPy turns Python into a complete scientific computing environment: optimization and model fitting with optimize, statistical inference with stats, differential equation solving with integrate, function reconstruction with interpolate, frequency analysis and filtering with signal and fft, and scaling up with sparse. Each submodule encapsulates reference algorithms validated over decades, saving you from reimplementing delicate numerical methods. Combined with NumPy for data structures, Pandas for labeled tables, and Matplotlib for visualization, SciPy completes the standard toolbox for any scientific or engineering work in Python.