[FEAT] Deprecate Numba

[Olivier Sprangers (elephaint)]

This PR removes numba as a runtime dependency by migrating all @njit-decorated hot loops to C++ via pybind11 and Eigen. The affected modules are CES, GARCH, MFLES, SES (intermittent demand helpers), and TBATS. Each module gets its own .cpp file under src/, compiled into the existing statsforecast._lib extension. A Nelder-Mead correctness fix and several robustness improvements are included.

Performance

The C++ implementations replace numba JIT compilation, eliminating:

• ~2-5 second cold-start penalty on first call (numba compilation)
• Runtime dependency on numba (22.8 MB) and llvmlite (161.0 MB), saving ~184 MB of install size
• The new C++ modules add only ~239 KB to the existing _lib.so (577 KB → 816 KB)

Benchmark Results

See experiments/numba_removal/benchmark.py for detailed comparisons.

All benchmarks run on series of length 500 (seed=42). Output equivalence: 88/88 PASS across all functions and series lengths.

Hot-path speedups (n=500, numba JIT-warmed):

Function	Series type	Numba (ms)	C++ (ms)	Speedup
garch_sigma2	garch	0.033	0.002	15.2x
garch_loglik	garch	0.035	0.005	7.8x
get_basis	seasonal	0.042	0.002	18.0x
cescalc	seasonal	0.124	0.035	3.5x
calcFaster (TBATS)	seasonal	0.086	0.028	3.1x
makeFMatrix (TBATS)	seasonal	0.005	0.002	3.1x
siegel_medians	seasonal	1.617	1.276	1.3x
ses_sse	all	0.001	0.001	~1.0x
ses_forecast	all	0.001	0.002	~0.7x
expand_demand/intervals	intermittent	0.001	0.001	~1.0x

SES and expand functions show parity — the operations are too lightweight for C++ call overhead to pay off at n=500. GARCH and MFLES get_basis show the largest gains (15-20x).

Full-pipeline speedups (optimizer + fitted values):

Pipeline	Numba (ms)	C++ (ms)	Speedup
GARCH fit (SLSQP + sigma2)	13.28	7.29	1.8x
CES fit (cescalc + fitted)	0.125	0.034	3.7x
TBATS fit (calcFaster + fitted)	0.084	0.026	3.3x
SES fit (optimize alpha + fitted)	0.064	0.066	~1.0x
Croston (SES + expand)	0.014	0.014	~1.0x

Scaling with series length (speedup at various n):

Function	n=50	n=100	n=500	n=1000	n=5000	n=10000
garch_sigma2	5.1x	7.8x	14.9x	17.3x	20.0x	20.4x
garch_loglik	3.5x	5.1x	8.5x	9.0x	9.7x	9.8x
cescalc	2.5x	3.0x	3.8x	3.9x	4.0x	4.0x
calcFaster	2.4x	2.8x	3.2x	3.3x	3.4x	3.2x
get_basis	50.2x	43.4x	16.4x	10.0x	3.4x	2.6x
siegel_medians	1.9x	1.4x	1.3x	1.3x	1.3x	1.3x
ses_sse	0.8x	0.8x	0.9x	1.0x	1.0x	1.0x

GARCH speedup increases with series length (up to 20x at n=10K) due to reduced Python↔numba dispatch overhead amortized over more iterations. get_basis shows the inverse: the numba JIT overhead dominates at small n (50x at n=50), while the algorithmic work converges at larger n. All scaling benchmarks pass equivalence checks (56/56).

Breaking Changes

numba removed from production dependencies. All @njit-decorated functions in ces.py, garch.py, mfles.py, models.py, tbats.py, and utils.py have been replaced with C++ (pybind11 + Eigen) implementations. Users who import numba transitively through statsforecast will need to add it to their own dependencies.

Removed public symbols from statsforecast.utils:

• CACHE, NOGIL, restrict_to_bounds — these were numba infrastructure. Any downstream code importing them will break.

Environment variables no longer recognized:

• NIXTLA_NUMBA_RELEASE_GIL — previously controlled @njit(nogil=...). GIL release is now handled at the pybind11 binding level.
• NIXTLA_NUMBA_CACHE — previously controlled @njit(cache=...). No longer applicable; C++ is compiled at install time.

Behavioral Changes

GARCH optimizer lower bounds changed from 0 to 1e-8. In garch.py, the SLSQP bounds were tightened from (0, None) to (1e-8, None). This prevents degenerate models where coefficients are exactly zero, but may produce slightly different fitted parameters for series where coefficients previously converged near zero.

