SciPy Tutorial — Scientific Computing in Python

SciPy (Scientific Python) is an open-source library that builds on NumPy and adds a large collection of algorithms for mathematics, science, and engineering. Where NumPy gives you the ndarray and basic operations, SciPy gives you the specialist solvers: numerical integration, optimization, interpolation, signal processing, statistics, spatial algorithms, and more.

This chapter covers:

Installing SciPy and the import convention
Linear algebra (scipy.linalg)
Numerical integration and differentiation (scipy.integrate)
Optimization — finding minima and roots (scipy.optimize)
Interpolation (scipy.interpolate)
Statistics (scipy.stats)
N-dimensional image processing (scipy.ndimage)

Installing SciPy

SciPy is included in the Anaconda distribution. To install it with pip:

pip install scipy

SciPy depends on NumPy, which pip installs automatically if it is not already present.

Import convention

SciPy is organized into sub-packages. Import only the sub-packages you need rather than the entire library:

import numpy as np
from scipy import linalg, integrate, optimize, interpolate, stats, ndimage

You can also check the installed version:

import scipy
print(scipy.__version__)  # e.g. 1.13.0

Linear Algebra

scipy.linalg extends NumPy's linear algebra routines with additional decompositions and solvers. All functions work on regular NumPy arrays.

Determinant and inverse

import numpy as np
from scipy import linalg

a = np.array([[1, 2],
              [3, 4]])

det = linalg.det(a)
print(det)          # -2.0

inv = linalg.inv(a)
print(inv)
# [[-2.   1. ]
#  [ 1.5 -0.5]]

det([[1,2],[3,4]]) = 1×4 − 2×3 = −2. The inverse satisfies a @ inv == I.

Eigenvalues and eigenvectors

Eigenvalues describe how a matrix stretches space; eigenvectors give the directions that do not rotate.

import numpy as np
from scipy import linalg

a = np.array([[1, 2],
              [3, 4]])

eigenvalues, eigenvectors = linalg.eig(a)
print(eigenvalues)
# [-0.37228132+0.j  5.37228132+0.j]
print(eigenvectors)
# [[-0.82456484 -0.41597356]
#  [ 0.56576746 -0.90937671]]

Each column of eigenvectors corresponds to the matching eigenvalue. The eigenvalues are returned as complex numbers even when the imaginary part is zero.

Singular Value Decomposition (SVD)

SVD factorizes a matrix A into three matrices U, s, Vt such that A = U @ diag(s) @ Vt. It is the foundation of principal component analysis (PCA) and many dimensionality-reduction techniques.

import numpy as np
from scipy import linalg

a = np.array([[1, 2],
              [3, 4]])

u, s, vt = linalg.svd(a)
print(u)
# [[-0.40455358 -0.9145143 ]
#  [-0.9145143   0.40455358]]
print(s)    # [5.4649857  0.36596619]
print(vt)
# [[-0.57604844 -0.81741556]
#  [ 0.81741556 -0.57604844]]

Solving a linear system

linalg.solve is the correct way to solve Ax = b. It is faster and more numerically stable than computing the inverse and multiplying.

import numpy as np
from scipy import linalg

# Solve: 1x + 2y = 5
#        3x + 4y = 11
A = np.array([[1, 2],
              [3, 4]])
b = np.array([5, 11])

x = linalg.solve(A, b)
print(x)   # [1. 2.]

# Verify: A @ x should equal b
print(np.allclose(A @ x, b))  # True

Numerical Integration

scipy.integrate provides routines for computing definite integrals when an analytical solution is impractical.

Single-variable integration with `quad`

integrate.quad uses adaptive quadrature to integrate a function over an interval. It returns the result and an estimate of the absolute error.

import numpy as np
from scipy import integrate

# Integrate f(x) = x^2 + 2x + 1 from 0 to 1
# Analytical result: [x^3/3 + x^2 + x] from 0 to 1 = 1/3 + 1 + 1 = 7/3
def f(x):
    return x**2 + 2*x + 1

