psd-covariance 1.0.5
Computes several covariance matrix estimators that ensure positive semi-definiteness (PSD).
License: MIT
Requires Python: >=3.7

# A Package for Postive Definite Covariance Matrix Estimation *(psd-covariance)*

[![PyPI](https://img.shields.io/pypi/v/psd-covariance.svg)](https://pypi.org/project/psd-covariance/)
![Python versions](https://img.shields.io/pypi/pyversions/psd-covariance.svg)
![License](https://img.shields.io/pypi/l/psd-covariance.svg)
[![Downloads](https://pepy.tech/badge/psd-covariance)](https://pepy.tech/project/psd-covariance)

## Introduction
`psd-covariance` provides fast and reliable tools for estimating **positive semidefinite (PSD) covariance matrices**, even in challenging high-dimensional or noisy settings. It is designed for practitioners and researchers who frequently encounter unstable, ill-conditioned, or non-PSD covariance matrices.
At the core of the package are the **Posterior Mean (PM)** and **Fixed-Trace Posterior Mean (PM-FT)** estimators based on Bayesian eigenvalue regularization, as introduced in Boudt (2025). 

Whether you are building portfolio optimization models, performing principal component analysis, or simply need a numerically robust PSD covariance estimate, `psd-covariance` provides a practical and user-friendly implementation of these state-of-the-art estimators.

## Quick Start
#### Installation
```python
pip install psd-covariance
```
#### Imports
```python
import pandas as pd
import numpy as np
from psd_covariance import PosteriorMeanEstimator
```

### Example 1: Transforming non-PSD matrices
We construct a non-PSD estimated covariance matrix with d=10, such that the smallest two eigenvalues are negative. We update these eigenvalues using the (fixed-trace) posterior mean estimators.
```python
eigvals = [-0.50014538, -0.47905291,  0.59185356,  0.69919373,  0.8262933,   0.89063562, 1.005641,    1.39291291,  1.66503103,  1.95878406]
# Create matrix
d = len(eigvals)
A = np.random.randn(d, d)
Q, _ = np.linalg.qr(A)
Sigma_tilde = Q @ np.diag(eigvals) @ Q.T
```
```python
# 1. Posterior Mean 
pm = PosteriorMeanEstimator(fixed_trace=False)
pm_result = pm.fit(Sigma_tilde, sigma=0.5)  
eigvals_pm = np.linalg.eigvalsh(pm_result.cov)

# 2. Posterior Mean Fixed Trace
ft = PosteriorMeanEstimator(fixed_trace=True)
ft_result = ft.fit(Sigma_tilde, sigma=0.5)
eigvals_ft = np.linalg.eigvalsh(ft_result.cov)
```

``` python
print("\nEigenvalues of posterior mean matrices:")
print("PM Estimator            :", eigvals_pm)
print("PM Estimator (Fixed Tr.):", eigvals_ft)
# Eigenvalues of cleaned matrices:
# PM Estimator            : [0.2625387  0.26678995 0.70412859 0.78083976 0.87984287 0.93304454 1.03263025 1.39704147 1.66581096 1.9588768 ]
# PM Estimator (Fixed Tr.): [0.21390763 0.21737141 0.5737001  0.63620176 0.71686614 0.76021306 0.84135212 1.13826203 1.35724629 1.5960264 ]
```

### Example 2: PM and PM-FT Estimation using Cross-Validation
We generate a covariance matrix with a Toeplitz structure with d=10 and we draw n=20 observations from a Normal distribution with mean 0. We compute PM and PM-FT estimators using 10-fold cross validation.
``` python
# Generate data
np.random.seed(0)
d = 10
n = 20
rho = 0.8
cov_matrix = np.fromfunction(lambda i, j: rho ** np.abs(i - j), (d, d))
X = np.random.multivariate_normal(np.zeros(d), cov_matrix, size=n)
S = sample_cov(X)

reg_range = np.linspace(0.01, 2.0, 150)

# PM cross validation
pm = PosteriorMeanEstimator(fixed_trace=False)
reg_pm = pm.cross_validate_sigma(X, reg_range)
print(reg_pm)
# 0.2504026845637584
result_pm = pm.fit(S, reg_pm)

# FT cross validation
ft = PosteriorMeanEstimator(fixed_trace=True)
reg_ft = ft.cross_validate_sigma(X, reg_range)
print(reg_ft)
# 0.2771140939597316
result_ft = ft.fit(S, reg_ft)
```

## Available Classes with Methods

### PosteriorMeanEstimator
Input for posterior mean estimation: a covariance matrix `Sigma_tilde : ndarray (dxd)` together with the regularization parameter `sigma : float`, and the optional argument `fixed_trace : bool` indicating whether the fixed-trace variant (PM-FT) should be applied.

- `PosteriorMeanEstimator(fixed_trace=False)`: constructor that selects PM (`False`) or PM-FT (`True`).
- `fit(Sigma_tilde, sigma)`: computes the posterior mean covariance estimator using the chosen regularization parameter. Returns the named tuple `PosteriorMeanResult(cov, inv, sigma)`, where `cov` is the estimated covariance matrix, `inv` its precision matrix, and `sigma` the regularization parameter used in the estimator.
- `cross_validate_sigma(X, regularization_range, n_splits=10, n_jobs=-1)`: selects the optimal regularization parameter by evaluating each value in `regularization_range` using K-fold predictive log-likelihood cross-validation. Here `X` is a data matrix, `regularization_range` is an array of candidate regularization values, `n_splits` sets the number of folds (default 10), and `n_jobs` controls parallel computing (`-1` uses all CPU cores). Method returns optimal regularization value.

### EigenvalueCleaning
Input for eigenvalue cleaning methods: a covariance matrix `cov : ndarray (dxd)` together with the optional parameters `epsilon : float` (replacement value for negative eigenvalues), `fixed_trace : bool` (whether to preserve the trace of the original matrix), and `positive_definite : bool` (whether to enforce strict positive definiteness).

- `threshold_negative(cov, fixed_trace=False)`. Sets all negative eigenvalues to zero.
- `replace_negative(cov, epsilon=1e-4, fixed_trace=False, positive_definite=False)` Replaces negative eigenvalues with `epsilon`. 
- `absolute_negative(cov, fixed_trace=False, positive_definite=False)`. Takes the absolute value of all eigenvalues.  

All eigenvalue cleaning methods return return the named tuple `CleanedCovariance(cov, inv)`, with `cov` the estimated covariance matrix and `inv` its precision matrix.

### ShrinkageEstimator
Input for shrinkage estimators: data matrix of observations `X : pandas.DataFrame (nxd)`.

- `linear_shrinkage(X)`: Ledoit & Wolf (2002) linear shrinkage. Returns covariance, precision, and shrinkage intensity `alpha`.
- `quadratic_inverse_shrinkage(X)`: Ledoit & Wolf (2022) quadratic inverse shrinkage. Returns covariance and precision matrices.

All shrinkage estimators return `ShrinkageResult(cov, inv, shrinkage)`, with `cov` the estimated covariance matrix, `inv` its precision matrix, and `shrinkage` the shrinkage intensity (only for linear shrinkage; `None` for QIS).


## Authors
This package is based on the paper *'Well-Conditioned Covariance Estimation via Bayesian
Eigenvalue Regularization'*, by Kris Boudt, Jesper Cremers, Kirill Dragun & Steven Vanduffel. The *psd-covariance* package is developed and maintained by Jesper Cremers.

## References
- Boudt, K., J. Cremers, K. Dragun, and S. Vanduffel (2025). Well-conditioned covariance estimation via bayesian eigenvalue regularization. *Working paper.*
- Ledoit, O. and M. Wolf (2004). Honey, I shrunk the sample covariance matrix. *The Journal of Portfolio Management 30 (4), 110-119.*
- Ledoit, O. and M. Wolf (2022). Quadratic shrinkage for large covariance matrices. *Bernoulli 28 (3), 1519-1547.*
- Rousseeuw, P. J. and G. Molenberghs (1993). Transformation of non positive semidefinite correlation matrices. *Communications in Statistics - Theory and Methods 22 (4), 965-984.*
- Parts of the ShrinkageEstimator class contain adapted code from Michael Wolf's reference implementation: `https://github.com/pald22/covShrinkage`.