CES optimization now uses float64 initial parameters. In ces.py, x0 changed from np.float32 to np.float64. The Nelder-Mead optimizer now operates in double precision throughout, which may produce slightly different optimal parameters (generally more accurate).

Nelder-Mead shrink step corrected. In nelder_mead.h, the shrink step previously always skipped row 0 (for j in range(1, ...)), which was only correct when the best point happened to be at index 0. The fix skips the actual best_idx. This is a correctness improvement that may change optimization results for any model using Nelder-Mead (CES, ETS, Theta).

Numerical Precision Notes

CES results have small float32 rounding differences. The C++ compilation may reorder or fuse float32 operations differently from numba's LLVM backend (e.g. FMA instructions), producing diffs at the last bit of float32 precision. Measured on AirPassengers with fixed parameters:

Path	Max absolute diff
Nonseasonal / Simple seasonal states	0.0 (exact)
Partial seasonal states	~6e-5
Full seasonal states	~2e-4
Forecasts (all season types)	up to ~3e-3

The forecast diffs are amplified by the float64→float32 round-trip in the cesforecast Python wrapper. All tests pass at decimal=2 or better.

Test Coverage Gaps

The following C++ modules lack dedicated unit tests and are only exercised indirectly through integration tests:

• src/garch.cpp — tested via test_garch.py (11 tests, model-level)
• src/mfles.cpp — tested via test_mfles.py (2 tests, model-level)
• src/ses.cpp — tested via test_models.py (167 tests, model-level)
• src/tbats.cpp — tested via test_tbats.py (1 test, model-level)

The experiments/numba_removal/benchmark.py contains correctness checks comparing old numba vs new C++ output, but these don't run in CI.

[codspeed]

Merging this PR will improve performance by ×4.6

⚡ 2 improved benchmarks
✅ 4 untouched benchmarks
🆕 32 new benchmarks

Performance Changes

	Benchmark	BASE	HEAD	Efficiency
🆕	test_efficiency[NaNModel]	N/A	57.5 µs	N/A
🆕	test_efficiency[MFLES]	N/A	7.9 ms	N/A
🆕	test_efficiency[Naive]	N/A	65.1 µs	N/A
🆕	test_efficiency[ARIMA]	N/A	3.9 ms	N/A
🆕	test_efficiency[OptimizedTheta]	N/A	5.4 ms	N/A
🆕	test_efficiency[RandomWalkWithDrift]	N/A	94.3 µs	N/A
🆕	test_efficiency[MSTL]	N/A	1.1 s	N/A
🆕	test_efficiency[Theta]	N/A	4.8 ms	N/A
🆕	test_efficiency[SeasonalExponentialSmoothingOptimized]	N/A	15.8 ms	N/A
🆕	test_efficiency[AutoMFLES]	N/A	130.5 ms	N/A
🆕	test_efficiency[AutoRegressive]	N/A	154.6 ms	N/A
⚡	test_efficiency[AutoCES]	13.8 ms	3 ms	×4.6
🆕	test_efficiency[ARCH]	N/A	19.9 ms	N/A
🆕	test_efficiency[AutoTBATS]	N/A	8.6 s	N/A
🆕	test_efficiency[SimpleExponentialSmoothing]	N/A	92.7 µs	N/A
🆕	test_efficiency[CrostonClassic]	N/A	217.5 µs	N/A
🆕	test_efficiency[ADIDA]	N/A	909.6 µs	N/A
🆕	test_efficiency[TSB]	N/A	157.9 µs	N/A
🆕	test_efficiency[CrostonOptimized]	N/A	1.2 ms	N/A
🆕	test_efficiency[IMAPA]	N/A	932 µs	N/A
...	...	...	...	...

ℹ️ Only the first 20 benchmarks are displayed. Go to the app to view all benchmarks.

---

_{Comparing feat/deprecate_numba (68314af) with main (65569e1)}

[Saul (nasaul)]

Overall this is an excellent PR which get's an increase speed, great work!
However before merging there are things that need to be addressed:

• Remove numba_cache.ipynb from the repo and from the documentation.
• Remove numba references from
  • README.md
  • nbs/docs/experiments/AmazonStatsForecast.ipynb
  • nbs/docs/experiments/AutoArima_vs_Prophet.ipynb
  • nbs/docs/getting-started/*
  • nbs/docs/tutorials/MultipleSeasonalities.ipynb
• Decide what to do with the experimetns/numba_removal folder, my suggestion would be to remove it, but we could also add a proper experiment folder.
• Decide what to do with the experiments that cite numba.

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[Saul (nasaul)]

LGTM. Should we re run some of the experiments with this new improvement?

[Olivier Sprangers (elephaint)]

LGTM. Should we re run some of the experiments with this new improvement?

Which experiments?