Run this notebook: Open in Colab Open in Kaggle

Plot Lw Vs Oas

============================= Ledoit-Wolf vs OAS estimation

The usual covariance maximum likelihood estimate can be regularized using shrinkage. Ledoit and Wolf proposed a close formula to compute the asymptotically optimal shrinkage parameter (minimizing a MSE criterion), yielding the Ledoit-Wolf covariance estimate.

Chen et al. [1]_ proposed an improvement of the Ledoit-Wolf shrinkage parameter, the OAS coefficient, whose convergence is significantly better under the assumption that the data are Gaussian.

This example, inspired from Chen’s publication [1]_, shows a comparison of the estimated MSE of the LW and OAS methods, using Gaussian distributed data.

… rubric :: References

… [1] “Shrinkage Algorithms for MMSE Covariance Estimation” Chen et al., IEEE Trans. on Sign. Proc., Volume 58, Issue 10, October 2010.

Imports for Ledoit-Wolf vs OAS Shrinkage Comparison

Both LedoitWolf and OAS compute optimal shrinkage coefficients analytically, but OAS converges faster for Gaussian data: Shrinkage covariance estimation regularizes the sample covariance by blending it with a scaled identity matrix: Sigma_shrunk = (1-alpha)Sigma_empirical + alphatrace(Sigma_empirical)/p * I. The Ledoit-Wolf formula minimizes the expected squared Frobenius-norm error to the true covariance asymptotically, while the OAS formula by Chen et al. provides a finite-sample improvement under the Gaussian assumption. For very small samples (n_samples << n_features), OAS typically achieves lower MSE because its estimator accounts for the fourth-moment structure of Gaussian distributions.

The AR(1) Toeplitz covariance structure tests estimator behavior with realistic correlation decay: The true covariance toeplitz(r ** np.arange(n_features)) with r=0.1 produces a nearly diagonal matrix with weak correlations – a favorable case for shrinkage toward identity since the target is close to the truth. The Cholesky decomposition cholesky(real_cov) serves as the coloring matrix that transforms iid standard Gaussian samples into samples with the desired covariance. Running 100 Monte Carlo repetitions at each sample size (6 to 30) and computing error_norm(real_cov) (Frobenius-norm distance to the true covariance) produces reliable MSE estimates with error bars, showing how both methods’ accuracy improves with more samples and how their shrinkage coefficients decrease as the empirical covariance becomes more reliable.

# Authors: The scikit-learn developers
# SPDX-License-Identifier: BSD-3-Clause

import matplotlib.pyplot as plt
import numpy as np
from scipy.linalg import cholesky, toeplitz

from sklearn.covariance import OAS, LedoitWolf

np.random.seed(0)
# %%
n_features = 100
# simulation covariance matrix (AR(1) process)
r = 0.1
real_cov = toeplitz(r ** np.arange(n_features))
coloring_matrix = cholesky(real_cov)

n_samples_range = np.arange(6, 31, 1)
repeat = 100
lw_mse = np.zeros((n_samples_range.size, repeat))
oa_mse = np.zeros((n_samples_range.size, repeat))
lw_shrinkage = np.zeros((n_samples_range.size, repeat))
oa_shrinkage = np.zeros((n_samples_range.size, repeat))
for i, n_samples in enumerate(n_samples_range):
    for j in range(repeat):
        X = np.dot(np.random.normal(size=(n_samples, n_features)), coloring_matrix.T)

        lw = LedoitWolf(store_precision=False, assume_centered=True)
        lw.fit(X)
        lw_mse[i, j] = lw.error_norm(real_cov, scaling=False)
        lw_shrinkage[i, j] = lw.shrinkage_

        oa = OAS(store_precision=False, assume_centered=True)
        oa.fit(X)
        oa_mse[i, j] = oa.error_norm(real_cov, scaling=False)
        oa_shrinkage[i, j] = oa.shrinkage_

# plot MSE
plt.subplot(2, 1, 1)
plt.errorbar(
    n_samples_range,
    lw_mse.mean(1),
    yerr=lw_mse.std(1),
    label="Ledoit-Wolf",
    color="navy",
    lw=2,
)
plt.errorbar(
    n_samples_range,
    oa_mse.mean(1),
    yerr=oa_mse.std(1),
    label="OAS",
    color="darkorange",
    lw=2,
)
plt.ylabel("Squared error")
plt.legend(loc="upper right")
plt.title("Comparison of covariance estimators")
plt.xlim(5, 31)

# plot shrinkage coefficient
plt.subplot(2, 1, 2)
plt.errorbar(
    n_samples_range,
    lw_shrinkage.mean(1),
    yerr=lw_shrinkage.std(1),
    label="Ledoit-Wolf",
    color="navy",
    lw=2,
)
plt.errorbar(
    n_samples_range,
    oa_shrinkage.mean(1),
    yerr=oa_shrinkage.std(1),
    label="OAS",
    color="darkorange",
    lw=2,
)
plt.xlabel("n_samples")
plt.ylabel("Shrinkage")
plt.legend(loc="lower right")
plt.ylim(plt.ylim()[0], 1.0 + (plt.ylim()[1] - plt.ylim()[0]) / 10.0)
plt.xlim(5, 31)

plt.show()