result, error = integrate.quad(f, 0, 1)
print(result)   # 2.3333333333333335
print(error)    # ~2.6e-14  (absolute error estimate)

Numerical differentiation

SciPy's approx_fprime computes a finite-difference gradient. For scalar functions the central-difference derivative from scipy.misc is simpler:

import numpy as np
from scipy.optimize import approx_fprime

# Derivative of sin(x) at x = 0 should be cos(0) = 1
result = approx_fprime([0.0], lambda x: np.sin(x[0]), 1e-8)
print(result[0])   # ~1.0

Optimization

scipy.optimize finds minima, maxima (via minimizing the negative), and roots of functions.

Minimizing a multivariate function

optimize.minimize supports many methods (Nelder-Mead, BFGS, L-BFGS-B, …). The default method is chosen automatically.

from scipy import optimize

# Minimize f(x) = x^2 + 2x + 1 = (x + 1)^2  — minimum at x = -1
def f(x):
    return x**2 + 2*x + 1

result = optimize.minimize(f, x0=0)   # x0 is the starting guess
print(result.success)   # True
print(result.x)         # [-1.00000001]  (near -1)
print(result.fun)       # ~0.0           (minimum value)

Minimizing a scalar function over a bounded interval

optimize.minimize_scalar is simpler for single-variable problems. Always supply bounds with method='bounded' when the function has no global minimum (i.e., is unbounded):

from scipy import optimize

# Minimize h(x) = x^2 - 4x + 3 over [0, 4]  — minimum at x = 2, h(2) = -1
def h(x):
    return x**2 - 4*x + 3

result = optimize.minimize_scalar(h, bounds=(0, 4), method='bounded')
print(result.x)     # ~2.0
print(result.fun)   # -1.0

Finding roots

optimize.root_scalar finds where a function crosses zero:

from scipy.optimize import root_scalar

# Solve x^2 - 4 = 0 in the interval [0, 3]  — root at x = 2
res = root_scalar(lambda x: x**2 - 4, bracket=[0, 3])
print(res.root)     # 2.0

Interpolation

scipy.interpolate fits a smooth curve through data points so you can estimate values between them.

1-D interpolation

interp1d creates a callable interpolation function from discrete (x, y) pairs. The kind parameter selects the method: 'linear' (default), 'quadratic', or 'cubic'.

import numpy as np
from scipy.interpolate import interp1d

# Sample points from y = x^2
x = np.array([0, 1, 2, 3, 4])
y = np.array([0, 1, 4, 9, 16])

f_linear = interp1d(x, y)           # piecewise linear
f_cubic  = interp1d(x, y, kind='cubic')  # cubic spline

# Estimate y at x = 2.5  (exact value: 2.5^2 = 6.25)
print(float(f_linear(2.5)))   # 6.5   (linear — slightly off)
print(float(f_cubic(2.5)))    # 6.25  (cubic — matches exactly for polynomials)

Cubic interpolation recovers the exact result because the data comes from a quadratic polynomial, and the cubic spline is flexible enough to fit it perfectly.

Statistics

scipy.stats contains over 80 continuous and discrete probability distributions plus a collection of statistical tests.

Descriptive statistics and normal distribution

import numpy as np
from scipy.stats import norm

# Standard normal distribution (mean=0, std=1)
print(norm.pdf(0))     # 0.3989422804014327  — probability density at x = 0
print(norm.cdf(1.96))  # 0.9750021048517795  — P(X <= 1.96)
print(norm.ppf(0.975)) # 1.9599639845400536  — inverse CDF (quantile function)

# Fit a normal distribution to data
data = np.array([2.1, 3.3, 2.8, 3.1, 2.5, 3.0, 2.7])
mu, sigma = norm.fit(data)
print(f'Fitted mean: {mu:.4f}, std: {sigma:.4f}')

Hypothesis testing

scipy.stats includes t-tests, chi-squared tests, ANOVA, Kolmogorov-Smirnov tests, and more.

import numpy as np
from scipy.stats import ttest_1samp, ttest_ind

