Performance#

The filter in statsmodels is compiled Cython calling BLAS. The time around it is spent in Python: building the matrices at each evaluation, finite differences for the gradient, storing the full filter history when only the likelihood is needed, and fitting series one at a time. The library moves all of this into Fortran and adds four changes to the core.

Small matrices#

Decision. gemm, gemv and chol_inv use inline loops when the operands are small (below about 10 × 10 × 10 for products and order 16 for the inverse), and BLAS or LAPACK otherwise.

Why. State space models are small; m and p are often below 10. At that size the cost of a call into OpenBLAS, which checks its threading and dispatches on the operation, exceeds the arithmetic. The inline loops halved the time of the local level and ARMA filters.

A copy on every call#

Decision. The array arguments of the numerical routines are declared contiguous.

Why. gfortran 15 copies an assumed-shape argument when it passes it on to a contiguous dummy, even when the array is contiguous at run time; a debug build reported millions of such copies during the tests. Declaring the dummies contiguous all the way down removed them and halved the time per filter step again. The interfaces users implement (update, fill and the parameter transforms) are left unchanged, so that existing models still compile.

The steady state and the likelihood#

The log likelihood is computed without storing the filter history, and in time-invariant models the filter switches to its steady state (The steady state). For estimation, which evaluates only the likelihood and the score, these two changes matter most.

Parallel fits#

Decision. fit_many fits independent models on OpenMP threads, with the GIL released in Python.

Why. Fitting many short series is a common use and parallelizes without coordination. L-BFGS-B keeps its state in the caller’s arrays and the library has no global state, so the routines are reentrant; each model is touched by one thread only. BLAS should then run single-threaded (OPENBLAS_NUM_THREADS=1).

Numbers#

From example/bench.f90 and bench/bench_statsmodels.py, time for the log likelihood at n = 100,000 in milliseconds, before and after these changes, against statsmodels’ Cython filter:

model

before

after

statsmodels

local level (n = 1,000,000)

728

49

415

level and trigonometric seasonal

141

9.4

251

8 series, 3 factors

180

8.5

118

ARMA(2, 1)

79

5.2

44

Problems found#

  • An address sanitizer run found a heap overflow in collapse_observations: gfortran at -O2 sized an inlined matmul wrongly in an assignment whose shape changed between periods. The array is now allocated explicitly. The tests run under debug, release, OpenMP and address-sanitizer builds.

  • gfortran 15 at -O3 fails with an internal compiler error on names = mb%model%param_names() for a polymorphic model; a helper routine avoids it (model_names in statespace_capi_models).