Performance and SciPy internals

In optimisation and integration the Python callback often dominates runtime, and knowing the internals lets you remove it.

SciPy's heavy lifting runs in compiled Fortran, C and C++ (LAPACK, QUADPACK, HiGHS, PROPACK and others), but routines like quad, minimize and solve_ivp call your Python function repeatedly. With cheap integrands the interpreter overhead of those calls dominates.

Ways to cut it:

Vectorise callbacks: solve_ivp(..., vectorized=True), differential_evolution(..., vectorized=True), stats.bootstrap with an axis-aware statistic, quad_vec instead of looping over quad. One call then handles many points.
Parallel workers: differential_evolution, scipy.fft, cKDTree.query and others take workers= (processes or threads, or a map-like callable).
LowLevelCallable: quad, ndimage.generic_filter and a few others accept a compiled C function (from ctypes, cffi or Numba cfunc) wrapped in scipy.LowLevelCallable, skipping Python entirely.
Jacobians: analytic jac for minimize/root, and a sparse jac_sparsity for least_squares and stiff solve_ivp methods, reduce function calls dramatically on large problems.

For linear algebra, scipy.linalg.get_blas_funcs and get_lapack_funcs return the raw routine for a dtype, and options like overwrite_a=True and check_finite=False skip copies and checks in hot loops (only when the inputs are known to be clean).

Thread counts of the underlying BLAS are controlled by environment variables or threadpoolctl. Many SciPy functions also support the Python array API standard, accepting CuPy or PyTorch arrays when SCIPY_ARRAY_API=1 is set.

Example

import numpy as np
from scipy import integrate, linalg

a = np.linspace(1, 5, 200)
calls = 0
def f(x, a):
    global calls
    calls += 1
    return np.exp(-a * x ** 2)

loop = [integrate.quad(f, 0, 1, args=(ai,))[0] for ai in a]   # 200 separate integrals
print(calls)                                   # 4200 Python calls
calls = 0
vec, err = integrate.quad_vec(lambda x: f(x, a), 0, 1)        # one vector-valued integral
print(calls, np.allclose(loop, vec))           # 63 True

A = np.random.default_rng(0).normal(size=(200, 200))
gesv, = linalg.get_lapack_funcs(("gesv",), (A,))
lu, piv, x, info = gesv(A, np.ones(200))       # raw LAPACK call
print(info, np.allclose(A @ x, 1))             # 0 True

Remember

Python callbacks often dominate; vectorise them or use LowLevelCallable
Use workers= and analytic or sparse Jacobians on large problems
check_finite=False and overwrite_a=True cut overhead in hot loops

Watch out for

Parallelising with many processes while BLAS also uses all cores, causing oversubscription.

Check yourself

Why can quad on a cheap integrand be dominated by Python overhead?

Official docs: docs.scipy.org/doc/scipy/reference/ccallback.html

Example

Remember

Watch out for

Check yourself

Builds on

Unlocks