# One-sample t-test: is the sample mean significantly different from 5?
sample = np.array([4.8, 5.1, 4.9, 5.3, 5.2, 4.7, 5.0])
t_stat, p_value = ttest_1samp(sample, popmean=5.0)
print(f't = {t_stat:.4f}, p = {p_value:.4f}')
# p > 0.05 → no significant difference from 5

# Two-sample t-test: are these two groups different?
group_a = np.array([5.1, 5.3, 4.9, 5.2, 5.0])
group_b = np.array([6.1, 6.3, 5.8, 6.0, 5.9])
t_stat2, p_value2 = ttest_ind(group_a, group_b)
print(f't = {t_stat2:.4f}, p = {p_value2:.6f}')
# Very small p → groups are significantly different

Random variables and sampling

Every distribution in scipy.stats exposes the same interface: pdf, cdf, ppf, rvs (random variates), and fit.

from scipy.stats import norm, poisson

# Draw 5 samples from a normal distribution with mean=10, std=2
samples = norm.rvs(loc=10, scale=2, size=5, random_state=42)
print(samples.round(2))   # [10.99  9.72 11.3  13.05  9.53]

# Poisson distribution: P(X = k) for mean lambda=3
for k in range(6):
    print(f'P(X={k}) = {poisson.pmf(k, mu=3):.4f}')

N-Dimensional Image Processing

scipy.ndimage operates on arrays of any dimensionality (images, volumes, time-series cubes). Below is a self-contained example that uses a synthetic 2-D array so no external image file is needed.

Gaussian blur and connected-region labeling

import numpy as np
from scipy import ndimage

# Create a synthetic 5x5 "image" with a bright spot in the centre
image = np.array([
    [0, 0, 0, 0, 0],
    [0, 1, 1, 1, 0],
    [0, 1, 5, 1, 0],
    [0, 1, 1, 1, 0],
    [0, 0, 0, 0, 0],
], dtype=float)

# Smooth the image with a Gaussian filter (sigma controls the blur radius)
blurred = ndimage.gaussian_filter(image, sigma=1)
print('Center value after blur:', round(blurred[2, 2], 4))
# 1.4166  — the bright peak is spread across neighbouring pixels

# Label connected non-zero regions (like counting distinct objects)
binary = image > 0
labeled, num_regions = ndimage.label(binary)
print('Number of connected regions:', num_regions)  # 1
print(labeled)
# [[0 0 0 0 0]
#  [0 1 1 1 0]
#  [0 1 1 1 0]
#  [0 1 1 1 0]
#  [0 0 0 0 0]]

Common `ndimage` operations

Function	What it does
`gaussian_filter(a, sigma)`	Smooth with a Gaussian kernel
`sobel(a)`	Detect edges (Sobel gradient)
`label(a)`	Label connected regions in a binary array
`binary_dilation(a)`	Expand foreground regions
`zoom(a, factor)`	Resize an array
`rotate(a, angle)`	Rotate an array (in degrees)

When to Use SciPy vs. Other Libraries

Task	Preferred tool
Array creation and basic math	NumPy
DataFrames, time-series, I/O	Pandas
Scientific algorithms (integration, optimization)	SciPy
Machine learning	scikit-learn (builds on SciPy)
Plotting	Matplotlib

SciPy is not a replacement for NumPy — it depends on NumPy and extends it. In practice you will import both.

Quick Reference

Sub-package	Key functions
`scipy.linalg`	`det`, `inv`, `eig`, `svd`, `solve`
`scipy.integrate`	`quad`, `dblquad`, `solve_ivp`
`scipy.optimize`	`minimize`, `minimize_scalar`, `root_scalar`
`scipy.interpolate`	`interp1d`, `CubicSpline`, `griddata`
`scipy.stats`	`norm`, `ttest_1samp`, `ttest_ind`, `chi2_contingency`
`scipy.ndimage`	`gaussian_filter`, `label`, `sobel`, `zoom`
`scipy.signal`	`butter`, `lfilter`, `find_peaks`
`scipy.spatial`	`distance`, `KDTree`, `ConvexHull`