PYTHON / SCIPY
What SciPy adds on top of NumPy
Explain how SciPy layers domain algorithms on NumPy's ndarray, import its subpackages correctly, and choose between scipy.linalg and numpy.linalg.
What you will learn
- Import SciPy by subpackage: from scipy import linalg, special
- Know that SciPy takes and returns plain ndarrays, so NumPy ops chain onto its results
- Pick scipy.linalg for LAPACK extras like expm, assume_a, solve_triangular
- Use scipy.special (gammaln, logsumexp) instead of formulas that overflow float64
Understanding What SciPy adds on top of NumPy
NumPy defines the data structure and the elementwise vocabulary: the ndarray, dtypes, shapes, broadcasting, ufuncs, and a small linear algebra core. SciPy introduces no new array type at all; it is a library of numerical algorithms whose inputs are ndarrays and whose outputs are ndarrays. Most of those algorithms are thin, careful wrappers around long-established compiled libraries such as LAPACK, QUADPACK, MINPACK, FITPACK and ARPACK. The mental model is two layers: NumPy is the container plus arithmetic, SciPy is the shelf of solvers that operate on that container.
SciPy is organised by problem domain into subpackages: linalg, special, optimize, integrate, interpolate, stats, sparse, signal, spatial, fft, ndimage, io, constants, cluster. You import the one you need, for example `from scipy import linalg`, rather than reaching for a flat top-level namespace. That layout is deliberate: each subpackage pulls in its own compiled extension modules, and importing all of them for every script would make startup noticeably slower. Nothing useful for day-to-day work is exposed at the top level except version metadata.
Two areas overlap with NumPy on purpose: scipy.linalg versus numpy.linalg, and scipy.fft versus numpy.fft. SciPy's linalg exposes far more of LAPACK — expm, solve_triangular, lu, schur, polar, structure hints such as solve(..., assume_a='pos'), and performance switches like overwrite_a and check_finite=False — because it always links against a real LAPACK, while numpy.linalg is the dependency-light baseline whose routines broadcast over stacks of matrices. Beyond linear algebra, scipy.special supplies numerically stable formulations (gammaln, logsumexp, expit) that you cannot safely reassemble from raw ufuncs, and scipy.constants ships reference physical values instead of algorithms.
import numpy as np
from scipy import linalg, special
A = np.array([[4.0, 1.0],
[1.0, 3.0]])
b = np.array([1.0, 2.0])
# Same problem, two libraries: NumPy assumes nothing, SciPy can be told the structure.
x_numpy = np.linalg.solve(A, b)
x_scipy = linalg.solve(A, b, assume_a="pos") # Cholesky path, no NumPy equivalent
print("agree:", np.allclose(x_numpy, x_scipy))
print("plain ndarray:", type(x_scipy) is np.ndarray)
# Algorithms NumPy simply does not ship.
E = linalg.expm(np.diag([0.0, 1.0]))
print("expm diagonal:", [round(float(v), 6) for v in np.diag(E)])
print("log-gamma(10):", round(float(special.gammaln(10.0)), 6))
print("erf on an array:", [round(float(v), 6) for v in special.erf(np.array([0.0, 1.0]))])SciPy is not an alternative array library; it is a domain-organised collection of algorithms that consume and return NumPy ndarrays, filling in what NumPy deliberately leaves out.
Worked examples
Stable formulations in scipy.special
Shows why SciPy ships gammaln and logsumexp rather than leaving you to compose them from NumPy ufuncs.
import math
import numpy as np
from scipy import special
n = 200.0
print("gamma(200) :", float(special.gamma(n)))
print("log of that :", float(np.log(special.gamma(n))))
print("gammaln(200) :", round(float(special.gammaln(n)), 2))
try:
math.log(math.exp(1000.0) + math.exp(1000.0))
except OverflowError:
print("naive logsumexp : OverflowError")
print("special.logsumexp:", round(float(special.logsumexp([1000.0, 1000.0])), 4))Example explained
Line 1Gamma(200) is about 10**372, past the float64 ceiling, so special.gamma returns inf and the logarithm of it is unrecoverable.
Line 2special.gammaln evaluates log-Gamma directly, so the answer 857.93 never leaves float range.
Line 3math.exp(1000.0) raises OverflowError, which is exactly what a hand-written log(sum(exp(x))) hits.
Line 4special.logsumexp subtracts the maximum before exponentiating, returning 1000 + log(2) = 1000.6931.
LAPACK routines that numpy.linalg does not expose
Uses a structure-aware triangular solver and an explicit LU factorisation, neither of which exists in numpy.linalg.
import numpy as np
from scipy import linalg
L = np.array([[2.0, 0.0, 0.0],
[1.0, 3.0, 0.0],
[4.0, 1.0, 5.0]])
b = np.array([2.0, 7.0, 21.0])
x = linalg.solve_triangular(L, b, lower=True)
print(x)
print("matches general solver:", np.allclose(x, np.linalg.solve(L, b)))
print("numpy has solve_triangular:", hasattr(np.linalg, "solve_triangular"))
P, Lf, U = linalg.lu(L)
print("P @ L @ U rebuilds L:", np.allclose(P @ Lf @ U, L))Example explained
Line 1solve_triangular runs forward substitution only, O(n**2), because you promised the matrix is triangular; np.linalg.solve must factor first, O(n**3).
Line 2lower=True makes LAPACK read only the lower triangle — anything above the diagonal is ignored rather than validated.
Line 3hasattr returns False: structure-aware drivers are part of what SciPy adds, not part of numpy.linalg.
Line 4linalg.lu hands you the factorisation itself, which numpy.linalg computes internally and throws away.
SciPy output is ordinary NumPy data
Feeds a SciPy eigensolver result straight into NumPy broadcasting and ndarray methods with no conversion.
import numpy as np
from scipy import linalg
A = np.array([[2.0, 1.0],
[1.0, 2.0]])
w, V = linalg.eigh(A)
print("eigenvalues:", w)
print("returned by:", type(w).__module__, type(w).__name__)
rebuilt = (V * w) @ V.T
print("rebuilt A:", np.allclose(rebuilt, A))
print("mean eigenvalue:", float(w.mean()))Example explained
Line 1linalg.eigh wraps LAPACK's symmetric eigensolver and returns eigenvalues in ascending order as a float64 ndarray.
Line 2type(w).__module__ is 'numpy', which is why there is never a conversion step between the two libraries.
Line 3V * w broadcasts the eigenvalues across the columns of V: a pure NumPy operation applied to SciPy output.
Line 4w.mean() works because w is an ndarray with the full NumPy method set, not a SciPy-specific object.
Important notes
scipy.fftpack is a legacy interface kept only for backwards compatibility; new code should use scipy.fft, or numpy.fft for a plain transform.
SciPy is compiled against a specific NumPy ABI range, so upgrade the pair together; a mismatch surfaces as an ImportError complaining about binary incompatibility rather than a missing function.
Common mistakes
Writing `import scipy` and then `scipy.optimize.minimize(...)`: on many versions the subpackage is never loaded, and you get AttributeError: module 'scipy' has no attribute 'optimize'.
Passing a stack such as shape (500, 3, 3) to a scipy.linalg routine expecting one result per matrix the way np.linalg.det gives; two-dimensional-only routines raise ValueError instead of looping for you.
Hunting for array basics in SciPy (scipy.array, scipy.mean); those NumPy aliases were removed, so you get AttributeError — creation and reductions stay in NumPy.
Try it yourself
Change, predict, then run
Build A = [[4.0, 1.0, 0.0], [1.0, 4.0, 1.0], [0.0, 1.0, 4.0]] and b = [1.0, 2.0, 3.0], solve Ax = b with both np.linalg.solve and scipy.linalg.solve(A, b, assume_a='pos'), and print np.allclose of the two results. Then print scipy.linalg.expm(np.zeros((2, 2))) and explain why NumPy alone cannot produce it.
Open the Python workspaceCheck your understanding
You have a 2-D float ndarray A and need its matrix exponential. Which statement is accurate?
- scipy.linalg.expm(A) computes it and hands back a NumPy ndarray, because NumPy ships no matrix exponential
- np.exp(A) already computes it, since the exponential of a matrix is defined elementwise
- scipy.linalg.expm returns a SciPy matrix object, so you must call np.asarray on it before using NumPy operations
- You must first convert A with scipy.asarray, because SciPy routines reject NumPy arrays
Show answer
expm evaluates the series sum A**k / k! (in practice via scaling-and-squaring with Pade approximants), a whole-matrix algorithm that NumPy does not provide. np.exp is a ufunc that applies exp to each entry independently, which is a different matrix entirely except for diagonal input. And there is no separate SciPy array type in either direction, so no conversion is ever needed.