A scikit-learn-compatible implementation, with a JAX-compiled numerical core, of the analytical nonlinear shrinkage covariance estimator from Ledoit & Wolf (2020), Analytical nonlinear shrinkage of large-dimensional covariance matrices, The Annals of Statistics 48(5):3043-3065. The method and its theoretical guarantees are the original authors' work; see Credits.
A sample covariance matrix estimated from limited data is a poor input for downstream optimization. Its largest eigenvalues are biased upward and its smallest are biased downward, and portfolio or risk models built directly on it inherit that distortion. Shrinkage estimators correct this by pulling sample eigenvalues toward a more stable target.
sklearn.covariance.LedoitWolf implements the 2004 linear shrinkage
estimator: every eigenvalue moves toward the same grand mean, with one shared
shrinkage intensity for the whole spectrum. The 2020 analytical nonlinear
estimator implemented here computes a separate shrinkage intensity for each
eigenvalue, derived from a kernel estimate of the sample spectral density and
its Hilbert transform. The paper proves this is asymptotically optimal, in
the minimum-variance sense, among rotation-equivariant estimators as the
number of features grows with the number of samples. The authors recommend it
starting around 12 features; below that, results typically match linear
shrinkage.
pip install nlshrink-jaximport numpy as np
from nlshrink_jax import NonlinearShrinkage
X = np.random.default_rng(0).normal(size=(500, 20))
est = NonlinearShrinkage().fit(X)
est.covariance_ # (20, 20) arrayThe estimator follows sklearn.covariance.LedoitWolf's interface (fit,
covariance_, get_precision, score, get_params/set_params), so it
substitutes directly wherever a scikit-learn covariance estimator is
accepted:
from sklearn.covariance import LedoitWolf
from nlshrink_jax import NonlinearShrinkage
for estimator in [LedoitWolf(), NonlinearShrinkage()]:
estimator.fit(X)
print(type(estimator).__name__, estimator.score(X_holdout))A functional form is also available for callers that want a plain array back without the estimator object:
from nlshrink_jax import analytical_nonlinear_shrinkage
cov = analytical_nonlinear_shrinkage(X, ridge=1e-6)n_samples >= n_features, and the sample covariance must be as full-rank
as that many observations can support (rank min(n_samples - 1, n_features)).
Duplicated or exactly collinear rows in X violate the second condition
even when the first holds. Both cases raise ValueError rather than
returning a silently incorrect result. The paper's n_features > n_samples
regime uses a different branch of the estimator that this implementation
does not provide.
The Hilbert transform of the Epanechnikov kernel has two terms that cancel
for eigenvalues far apart relative to the shrinkage bandwidth. Measured
directly against a float64 reference, a naive float32 evaluation loses
essentially all useful precision there (relative error of order 1, not a
small correction). This implementation switches to a cancellation-free
asymptotic series above that threshold, holding float32 relative error near
1e-6 on the same inputs. tests/test_estimator.py checks the shipped
series against a float64 reference oracle on a matrix constructed to
trigger the cancellation, at the coarser, end-to-end covariance-estimate
level (bounded at 1e-3 to 3e-3 Frobenius-relative there, since eigenvalue
recomposition adds its own error on top of the transform).
The numerical core is JIT-compiled with a fixed matrix shape, so repeated calls with the same number of features do not retrigger compilation. Cost is O(p²) per call for p features.
benchmarks/compare_baselines.py compares this estimator against the
sample covariance and sklearn.covariance.LedoitWolf on synthetic data
with a known ground-truth covariance, reporting estimation error, realized
minimum-variance portfolio risk, and fit time.
benchmarks/dimension_sweep.py runs the same comparison on real market
data (Binance spot daily returns), sweeping the number of assets across the
paper's stated 12-16-feature threshold. It measures two things: whole-matrix
estimation error (where this estimator shows a stable ~70-75% win rate
against linear shrinkage across the whole range) and realized minimum-variance
portfolio risk (where, on this dataset, linear shrinkage wins instead, a
genuine and explained divergence between a whole-matrix accuracy metric and
an inverse-sensitive downstream one). See the script's own docstring for
the full result, the reasoning behind the two metrics disagreeing, and how
to reproduce it (pip install -e ".[benchmark]" first; it fetches data
over HTTP from Binance's public API on first run).
git clone https://github.com/f-lndt/nlshrink-jax
cd nlshrink-jax
pip install -e ".[test]"
pytestThe estimator and its theoretical results (the analytical shrinkage formula, the Hilbert-transform kernel construction, the optimality proof) are entirely the work of the original paper's authors:
Ledoit, O. & Wolf, M. (2020). Analytical nonlinear shrinkage of large-dimensional covariance matrices. The Annals of Statistics, 48(5), 3043-3065. https://doi.org/10.1214/19-AOS1921
@article{ledoit2020analytical,
author = {Ledoit, Olivier and Wolf, Michael},
title = {Analytical nonlinear shrinkage of large-dimensional covariance matrices},
journal = {The Annals of Statistics},
volume = {48},
number = {5},
pages = {3043--3065},
year = {2020},
doi = {10.1214/19-AOS1921}
}This package is an independent implementation: Frederik Landt built the JAX-backed numerical core, adapted it to the scikit-learn estimator API, packaged it, and maintains this repository. It is not affiliated with or endorsed by the paper's authors.
MIT