API Reference

Core Classes

shrinkr

Shrinkr package.

CovarianceEstimator

Bases: BaseEstimator

Covariance matrix estimator with optional shrinkage.

Wraps several shrinkage methods behind a scikit-learn-compatible fit / predict interface.

Parameters:

Name	Type	Description	Default
method	str	Shrinkage method to apply. One of 'empirical', 'lw_linear', 'lw_analytical', 'oas', or their 'ref_*' reference counterparts. Default is 'empirical'.	'empirical'
tol	float	Eigenvalue threshold passed to eigenvalue-based methods. Default is 1e-8.	1e-08

See Also

shrinkr.functional.lw_linear Ledoit-Wolf Linear Shrinkage

shrinkr.functional.lw_analytical Ledoit-Wolf Analytical Shrinkage

shrinkr.functional.oas Oracle Approximating Shrinkage

Attributes:

Name	Type	Description
is_fitted_	bool	True after fit has been called successfully.
data_	ndarray	The data passed to fit.
cov_	ndarray or None	Raw sample covariance matrix. Set only for methods that require it.
shrunk_cov_	ndarray or None	Shrinkage-regularized covariance matrix produced by fit.

Source code in shrinkr/cov.py

class CovarianceEstimator(BaseEstimator):
    """Covariance matrix estimator with optional shrinkage.

    Wraps several shrinkage methods behind a scikit-learn-compatible
    ``fit`` / ``predict`` interface.

    Parameters
    ----------
    method : str, optional
        Shrinkage method to apply. One of ``'empirical'``, ``'lw_linear'``,
        ``'lw_analytical'``, ``'oas'``, or their ``'ref_*'`` reference
        counterparts. Default is ``'empirical'``.
    tol : float, optional
        Eigenvalue threshold passed to eigenvalue-based methods. Default is 1e-8.


    See Also
    --------
    [`shrinkr.functional.lw_linear`][]
        Ledoit-Wolf Linear Shrinkage

    [`shrinkr.functional.lw_analytical`][]
        Ledoit-Wolf Analytical Shrinkage

    [`shrinkr.functional.oas`][]
        Oracle Approximating Shrinkage


    Attributes
    ----------
    is_fitted_ : bool
        True after [`fit`][shrinkr.cov.CovarianceEstimator.fit] has been called successfully.
    data_ : np.ndarray
        The data passed to [`fit`][shrinkr.cov.CovarianceEstimator.fit].
    cov_ : np.ndarray or None
        Raw sample covariance matrix. Set only for methods that require it.
    shrunk_cov_ : np.ndarray or None
        Shrinkage-regularized covariance matrix produced by [`fit`][shrinkr.cov.CovarianceEstimator.fit].
    """

    def __init__(self, *, param=1, method: str = "empirical", tol: float = 1e-8):
        self.param = param
        self.tol = tol

        if method not in METHODS:
            raise ValueError(f"Method '{method}' must be one of: {METHODS}")
        self.method = method

        self.is_fitted_: bool = False
        self.data_: np.ndarray | None = None
        self.cov_: np.ndarray | None = None
        self.shrunk_cov_: np.ndarray | None = None

    def fit(self, X: np.ndarray, y=None):
        """Compute the (shrunk) covariance matrix from data.

        Parameters
        ----------
        X : np.ndarray
            Data matrix of shape (n_samples, n_features).
        y : ignored
            Present for API compatibility.

        Returns
        -------
        self : CovarianceEstimator
            Fitted estimator.
        """
        if X.ndim != 2:
            raise ValueError(f"Expected 2D array, got {X.ndim}D array instead.")

        self.data_ = X
        n, p = X.shape

        ## Methods requiring only data
        if self.method == "lw_linear":
            self.shrunk_cov_, _ = lw_linear(X, assume_centered=False)
            self.is_fitted_ = True
            return self

        if self.method == "ref_lw_linear":
            self.shrunk_cov_, _ = ref_lw_linear(X, assume_centered=False)
            self.is_fitted_ = True
            return self

        ## Methods requiring standard covariance
        self.cov_ = np.cov(X, rowvar=False, bias=True)

        if self.method == "empirical":
            self.shrunk_cov_ = self.cov_  # Empirical usually means unshrunk
            self.is_fitted_ = True
            return self

        if self.method == "oas":
            self.shrunk_cov_, _ = oas(self.cov_, n, p)
            self.is_fitted_ = True
            return self

        if self.method == "ref_oas":
            self.shrunk_cov_, _ = ref_oas(self.cov_, n, p)
            self.is_fitted_ = True
            return self

        ## Methods requiring eigenvalue decomposition
        evals, U = np.linalg.eigh(self.cov_)

        if self.method == "lw_analytical":
            shrunk_evals = lw_analytical(evals, n, p, self.tol)
            self.shrunk_cov_ = U @ np.diag(shrunk_evals) @ U.T
        elif self.method == "ref_lw_analytical":
            shrunk_evals = ref_lw_analytical(evals, n, self.tol)
            self.shrunk_cov_ = U @ np.diag(shrunk_evals) @ U.T
        else:
            raise ValueError(f"Implementation for method '{self.method}' not found.")

        self.is_fitted_ = True
        return self

    def predict(self, X):
        """Return the fitted covariance matrix.

        Parameters
        ----------
        X : ignored
            Present for API compatibility.

        Returns
        -------
        np.ndarray
            The shrinkage-regularized covariance matrix, or the raw sample
            covariance if no shrinkage method produced a result.
        """
        if not self.is_fitted_:
            raise ValueError("This CovarianceEstimator instance is not fitted yet.")

        if self.shrunk_cov_ is not None:
            return self.shrunk_cov_
        return self.cov_

    def fit_predict(self, X, y=None):
        """Fit and immediately return the covariance matrix.

        Equivalent to calling [`fit`][shrinkr.cov.CovarianceEstimator.fit] followed by [`predict`][shrinkr.cov.CovarianceEstimator.predict].

        Parameters
        ----------
        X : np.ndarray
            Data matrix of shape (n_samples, n_features).
        y : ignored
            Present for API compatibility.

        Returns
        -------
        np.ndarray
            The shrinkage-regularized covariance matrix.
        """
        return self.fit(X, y).predict(X)

fit(X, y=None)

Compute the (shrunk) covariance matrix from data.

Parameters:

Name	Type	Description	Default
X	ndarray	Data matrix of shape (n_samples, n_features).	required
y	ignored	Present for API compatibility.	None

Returns:

Name	Type	Description
self	CovarianceEstimator	Fitted estimator.

Source code in shrinkr/cov.py

def fit(self, X: np.ndarray, y=None):
    """Compute the (shrunk) covariance matrix from data.

    Parameters
    ----------
    X : np.ndarray
        Data matrix of shape (n_samples, n_features).
    y : ignored
        Present for API compatibility.

    Returns
    -------
    self : CovarianceEstimator
        Fitted estimator.
    """
    if X.ndim != 2:
        raise ValueError(f"Expected 2D array, got {X.ndim}D array instead.")

    self.data_ = X
    n, p = X.shape

    ## Methods requiring only data
    if self.method == "lw_linear":
        self.shrunk_cov_, _ = lw_linear(X, assume_centered=False)
        self.is_fitted_ = True
        return self

    if self.method == "ref_lw_linear":
        self.shrunk_cov_, _ = ref_lw_linear(X, assume_centered=False)
        self.is_fitted_ = True
        return self

    ## Methods requiring standard covariance
    self.cov_ = np.cov(X, rowvar=False, bias=True)

    if self.method == "empirical":
        self.shrunk_cov_ = self.cov_  # Empirical usually means unshrunk
        self.is_fitted_ = True
        return self

    if self.method == "oas":
        self.shrunk_cov_, _ = oas(self.cov_, n, p)
        self.is_fitted_ = True
        return self

    if self.method == "ref_oas":
        self.shrunk_cov_, _ = ref_oas(self.cov_, n, p)
        self.is_fitted_ = True
        return self

    ## Methods requiring eigenvalue decomposition
    evals, U = np.linalg.eigh(self.cov_)

    if self.method == "lw_analytical":
        shrunk_evals = lw_analytical(evals, n, p, self.tol)
        self.shrunk_cov_ = U @ np.diag(shrunk_evals) @ U.T
    elif self.method == "ref_lw_analytical":
        shrunk_evals = ref_lw_analytical(evals, n, self.tol)
        self.shrunk_cov_ = U @ np.diag(shrunk_evals) @ U.T
    else:
        raise ValueError(f"Implementation for method '{self.method}' not found.")

    self.is_fitted_ = True
    return self

fit_predict(X, y=None)

Fit and immediately return the covariance matrix.

Equivalent to calling fit followed by predict.

Parameters:

Name	Type	Description	Default
X	ndarray	Data matrix of shape (n_samples, n_features).	required
y	ignored	Present for API compatibility.	None

Returns:

Type	Description
ndarray	The shrinkage-regularized covariance matrix.

Source code in shrinkr/cov.py

def fit_predict(self, X, y=None):
    """Fit and immediately return the covariance matrix.

    Equivalent to calling [`fit`][shrinkr.cov.CovarianceEstimator.fit] followed by [`predict`][shrinkr.cov.CovarianceEstimator.predict].

    Parameters
    ----------
    X : np.ndarray
        Data matrix of shape (n_samples, n_features).
    y : ignored
        Present for API compatibility.

    Returns
    -------
    np.ndarray
        The shrinkage-regularized covariance matrix.
    """
    return self.fit(X, y).predict(X)

predict(X)

Return the fitted covariance matrix.

Parameters:

Name	Type	Description	Default
X	ignored	Present for API compatibility.	required

Returns:

Type	Description
ndarray	The shrinkage-regularized covariance matrix, or the raw sample covariance if no shrinkage method produced a result.

Source code in shrinkr/cov.py

def predict(self, X):
    """Return the fitted covariance matrix.

    Parameters
    ----------
    X : ignored
        Present for API compatibility.

    Returns
    -------
    np.ndarray
        The shrinkage-regularized covariance matrix, or the raw sample
        covariance if no shrinkage method produced a result.
    """
    if not self.is_fitted_:
        raise ValueError("This CovarianceEstimator instance is not fitted yet.")

    if self.shrunk_cov_ is not None:
        return self.shrunk_cov_
    return self.cov_

LinearDiscriminantAnalysis

Bases: BaseEstimator

Binary Linear Discriminant Analysis with pluggable covariance shrinkage.

Fits a two class LDA model and classifies by the log-posterior ratio. The pooled covariance can be estimated with any method supported by CovarianceEstimator or specialized LDA shrinkage deal.

Parameters:

Name	Type	Description	Default
method	str	Estimator used to shrink (compute) the pooled within-class covariance. See shrinkr.CovarianceEstimator for allowed values. Additionally we allow deal and ref_deal. Default is empirical.	'empirical'

See Also

shrinkr.CovarianceEstimator CovarianceEstimator class with non directional shrinkage methods

shrinkr.functional.lw_analytical Deterministic Equivalent Adjusted LDA (directional shrinakge)

Attributes:

Name	Type	Description
covariance_estimator_	CovarianceEstimator	Instance of the Covariance Estimator for shrinkage LDA.
classes_	np.ndarray of shape (2,)	The two class labels seen during fit.
priors_	np.ndarray of shape (2,)	Class prior probabilities estimated from the training data.
means_	np.ndarray of shape (2, n_features)	Per-class sample means.
covariance_	np.ndarray of shape (n_features, n_features)	Pooled within-class covariance matrix (after shrinkage if applicable).
precision_	np.ndarray of shape (n_features, n_features)	Pseudo-inverse of covariance_.
coef_	np.ndarray of shape (1, n_features)	Linear discriminant direction.
intercept_	np.ndarray of shape (1,)	Decision boundary offset.
is_fitted_	bool	True after fit has been called successfully.

Source code in shrinkr/lda.py

class LinearDiscriminantAnalysis(BaseEstimator):
    """Binary Linear Discriminant Analysis with pluggable covariance shrinkage.

    Fits a two class LDA model and classifies by the
    log-posterior ratio. The pooled covariance can be estimated with any
    method supported by [`CovarianceEstimator`][shrinkr.CovarianceEstimator]
    or specialized LDA shrinkage ``deal``.

    Parameters
    ----------
    method : str, optional
        Estimator used to shrink (compute) the pooled within-class covariance.
        See [`shrinkr.CovarianceEstimator`][] for allowed values.
        Additionally we allow ``deal`` and ``ref_deal``.
        Default is ``empirical``.

    See Also
    --------
    [`shrinkr.CovarianceEstimator`][]
        CovarianceEstimator class with non directional shrinkage methods

    [`shrinkr.functional.lw_analytical`][]
        Deterministic Equivalent Adjusted LDA (directional shrinakge)


    Attributes
    ----------
    covariance_estimator_ : CovarianceEstimator
        Instance of the Covariance Estimator for shrinkage LDA.
    classes_ : np.ndarray of shape (2,)
        The two class labels seen during [`fit`][shrinkr.lda.LinearDiscriminantAnalysis.fit].
    priors_ : np.ndarray of shape (2,)
        Class prior probabilities estimated from the training data.
    means_ : np.ndarray of shape (2, n_features)
        Per-class sample means.
    covariance_ : np.ndarray of shape (n_features, n_features)
        Pooled within-class covariance matrix (after shrinkage if applicable).
    precision_ : np.ndarray of shape (n_features, n_features)
        Pseudo-inverse of `covariance_`.
    coef_ : np.ndarray of shape (1, n_features)
        Linear discriminant direction.
    intercept_ : np.ndarray of shape (1,)
        Decision boundary offset.
    is_fitted_ : bool
        True after [`fit`][shrinkr.lda.LinearDiscriminantAnalysis.fit] has been called successfully.
    """

    def __init__(self, method: str = "empirical"):

        self.method = method

        if self.method in LDA_ONLY_METHODS:
            self.covariance_estimator_ = CovarianceEstimator(method="empirical")
        else:
            self.covariance_estimator_ = CovarianceEstimator(method=self.method)

        self.classes_: np.ndarray | None = None
        self.priors_: np.ndarray | None = None
        self.means_: np.ndarray | None = None
        self.covariance_: np.ndarray | None = None
        self.precision_: np.ndarray | None = None
        self.coef_: np.ndarray | None = None
        self.intercept_: np.ndarray | None = None
        self.is_fitted_: bool = False

    def fit(self, X: np.ndarray, y: np.ndarray):
        """Fit the LDA model on labelled training data.

        Parameters
        ----------
        X : np.ndarray of shape (n_samples, n_features)
            Training data.
        y : np.ndarray of shape (n_samples,)
            Binary class labels. Exactly two distinct values must be present.

        Returns
        -------
        self : LinearDiscriminantAnalysis
            Fitted estimator.

        Raises
        ------
        ValueError
            If ``X`` is not 2-D, ``X`` and ``y`` have different lengths, or
            ``y`` does not contain exactly two classes.
        """
        if X.ndim != 2:
            raise ValueError(f"Expected 2D array, got {X.ndim}D array instead.")
        if len(X) != len(y):
            raise ValueError("X and y must have the same number of samples.")

        self.classes_, y_indices = np.unique(y, return_inverse=True)

        if len(self.classes_) != 2:
            raise ValueError(f"Number of classes must be 2, but found {len(self.classes_)}.")

        n_samples, n_features = X.shape
        self.priors_ = np.bincount(y_indices) / n_samples
        self.means_ = np.zeros((2, n_features))

        X_centered = np.empty_like(X, dtype=float)
        for idx, cls in enumerate(self.classes_):
            class_mask = y == cls
            X_class = X[class_mask]
            self.means_[idx] = np.mean(X_class, axis=0)
            X_centered[class_mask] = X_class - self.means_[idx]

        self.covariance_estimator_.fit(X_centered)
        self.covariance_ = self.covariance_estimator_.predict(X_centered)

        if self.method in LDA_ONLY_METHODS:
            # TODO: Extend beyond just DEAL
            cov_evals, U = np.linalg.eigh(self.covariance_)

            Ut = U.T
            z1_vec = Ut @ self.means_[1]
            z0_vec = Ut @ self.means_[0]
            z_vec = z1_vec - z0_vec  # bcs linearity

            if self.method == "deal":
                adj_evals = deal(cov_evals, z_vec, n_samples - 2)
            elif self.method == "ref_deal":
                adj_evals = ref_deal(cov_evals, z_vec, n_samples - 2)
            else:
                raise ValueError("Method not implemented.")

            inv_adj_evals = 1.0 / adj_evals
            w = U @ (inv_adj_evals * z_vec)
            b = -0.5 * (
                np.sum(
                    z1_vec * z1_vec * inv_adj_evals
                )  # z1_vec.T @ np.diag(inv_adj_evals) @ z1_vec
                - np.sum(z0_vec * z0_vec * inv_adj_evals)
            ) + np.log(self.priors_[1] / self.priors_[0])
        else:
            self.precision_ = np.linalg.pinv(self.covariance_)
            mean_diff = self.means_[1] - self.means_[0]
            w = self.precision_ @ mean_diff
            b = -0.5 * (
                self.means_[1] @ self.precision_ @ self.means_[1]
                - self.means_[0] @ self.precision_ @ self.means_[0]
            ) + np.log(self.priors_[1] / self.priors_[0])

        self.coef_ = np.array([w])
        self.intercept_ = np.array([b])
        self.is_fitted_ = True

        return self

    def decision_function(self, X: np.ndarray):
        """Compute the log-odds score for the positive class.

        Parameters
        ----------
        X : np.ndarray of shape (n_samples, n_features)
            Samples to score.

        Returns
        -------
        np.ndarray of shape (n_samples,)
            Log-odds of the positive class for each sample.
        """
        if not self.is_fitted_:
            raise ValueError("This estimator is not fitted yet.")
        return (X @ self.coef_.T + self.intercept_).ravel()

    def predict(self, X: np.ndarray) -> np.ndarray:
        """Predict binary class labels.

        Parameters
        ----------
        X : np.ndarray of shape (n_samples, n_features)
            Samples to classify.

        Returns
        -------
        np.ndarray of shape (n_samples,)
            Predicted class label for each sample.
        """
        scores = self.decision_function(X)
        indices = (scores > 0).astype(int)
        return self.classes_[indices]

    def fit_predict(self, X, y) -> np.ndarray:
        """Fit and predict the predict classes.

        Equivalent to calling [`fit`][shrinkr.lda.LinearDiscriminantAnalysis.fit] followed by [`predict`][shrinkr.lda.LinearDiscriminantAnalysis.predict].

        Parameters
        ----------
        X : np.ndarray of shape (n_samples, n_features)
            Training data.
        y : np.ndarray of shape (n_samples,)
            Binary class labels. Exactly two distinct values must be present.

        Returns
        -------
        np.ndarray
            The shrinkage-regularized covariance matrix.
        """
        return self.fit(X, y).predict(X)

    def predict_proba(self, X: np.ndarray):
        """Estimate class probabilities using the logistic sigmoid.

        Parameters
        ----------
        X : np.ndarray of shape (n_samples, n_features)
            Samples to score.

        Returns
        -------
        np.ndarray of shape (n_samples, 2)
            Columns are ``[P(class 0), P(class 1)]`` for each sample.
        """
        scores = self.decision_function(X)
        prob_1 = 1 / (1 + np.exp(-scores))
        return np.vstack([1 - prob_1, prob_1]).T

decision_function(X)

Compute the log-odds score for the positive class.

Parameters:

Name	Type	Description	Default
X	np.ndarray of shape (n_samples, n_features)	Samples to score.	required

Returns:

Type	Description
np.ndarray of shape (n_samples,)	Log-odds of the positive class for each sample.

Source code in shrinkr/lda.py

def decision_function(self, X: np.ndarray):
    """Compute the log-odds score for the positive class.

    Parameters
    ----------
    X : np.ndarray of shape (n_samples, n_features)
        Samples to score.

    Returns
    -------
    np.ndarray of shape (n_samples,)
        Log-odds of the positive class for each sample.
    """
    if not self.is_fitted_:
        raise ValueError("This estimator is not fitted yet.")
    return (X @ self.coef_.T + self.intercept_).ravel()

fit(X, y)

Fit the LDA model on labelled training data.

Parameters:

Name	Type	Description	Default
X	np.ndarray of shape (n_samples, n_features)	Training data.	required
y	np.ndarray of shape (n_samples,)	Binary class labels. Exactly two distinct values must be present.	required

Returns:

Name	Type	Description
self	LinearDiscriminantAnalysis	Fitted estimator.

Raises:

Type	Description
ValueError	If X is not 2-D, X and y have different lengths, or y does not contain exactly two classes.

Source code in shrinkr/lda.py

def fit(self, X: np.ndarray, y: np.ndarray):
    """Fit the LDA model on labelled training data.

    Parameters
    ----------
    X : np.ndarray of shape (n_samples, n_features)
        Training data.
    y : np.ndarray of shape (n_samples,)
        Binary class labels. Exactly two distinct values must be present.

    Returns
    -------
    self : LinearDiscriminantAnalysis
        Fitted estimator.

    Raises
    ------
    ValueError
        If ``X`` is not 2-D, ``X`` and ``y`` have different lengths, or
        ``y`` does not contain exactly two classes.
    """
    if X.ndim != 2:
        raise ValueError(f"Expected 2D array, got {X.ndim}D array instead.")
    if len(X) != len(y):
        raise ValueError("X and y must have the same number of samples.")

    self.classes_, y_indices = np.unique(y, return_inverse=True)

    if len(self.classes_) != 2:
        raise ValueError(f"Number of classes must be 2, but found {len(self.classes_)}.")

    n_samples, n_features = X.shape
    self.priors_ = np.bincount(y_indices) / n_samples
    self.means_ = np.zeros((2, n_features))

    X_centered = np.empty_like(X, dtype=float)
    for idx, cls in enumerate(self.classes_):
        class_mask = y == cls
        X_class = X[class_mask]
        self.means_[idx] = np.mean(X_class, axis=0)
        X_centered[class_mask] = X_class - self.means_[idx]

    self.covariance_estimator_.fit(X_centered)
    self.covariance_ = self.covariance_estimator_.predict(X_centered)

    if self.method in LDA_ONLY_METHODS:
        # TODO: Extend beyond just DEAL
        cov_evals, U = np.linalg.eigh(self.covariance_)

        Ut = U.T
        z1_vec = Ut @ self.means_[1]
        z0_vec = Ut @ self.means_[0]
        z_vec = z1_vec - z0_vec  # bcs linearity

        if self.method == "deal":
            adj_evals = deal(cov_evals, z_vec, n_samples - 2)
        elif self.method == "ref_deal":
            adj_evals = ref_deal(cov_evals, z_vec, n_samples - 2)
        else:
            raise ValueError("Method not implemented.")

        inv_adj_evals = 1.0 / adj_evals
        w = U @ (inv_adj_evals * z_vec)
        b = -0.5 * (
            np.sum(
                z1_vec * z1_vec * inv_adj_evals
            )  # z1_vec.T @ np.diag(inv_adj_evals) @ z1_vec
            - np.sum(z0_vec * z0_vec * inv_adj_evals)
        ) + np.log(self.priors_[1] / self.priors_[0])
    else:
        self.precision_ = np.linalg.pinv(self.covariance_)
        mean_diff = self.means_[1] - self.means_[0]
        w = self.precision_ @ mean_diff
        b = -0.5 * (
            self.means_[1] @ self.precision_ @ self.means_[1]
            - self.means_[0] @ self.precision_ @ self.means_[0]
        ) + np.log(self.priors_[1] / self.priors_[0])

    self.coef_ = np.array([w])
    self.intercept_ = np.array([b])
    self.is_fitted_ = True

    return self

fit_predict(X, y)

Fit and predict the predict classes.

Equivalent to calling fit followed by predict.

Parameters:

Name	Type	Description	Default
X	np.ndarray of shape (n_samples, n_features)	Training data.	required
y	np.ndarray of shape (n_samples,)	Binary class labels. Exactly two distinct values must be present.	required

Returns:

Type	Description
ndarray	The shrinkage-regularized covariance matrix.

Source code in shrinkr/lda.py

def fit_predict(self, X, y) -> np.ndarray:
    """Fit and predict the predict classes.

    Equivalent to calling [`fit`][shrinkr.lda.LinearDiscriminantAnalysis.fit] followed by [`predict`][shrinkr.lda.LinearDiscriminantAnalysis.predict].

    Parameters
    ----------
    X : np.ndarray of shape (n_samples, n_features)
        Training data.
    y : np.ndarray of shape (n_samples,)
        Binary class labels. Exactly two distinct values must be present.

    Returns
    -------
    np.ndarray
        The shrinkage-regularized covariance matrix.
    """
    return self.fit(X, y).predict(X)

predict(X)

Predict binary class labels.

Parameters:

Name	Type	Description	Default
X	np.ndarray of shape (n_samples, n_features)	Samples to classify.	required

Returns:

Type	Description
np.ndarray of shape (n_samples,)	Predicted class label for each sample.

Source code in shrinkr/lda.py

def predict(self, X: np.ndarray) -> np.ndarray:
    """Predict binary class labels.

    Parameters
    ----------
    X : np.ndarray of shape (n_samples, n_features)
        Samples to classify.

    Returns
    -------
    np.ndarray of shape (n_samples,)
        Predicted class label for each sample.
    """
    scores = self.decision_function(X)
    indices = (scores > 0).astype(int)
    return self.classes_[indices]

predict_proba(X)

Estimate class probabilities using the logistic sigmoid.

Parameters:

Name	Type	Description	Default
X	np.ndarray of shape (n_samples, n_features)	Samples to score.	required

Returns:

Type	Description
np.ndarray of shape (n_samples, 2)	Columns are [P(class 0), P(class 1)] for each sample.

Source code in shrinkr/lda.py

def predict_proba(self, X: np.ndarray):
    """Estimate class probabilities using the logistic sigmoid.

    Parameters
    ----------
    X : np.ndarray of shape (n_samples, n_features)
        Samples to score.

    Returns
    -------
    np.ndarray of shape (n_samples, 2)
        Columns are ``[P(class 0), P(class 1)]`` for each sample.
    """
    scores = self.decision_function(X)
    prob_1 = 1 / (1 + np.exp(-scores))
    return np.vstack([1 - prob_1, prob_1]).T

---

Functional API

shrinkr.functional

Functional C implementation of shrinkages and more.

accuracy(y, y_pred)

Classification accuracy.

Parameters:

Name	Type	Description	Default
y	ndarray	True class labels (1D integer array).	required
y_pred	ndarray	Predicted class labels (1D integer array).	required

Returns:

Type	Description
float	Fraction of correctly classified samples, in the range [0, 1].

Source code in shrinkr/functional/_losses.py

def accuracy(y: np.ndarray, y_pred: np.ndarray) -> float:
    """Classification accuracy.

    Parameters
    ----------
    y : np.ndarray
        True class labels (1D integer array).
    y_pred : np.ndarray
        Predicted class labels (1D integer array).

    Returns
    -------
    float
        Fraction of correctly classified samples, in the range [0, 1].
    """
    if y.ndim != 1:
        raise ValueError("y must be a 1D array")
    if y_pred.ndim != 1:
        raise ValueError("y_pred must be a 1D array")
    if y.shape != y_pred.shape:
        raise ValueError("y and y_pred must have the same shape")

    return float(np.sum(y == y_pred) / y.shape[0])

deal(evals, z_vec, n_eff, gamma_min=0.02, gamma_max=100, base_shrinkage='lw_analytical', surrogate_shrinkage='lw_analytical', eps=1e-08, **kwargs)

DEAL (Deterministic Equivalents for Adaptive LDA) shrinkage.

Parameters:

Name	Type	Description	Default
evals	ndarray	Eigenvalues of the empirical covariance matrix.	required
z_vec	ndarray	Vector of interest projected into the eigenvector space.	required
n_eff	int	Effective number of samples used to compute the empirical covariance matrix.	required
gamma_min	float	Minimum value for the gamma bounded search. Default is 0.02.	0.02
gamma_max	float	Maximum value for the gamma bounded search. Default is 100.	100
base_shrinkage	(lw_analytical, empirical)	Shrinkage method for the base eigenvalue estimation. Default is 'lw_analytical'.	'lw_analytical'
surrogate_shrinkage	(lw_analytical, empirical)	Shrinkage method for the surrogate eigenvalue estimation. Default is 'lw_analytical'.	'lw_analytical'
eps	float	Epsilon for numerical stability. Default is 1e-8.	1e-08

Notes

The DEAL method utilizes the Random Matrix Theory (RMT) [1] to construct an estimate and optimize an objective defined as an expectation over the data distribution of \(|| \hat\Sigma^{-1} \mu - \Sigma^{-1} \mu ||_\Sigma^2\) where \(\hat\Sigma\) is an optimal linear correction to any non-directional shrinkage, \(\Sigma\) is the True Population Covariance, \(\mu\) is a constant vector of interest and \(|| \cdot ||_\Sigma\) is the Mahalanobis distance based on the matrix \(\Sigma\). More details in a future paper.

Returns:

Type	Description
ndarray	Shrinkage-adjusted eigenvalues.

References

---

1. Hachem, W., Loubaton, P., Najim, J., & Vallet, P. (2013). On bilinear forms based on the resolvent of large random matrices. In Annales de l'IHP Probabilités et statistiques (Vol. 49, No. 1, pp. 36-63). https://www.numdam.org/article/AIHPB_2013__49_1_36_0.pdf ↩

Source code in shrinkr/functional/_deal.py

def deal(
    evals: np.ndarray,
    z_vec: np.ndarray,
    n_eff: int,
    gamma_min: float = 0.02,
    gamma_max: float = 100,
    base_shrinkage: EigenvalueOptions = "lw_analytical",
    surrogate_shrinkage: EigenvalueOptions = "lw_analytical",
    eps=1e-8,
    **kwargs,
):
    r"""DEAL (Deterministic Equivalents for Adaptive LDA) shrinkage.

    Parameters
    ----------
    evals : np.ndarray
        Eigenvalues of the empirical covariance matrix.
    z_vec : np.ndarray
        Vector of interest projected into the eigenvector space.
    n_eff : int
        Effective number of samples used to compute the empirical covariance matrix.
    gamma_min : float, optional
        Minimum value for the gamma bounded search. Default is 0.02.
    gamma_max : float, optional
        Maximum value for the gamma bounded search. Default is 100.
    base_shrinkage : {'lw_analytical', 'empirical'}, optional
        Shrinkage method for the base eigenvalue estimation. Default is 'lw_analytical'.
    surrogate_shrinkage : {'lw_analytical', 'empirical'}, optional
        Shrinkage method for the surrogate eigenvalue estimation. Default is 'lw_analytical'.
    eps : float, optional
        Epsilon for numerical stability. Default is 1e-8.

    Notes
    -----
    The DEAL method utilizes the Random Matrix Theory (RMT) [1] to
    construct an estimate and optimize an objective
    defined as an expectation over the data distribution of
    $|| \hat\Sigma^{-1} \mu - \Sigma^{-1} \mu ||_\Sigma^2$ where $\hat\Sigma$
    is an optimal linear correction to any non-directional shrinkage,
    $\Sigma$ is the True Population Covariance, $\mu$ is a constant vector
    of interest and $|| \cdot ||_\Sigma$ is the Mahalanobis distance
    based on the matrix $\Sigma$. More details in a future paper.

    Returns
    -------
    np.ndarray
        Shrinkage-adjusted eigenvalues.

    References
    ----------
    [^1]: Hachem, W., Loubaton, P., Najim, J., & Vallet, P. (2013).
        On bilinear forms based on the resolvent of large random matrices.
        In Annales de l'IHP Probabilités et statistiques (Vol. 49, No. 1, pp. 36-63).
        <https://www.numdam.org/article/AIHPB_2013__49_1_36_0.pdf>
    """
    # Rescale eigenvalues to Trace p
    evals /= np.mean(evals)
    p = evals.shape[0]

    # Compute (non-)linear shrinkages
    orig_evals = evals.copy()
    lw_nl_evals = py_lw_analytical(evals, n_eff, p, eps**2)

    # Select which eigenvalues to use for the base
    if base_shrinkage == "lw_analytical":
        evals = lw_nl_evals
    elif base_shrinkage == "empirical":
        evals = orig_evals
    else:
        raise ValueError("Unknown base shrinkage")

    # ... and the surrogate estimation
    if surrogate_shrinkage == "lw_analytical":
        surrogate_evals = lw_nl_evals
    elif surrogate_shrinkage == "empirical":
        surrogate_evals = evals
    else:
        raise ValueError("Unknown surrogate shrinkage")

    shrinkage = py_deal(evals, surrogate_evals, z_vec, gamma_min, gamma_max, n_eff, p)

    return evals + shrinkage

deal_objective(base_evals, surrogate_evals, z_vec, gamma, n, start_value=1.0)

Objective function of DEAL.

Computes the optimization objective using deterministic equivalents. Requires solving a fixed-point equation for delta at a given gamma.

Parameters:

Name	Type	Description	Default
base_evals	ndarray	First 1D array of eigenvalues for the objective (Those will be shrunk)	required
surrogate_evals	ndarray	Second 1D array of eigenvalues for the objective (Used to compute shrinkage paramters)	required
z_vec	ndarray	Vector of interest projected into the eigenvector space.	required
gamma	float	The value of gamma to evaluate. During optimization only this value changes.	required
n	int	Effective number of samples used to compute the empirical covariance matrix.	required
start_value	float	Starting value of delta for the fixed point iteration method used by for the objective.	1.0

See Also

shrinkr.functional.deal function for more information about the DEAL method.

Returns:

Type	Description
float	The DEAL objective estimate.

Source code in shrinkr/functional/_deal.py

def deal_objective(
    base_evals: np.ndarray,
    surrogate_evals: np.ndarray,
    z_vec: np.ndarray,
    gamma: float,
    n: int,
    start_value: float = 1.0,
):
    """Objective function of DEAL.

    Computes the optimization objective using deterministic equivalents.
    Requires solving a fixed-point equation for delta at a given gamma.

    Parameters
    ----------
    base_evals : np.ndarray
        First 1D array of eigenvalues for the objective (Those will be shrunk)
    surrogate_evals : np.ndarray
        Second 1D array of eigenvalues for the objective (Used to compute shrinkage paramters)
    z_vec : np.ndarray
        Vector of interest projected into the eigenvector space.
    gamma : float
        The value of gamma to evaluate. During optimization only this value changes.
    n : int
        Effective number of samples used to compute the empirical covariance matrix.
    start_value : float, optional
        Starting value of delta for the fixed point iteration method used by for the objective.

    See Also
    --------
    [`shrinkr.functional.deal`][]
        function for more information about the DEAL method.

    Returns
    -------
    float
        The DEAL objective estimate.
    """
    p = int(base_evals.shape[0])
    return py_deal_objective(base_evals, surrogate_evals, z_vec, gamma, start_value, n, p)

loss_fm(v, sigma, mu)

Fisher Margin (FM) loss.

Defined as \(FM(v) = -(v^T \mu)^2 / (v^T \Sigma v)\), where \(\mu\) is the true difference-in-means vector, \(\Sigma\) is the true Population Covariance matrix and the \(v\) is the considered LDA vector.

Minimizing the Fisher Margin leads to an optimal Bayesian Classifier on data which admits the LDA data assumptions. The loss is scale invariant.

Parameters:

Name	Type	Description	Default
v	ndarray	LDA vector computed from data.	required
sigma	ndarray	True covariance matrix.	required
mu	ndarray	True difference-in-means vector.	required

Notes

Practically the Fisher Margin is defined without the minus sign. It is there only to turn the maximization task in a minimization one making it a loss.

Returns:

Type	Description
float	Value of the FM loss.

Source code in shrinkr/functional/_losses.py

def loss_fm(v: np.ndarray, sigma: np.ndarray, mu: np.ndarray) -> float:
    r"""Fisher Margin (FM) loss.

    Defined as $FM(v) = -(v^T \mu)^2 / (v^T \Sigma v)$,
    where $\mu$ is the true difference-in-means vector,
    $\Sigma$ is the true Population Covariance matrix
    and the $v$ is the considered LDA vector.

    Minimizing the Fisher Margin leads to an optimal
    Bayesian Classifier on data which admits the LDA
    data assumptions. The loss is scale invariant.

    Parameters
    ----------
    v : np.ndarray
        LDA vector computed from data.
    sigma : np.ndarray
        True covariance matrix.
    mu : np.ndarray
        True difference-in-means vector.

    Notes
    -----
    Practically the Fisher Margin is defined
    without the minus sign. It is there only to turn
    the maximization task in a minimization one
    making it a `loss`.

    Returns
    -------
    float
        Value of the FM loss.
    """
    if len(v.shape) != 1:
        raise ValueError("v has to be a 1D vector")
    if len(mu.shape) != 1:
        raise ValueError("v has to be a 1D vector")
    if len(sigma.shape) != 2:
        raise ValueError("sigma has to be a matrix")

    A = np.dot(v, mu)
    B = v.T @ (sigma @ v)
    return -(A**2) / B

loss_fr(matrixA, matrixB)

Frobenius distance between two matrices.

Parameters:

Name	Type	Description	Default
matrixA	ndarray	First matrix.	required
matrixB	ndarray	Second matrix.	required

Returns:

Type	Description
float	Scaled squared Frobenius distance between the matrices.

Source code in shrinkr/functional/_losses.py

def loss_fr(matrixA: np.ndarray, matrixB: np.ndarray) -> float:
    """Frobenius distance between two matrices.

    Parameters
    ----------
    matrixA : np.ndarray
        First matrix.
    matrixB : np.ndarray
        Second matrix.

    Returns
    -------
    float
        Scaled squared Frobenius distance between the matrices.
    """
    # Checks
    if len(matrixA.shape) != 2:
        raise ValueError("matrixB hat has to be a matrix")
    if len(matrixB.shape) != 2:
        raise ValueError("matrixB has to be a matrix")
    if matrixA.shape != matrixB.shape:
        raise ValueError("matrixB hat has to have the same shape as matrixB")

    # Logic
    n, p = matrixB.shape
    delta = matrixA - matrixB
    return np.sum(delta.reshape(-1) ** 2) / p

loss_mv(sigma_hat, sigma)

Minimal Variance (MV) loss [1].

Parameters:

Name	Type	Description	Default
sigma_hat	ndarray	Estimated covariance matrix.	required
sigma	ndarray	True covariance matrix.	required

Returns:

Type	Description
float	Value of the MV loss.

References

---

1. Ledoit, O., & Wolf, M. (2020). Analytical nonlinear shrinkage of large-dimensional covariance matrices. The Annals of Statistics, 48(5), 3043-3065. http://www.ledoit.net/Analytical_AoS_2020.pdf ↩

Source code in shrinkr/functional/_losses.py

def loss_mv(sigma_hat: np.ndarray, sigma: np.ndarray) -> float:
    """Minimal Variance (MV) loss [1].

    Parameters
    ----------
    sigma_hat : np.ndarray
        Estimated covariance matrix.
    sigma : np.ndarray
        True covariance matrix.

    Returns
    -------
    float
        Value of the MV loss.

    References
    ----------
    [^1]: Ledoit, O., & Wolf, M. (2020).
        Analytical nonlinear shrinkage of large-dimensional covariance matrices.
        The Annals of Statistics, 48(5), 3043-3065.
        <http://www.ledoit.net/Analytical_AoS_2020.pdf>
    """
    # Checks
    if len(sigma_hat.shape) != 2:
        raise ValueError("sigma hat has to be a matrix")
    if len(sigma.shape) != 2:
        raise ValueError("sigma has to be a matrix")
    if sigma_hat.shape != sigma.shape:
        raise ValueError("sigma hat has to have the same shape as matrixB")

    # Logic
    n, p = sigma.shape
    omega_hat = np.linalg.inv(sigma_hat)
    num = np.trace(np.dot(np.dot(omega_hat, sigma), omega_hat)) / p
    denom = (np.trace(omega_hat) / p) ** 2
    alpha = np.trace(np.linalg.inv(sigma)) / p
    return num / denom - alpha

loss_prial(sample_cov, sigma_hat, sigma)

Percentage Relative Improvement in Average Loss (PRIAL) [1].

Parameters:

Name	Type	Description	Default
sample_cov	ndarray	Sample covariance matrix.	required
sigma_hat	ndarray	Estimated covariance matrix.	required
sigma	ndarray	True covariance matrix.	required

Returns:

Type	Description
float	Percentage improvement relative to the oracle, in the range [0, 1].

References

---

1. Ledoit, O., & Péché, S. (2011). Eigenvectors of some large sample covariance matrix ensembles. Probability Theory and Related Fields, 151(1), 233-264. https://link.springer.com/article/10.1007/s00440-010-0298-3 ↩

Source code in shrinkr/functional/_losses.py

def loss_prial(sample_cov: np.ndarray, sigma_hat: np.ndarray, sigma: np.ndarray) -> float:
    """Percentage Relative Improvement in Average Loss (PRIAL) [1].

    Parameters
    ----------
    sample_cov : np.ndarray
        Sample covariance matrix.
    sigma_hat : np.ndarray
        Estimated covariance matrix.
    sigma : np.ndarray
        True covariance matrix.

    Returns
    -------
    float
        Percentage improvement relative to the oracle, in the range [0, 1].

    References
    ----------
    [^1]: Ledoit, O., & Péché, S. (2011).
        Eigenvectors of some large sample covariance matrix ensembles.
        Probability Theory and Related Fields, 151(1), 233-264.
        <https://link.springer.com/article/10.1007/s00440-010-0298-3>
    """
    # Checks
    if len(sample_cov.shape) != 2:
        raise ValueError("Sigma hat has to be a matrix")
    if len(sigma_hat.shape) != 2:
        raise ValueError("Sigma has to be a matrix")
    if len(sigma.shape) != 2:
        raise ValueError("Sigma has to be a matrix")
    if sample_cov.shape != sigma.shape:
        raise ValueError("Sample cov has to have the same shape as Sigma")
    if sigma.shape != sigma_hat.shape:
        raise ValueError("Sigma hat has to have the same shape as Sigma")

    # Logic
    num = loss_mv(sample_cov, sigma) - loss_mv(sigma_hat, sigma)
    sigma_ast = mv_opt_cov(sample_cov, sigma)
    denom = loss_mv(sample_cov, sigma) - loss_mv(sigma_ast, sigma)
    return num / float(denom)

lw_analytical(eigenvalues, n, p=None, eps=1e-08)

Ledoit-Wolf Analytical (nonlinear) shrinkage of eigenvalues.

Based on a optimization free formula from Ledoit and Wolf (2020) [1]. Handles also the high-dimensional setting where \(p>n\).

Parameters:

Name	Type	Description	Default
eigenvalues	ndarray	1-D array of eigenvalues of the sample covariance matrix.	required
n	int	Number of observations used to compute the sample covariance.	required
p	int	Number of variables. If None, inferred as len(eigenvalues).	None
eps	float	Threshold below which eigenvalues are treated as numerically zero. Default is 1e-8.	1e-08

Returns:

Type	Description
ndarray	Analytically shrunk eigenvalues of the same shape as eigenvalues.

References

---

1. Ledoit, O., & Wolf, M. (2020). Analytical nonlinear shrinkage of large-dimensional covariance matrices. The Annals of Statistics, 48(5), 3043-3065. http://www.ledoit.net/Analytical_AoS_2020.pdf ↩

Source code in shrinkr/functional/_lw_analytical.py

def lw_analytical(
    eigenvalues: np.ndarray, n: int, p: int | None = None, eps: float = 1e-8
) -> np.ndarray:
    """Ledoit-Wolf Analytical (nonlinear) shrinkage of eigenvalues.

    Based on a optimization free formula from Ledoit and Wolf (2020) [1].
    Handles also the high-dimensional setting where $p>n$.

    Parameters
    ----------
    eigenvalues : np.ndarray
        1-D array of eigenvalues of the sample covariance matrix.
    n : int
        Number of observations used to compute the sample covariance.
    p : int, optional
        Number of variables. If None, inferred as ``len(eigenvalues)``.
    eps : float, optional
        Threshold below which eigenvalues are treated as numerically zero.
        Default is 1e-8.

    Returns
    -------
    np.ndarray
        Analytically shrunk eigenvalues of the same shape as ``eigenvalues``.

    References
    ----------
    [^1]: Ledoit, O., & Wolf, M. (2020).
        Analytical nonlinear shrinkage of large-dimensional covariance matrices.
        The Annals of Statistics, 48(5), 3043-3065.
        <http://www.ledoit.net/Analytical_AoS_2020.pdf>
    """
    # Checks
    if len(eigenvalues.shape) != 1:
        raise ValueError("eigenvalues have to be a flat ndarray")

    # Type conversion
    eigenvalues: np.ndarray = np.ascontiguousarray(eigenvalues)
    eigenvalues: NDArray[np.float64] = eigenvalues.astype(np.float64)

    if p is None:
        p = len(eigenvalues)

    # n correction
    num_not_effective = int(np.sum(eigenvalues[np.maximum(0, p - n) :] < eps))
    if num_not_effective > 0:
        n = n - num_not_effective  # Correct the number of effective samples

    return py_lw_analytical(eigenvalues, n, p, eps)

lw_linear(X, assume_centered=False)

Ledoit-Wolf linear shrinkage estimator.

The value of the shrinkage is constructed based on the Theorem 3.2 and Lemmata 3.2-3.5 from [1].

Parameters:

Name	Type	Description	Default
X	ndarray	Data matrix of shape (n_samples, n_features).	required
assume_centered	bool	If True, data is not mean-centered before computing the covariance. Default is False.	False

Notes

The regularized covariance is: \((1 - s) * S_c + s * \mu * I\), where \(\mu = Tr{(S_c)} / n_\text{features}\), \(s\) is the shrinkage value and \(S_c\) is the sample covariance matrix.

Returns:

Name	Type	Description
sample_cov_star	ndarray	Shrinkage-regularized covariance matrix of shape (n_features, n_features).
shrinkage	float	Optimal shrinkage coefficient.

References

---

1. Ledoit, O., & Wolf, M. (2004). A well-conditioned estimator for large-dimensional covariance matrices. Journal of multivariate analysis, 88(2), 365-411. http://www.ledoit.net/ole1a.pdf ↩

Source code in shrinkr/functional/_lw_linear.py

def lw_linear(X: np.ndarray, assume_centered: bool = False) -> tuple[np.ndarray, float]:
    r"""Ledoit-Wolf linear shrinkage estimator.

    The value of the shrinkage is constructed
    based on the Theorem 3.2 and Lemmata 3.2-3.5 from [1].

    Parameters
    ----------
    X : np.ndarray
        Data matrix of shape (n_samples, n_features).
    assume_centered : bool, optional
        If True, data is not mean-centered before computing the covariance.
        Default is False.

    Notes
    -----
    The regularized covariance is:
    $(1 - s) * S_c + s * \mu * I$,
    where $\mu = Tr{(S_c)} / n_\text{features}$, $s$ is the shrinkage value
    and $S_c$ is the sample covariance matrix.

    Returns
    -------
    sample_cov_star : np.ndarray
        Shrinkage-regularized covariance matrix of shape (n_features, n_features).
    shrinkage : float
        Optimal shrinkage coefficient.

    References
    ----------
    [^1]: Ledoit, O., & Wolf, M. (2004).
        A well-conditioned estimator for large-dimensional covariance matrices.
        Journal of multivariate analysis, 88(2), 365-411.
        <http://www.ledoit.net/ole1a.pdf>
    """
    # for only one feature, the result is the same whatever the shrinkage
    if len(X.shape) == 2 and X.shape[1] == 1:
        return np.var(X.reshape(-1)).reshape(1, 1)
    if X.ndim == 1:
        X = np.reshape(X, (1, -1))

    if X.shape[0] == 1:
        print("Only one sample available. You may want to reshape your data array")

    n, p = X.shape

    # optionally center data
    if not assume_centered:
        X = X - X.mean(0)

    sample_cov = np.dot(X.T, X) / n
    sample_cov_star, shrinkage = py_lw_linear(X, sample_cov, n, p)

    return sample_cov_star, shrinkage

mv_opt_cov(sample_cov, sigma)

Minimal variance optimal rotation equivariant estimator.

Oracle estimator derived in [1].

Parameters:

Name	Type	Description	Default
sample_cov	ndarray	Sample covariance matrix.	required
sigma	ndarray	True covariance matrix.	required

Returns:

Type	Description
ndarray	Oracle optimal rotation equivariant estimator under the MV loss.

References

---

1. Ledoit, O., & Wolf, M. (2020). Analytical nonlinear shrinkage of large-dimensional covariance matrices. The Annals of Statistics, 48(5), 3043-3065. http://www.ledoit.net/Analytical_AoS_2020.pdf ↩

Source code in shrinkr/functional/_losses.py

def mv_opt_cov(sample_cov: np.ndarray, sigma: np.ndarray) -> np.ndarray:
    """Minimal variance optimal rotation equivariant estimator.

    Oracle estimator derived in [1].

    Parameters
    ----------
    sample_cov : np.ndarray
        Sample covariance matrix.
    sigma : np.ndarray
        True covariance matrix.

    Returns
    -------
    np.ndarray
        Oracle optimal rotation equivariant estimator under the MV loss.

    References
    ----------
    [^1]: Ledoit, O., & Wolf, M. (2020).
        Analytical nonlinear shrinkage of large-dimensional covariance matrices.
        The Annals of Statistics, 48(5), 3043-3065.
        <http://www.ledoit.net/Analytical_AoS_2020.pdf>
    """
    # Checks
    if len(sample_cov.shape) != 2:
        raise ValueError("Sigma hat has to be a matrix")
    if len(sigma.shape) != 2:
        raise ValueError("Sigma has to be a matrix")
    if sample_cov.shape != sigma.shape:
        raise ValueError("Sigma hat has to have the same shape as Sigma")

    # Logic
    lam, u = np.linalg.eigh(sample_cov)
    d_start: np.ndarray = np.einsum("ji, jk, ki -> i", u, sigma, u)
    ud = np.dot(u, np.diag(d_start))
    return np.dot(ud, u.T)

oas(sample_cov, n, p=None)

Oracle Approximating Shrinkage (OAS) covariance estimator.

The formulation is based on [1].

Parameters:

Name	Type	Description	Default
sample_cov	ndarray	Sample covariance matrix of shape (p, p).	required
n	int	Number of observations used to compute the sample covariance.	required
p	int	Number of variables. If None, inferred from sample_cov.	None

Returns:

Name	Type	Description
sample_cov_star	ndarray	Shrinkage-regularized covariance matrix of shape (p, p).
shrinkage	float	Optimal shrinkage coefficient.

References

---

1. Chen, Y., Wiesel, A., Eldar, Y. C., & Hero, A. O. (2010). Shrinkage algorithms for MMSE covariance estimation. IEEE Transactions on Signal Processing, 58(10), 5016-5029. https://arxiv.org/pdf/0907.4698.pdf ↩

Source code in shrinkr/functional/_oas.py

def oas(sample_cov: np.ndarray, n: int, p: int | None = None) -> tuple[np.ndarray, float]:
    """Oracle Approximating Shrinkage (OAS) covariance estimator.

    The formulation is based on [1].

    Parameters
    ----------
    sample_cov : np.ndarray
        Sample covariance matrix of shape (p, p).
    n : int
        Number of observations used to compute the sample covariance.
    p : int, optional
        Number of variables. If None, inferred from ``sample_cov``.

    Returns
    -------
    sample_cov_star : np.ndarray
        Shrinkage-regularized covariance matrix of shape (p, p).
    shrinkage : float
        Optimal shrinkage coefficient.

    References
    ----------
    [^1]: Chen, Y., Wiesel, A., Eldar, Y. C., & Hero, A. O. (2010).
        Shrinkage algorithms for MMSE covariance estimation.
        IEEE Transactions on Signal Processing, 58(10), 5016-5029.
        <https://arxiv.org/pdf/0907.4698.pdf>
    """
    if len(sample_cov.shape) != 2:
        raise ValueError("Expected a square numpy matrix")
    if sample_cov.shape[0] != sample_cov.shape[1]:
        raise ValueError("Expected a square numpy matrix")

    if p is None:
        p = sample_cov.shape[0]

    sample_cov: np.ndarray = sample_cov.astype(np.float64)
    sample_cov = np.ascontiguousarray(sample_cov)

    sample_cov_star, shrinakge = py_oas(sample_cov, n, p)
    return sample_cov_star, shrinakge

---

Reference implementations

shrinkr.reference

Shrinkage reference implementations in numpy.

ref_deal(evals, z_vec, n_eff, gamma_min=0.02, gamma_max=100, base_shrinkage='lw_analytical', surrogate_shrinkage='lw_analytical', eps=1e-08, **kwargs)

DEAL (Deterministic Equivalents for Adaptive LDA) shrinkage (reference implementation).

Parameters:

Name	Type	Description	Default
evals	ndarray	Eigenvalues of the empirical covariance matrix.	required
z_vec	ndarray	Vector of interest projected into the eigenvector space.	required
n_eff	int	Effective number of samples used to compute the empirical covariance matrix.	required
gamma_min	float	Minimum value for the gamma bounded search. Default is 0.02.	0.02
gamma_max	float	Maximum value for the gamma bounded search. Default is 100.	100
base_shrinkage	(lw_analytical, empirical)	Shrinkage method for the base eigenvalue estimation. Default is 'lw_analytical'.	'lw_analytical'
surrogate_shrinkage	(lw_analytical, empirical)	Shrinkage method for the surrogate eigenvalue estimation. Default is 'lw_analytical'.	'lw_analytical'
eps	float	Epsilon for numerical stability. Default is 1e-8.	1e-08

Returns:

Type	Description
ndarray	Shrinkage-adjusted eigenvalues.

See Also

shrinkr.functional.deal Optimized implementation of this method. Go there for additional notes and references. Functions ref_* are reference implementations intended for validation.

Source code in shrinkr/reference/_deal.py

def ref_deal(
    evals: np.ndarray,
    z_vec: np.ndarray,
    n_eff: int,
    gamma_min: float = 0.02,
    gamma_max: float = 100,
    base_shrinkage: EigenvalueOptions = "lw_analytical",
    surrogate_shrinkage: EigenvalueOptions = "lw_analytical",
    eps=1e-8,
    **kwargs,
):
    """DEAL (Deterministic Equivalents for Adaptive LDA) shrinkage (reference implementation).

    Parameters
    ----------
    evals : np.ndarray
        Eigenvalues of the empirical covariance matrix.
    z_vec : np.ndarray
        Vector of interest projected into the eigenvector space.
    n_eff : int
        Effective number of samples used to compute the empirical covariance matrix.
    gamma_min : float, optional
        Minimum value for the gamma bounded search. Default is 0.02.
    gamma_max : float, optional
        Maximum value for the gamma bounded search. Default is 100.
    base_shrinkage : {'lw_analytical', 'empirical'}, optional
        Shrinkage method for the base eigenvalue estimation. Default is 'lw_analytical'.
    surrogate_shrinkage : {'lw_analytical', 'empirical'}, optional
        Shrinkage method for the surrogate eigenvalue estimation. Default is 'lw_analytical'.
    eps : float, optional
        Epsilon for numerical stability. Default is 1e-8.

    Returns
    -------
    np.ndarray
        Shrinkage-adjusted eigenvalues.

    See Also
    --------
    [`shrinkr.functional.deal`][]
        Optimized implementation of this method.
        Go there for additional notes and references.
        Functions ref_* are reference implementations intended for validation.

    """
    # Rescale eigenvalues to Trace p
    evals /= np.mean(evals)

    # Compute (non-)linear shrinkages
    orig_evals = evals.copy()
    lw_nl_evals = ref_lw_analytical(evals, n_eff, eps**2)

    # Select which eigenvalues to use for the base
    if base_shrinkage == "lw_analytical":
        evals = lw_nl_evals
    elif base_shrinkage == "empirical":
        evals = orig_evals
    else:
        raise ValueError("Unknown base shrinkage")

    # ... and the surrogate estimation
    if surrogate_shrinkage == "lw_analytical":
        surrogate_evals = lw_nl_evals
    elif surrogate_shrinkage == "empirical":
        surrogate_evals = evals
    else:
        raise ValueError("Unknown surrogate shrinkage")

    # Optimization
    def objective_fn(gamma, start_delta):
        return ref_deal_objective(
            surrogate_evals=surrogate_evals,
            base_evals=evals,
            z_vec=z_vec,
            n=n_eff,
            gamma=gamma,
            start_value=start_delta,
            eps=eps,
        )

    optimal_gamma, _, ratio = golden_section_search(
        objective_fn, gamma_min, gamma_max, initial_delta=1, tol=eps
    )

    return evals + optimal_gamma

ref_deal_objective(base_evals, surrogate_evals, z_vec, gamma, n, start_value=1, max_iters=200, eps=1e-08)

Objective function of DEAL (reference implementation).

Computes the optimization objective using deterministic equivalents. Requires solving a fixed-point equation for delta at a given gamma.

Parameters:

Name	Type	Description	Default
base_evals	ndarray	Eigenvalue estimates for the base covariance matrix (those that will be shrunk).	required
surrogate_evals	ndarray	Eigenvalue estimates used to compute the shrinkage parameters.	required
z_vec	ndarray	Vector of interest projected into the eigenvector space.	required
gamma	float	Scalar resolvent parameter for the matrix S.	required
n	int	Effective number of observations used to compute the sample covariance matrix.	required
start_value	float	Starting value for the fixed-point iteration (supports warm-starts). Default is 1.	1
max_iters	int	Maximum number of fixed-point iterations. Default is 200.	200
eps	float	Convergence tolerance for early stopping. Default is 1e-8.	1e-08

Returns:

Name	Type	Description
obj	float	The DEAL risk objective estimate.
delta	float	The converged delta value from the fixed-point iteration.
delta_prime	float	The derivative of delta with respect to gamma.

See Also

shrinkr.functional.deal_objective Optimized implementation of this method. Go there for additional notes and references. Functions ref_* are reference implementations intended for validation.

Source code in shrinkr/reference/_deal.py

def ref_deal_objective(
    base_evals: np.ndarray,
    surrogate_evals: np.ndarray,
    z_vec: np.ndarray,
    gamma: float,
    n: int,
    start_value: float = 1,
    max_iters: int = 200,
    eps: float = 1e-8,
):
    """Objective function of DEAL (reference implementation).

    Computes the optimization objective using deterministic equivalents.
    Requires solving a fixed-point equation for delta at a given gamma.

    Parameters
    ----------
    base_evals : np.ndarray
        Eigenvalue estimates for the base covariance matrix (those that will be shrunk).
    surrogate_evals : np.ndarray
        Eigenvalue estimates used to compute the shrinkage parameters.
    z_vec : np.ndarray
        Vector of interest projected into the eigenvector space.
    gamma : float
        Scalar resolvent parameter for the matrix S.
    n : int
        Effective number of observations used to compute the sample covariance matrix.
    start_value : float, optional
        Starting value for the fixed-point iteration (supports warm-starts). Default is 1.
    max_iters : int, optional
        Maximum number of fixed-point iterations. Default is 200.
    eps : float, optional
        Convergence tolerance for early stopping. Default is 1e-8.

    Returns
    -------
    obj : float
        The DEAL risk objective estimate.
    delta : float
        The converged delta value from the fixed-point iteration.
    delta_prime : float
        The derivative of delta with respect to gamma.

    See Also
    --------
    [`shrinkr.functional.deal_objective`][]
        Optimized implementation of this method.
        Go there for additional notes and references.
        Functions ref_* are reference implementations intended for validation.

    """
    # Compute delta (and beta) and diagT via fixed point iterations
    # diagT is the vector of the diagonal of the T (deterministic_equivalent) matrix
    # invDiagT is the vector of the inverse diagonal of the T matrix
    delta = start_value
    for _ in range(max_iters):
        beta = 1 / (1 + delta)
        invDiagT = gamma + surrogate_evals * beta
        new_delta = np.sum(surrogate_evals / invDiagT) / n
        if abs(new_delta - delta) < eps:
            break
        delta = new_delta
    else:
        print("WARN: Fixed point method did not coverge.")

    # Compute derivate of delta
    diagT = 1 / invDiagT
    a = np.sum((surrogate_evals) * (diagT**2))
    b = (beta**2) * np.sum((surrogate_evals * diagT) ** 2)
    delta_prime = (-a) / (n - b)

    # Compute intermediate values
    diagT_sq = diagT**2  # T^2
    z_vec_sq = z_vec**2  # z^2_i
    delta_inv_diagT = 1 - delta_prime * surrogate_evals * (
        beta**2
    )  # derivaite of the inverse of T wrt gamma

    # Compute both objective parts
    # ApN_gamma = -np.sum((z_vec_sq ) * diagT_sq * delta_inv_diagT)
    # B_gamma = np.sum(diagT * z_vec_sq / evals)
    # B_gamma_sq = B_gamma ** 2

    # Different formula for Mahalonobis based distance
    z_vec_sq_e = (z_vec_sq) * base_evals
    B_gamma = np.sum(z_vec_sq * diagT)
    B_gamma_sq = B_gamma**2
    ApN_gamma = -np.sum(z_vec_sq_e * diagT_sq * delta_inv_diagT)

    obj = B_gamma_sq / ApN_gamma

    return obj, delta, delta_prime

ref_lw_analytical(lam, n, eps=1e-08)

Ledoit-Wolf Analytical (nonlinear) shrinkage of eigenvalues (reference implementation).

Parameters:

Name	Type	Description	Default
lam	ndarray	1-D array of empirical eigenvalues.	required
n	int	Effective sample size.	required
eps	float	Threshold below which eigenvalues are treated as numerically zero. Default is 1e-8.	1e-08

Returns:

Type	Description
ndarray	Analytically shrunk eigenvalues.

See Also

shrinkr.functional.lw_analytical Optimized implementation of this method. This is a reference implementation intended for validation; prefer the functional version for performance-critical use.

Source code in shrinkr/reference/_lw_analytical.py

def ref_lw_analytical(lam: np.ndarray, n: int, eps: float = 1e-8) -> np.ndarray:
    """Ledoit-Wolf Analytical (nonlinear) shrinkage of eigenvalues (reference implementation).

    Parameters
    ----------
    lam : np.ndarray
        1-D array of empirical eigenvalues.
    n : int
        Effective sample size.
    eps : float, optional
        Threshold below which eigenvalues are treated as numerically zero.
        Default is 1e-8.

    Returns
    -------
    np.ndarray
        Analytically shrunk eigenvalues.

    See Also
    --------
    [`shrinkr.functional.lw_analytical`][]
        Optimized implementation of this method. This is a reference
        implementation intended for validation; prefer the functional version
        for performance-critical use.

    """
    assert len(lam.shape) == 1

    lam = lam.astype(np.float64)
    p = lam.shape[0]
    sqrt5 = np.sqrt(5)

    num_not_effective = int(np.sum(lam[np.maximum(0, p - n) :] < eps))
    if num_not_effective > 0:
        n = n - num_not_effective  # Correct the number of effective samples

    lam = lam[np.maximum(0, p - n) :]  # Use only the effective eigenvalues

    L = np.tile(lam[:, None], (1, np.minimum(p, n)))
    h = np.power(n, -1 / 3.0)
    # % Equation(4.9)

    H = h * L.T
    x = (L - L.T) / H
    abs_x = np.abs(x)

    ftilde = (3 / 4.0 / sqrt5) * np.mean(np.maximum(1 - (x**2) / 5.0, 0) / H, 1)
    # % Equation(4.7)

    Hftemp = np.zeros_like(x)
    # --- Evaluate Regime 1: Singularities ---
    # The log term multiplier (1 - x^2/5) becomes exactly 0, leaving only the linear term.
    singular_mask = np.isclose(abs_x, sqrt5, atol=eps)
    Hftemp[singular_mask] = (-3 / 10 / np.pi) * x[singular_mask]
    # --- Evaluate Regime 2: Large Tails ---
    # Use the asymptotic expansion to avoid subtracting massive numbers.
    # H ~= -1/(pi*x) * [1 + 1/x^2 + 15/(7x^4) + ...]
    large_mask = abs_x > 5.0
    xl = x[large_mask]
    xl2 = xl**2
    Hftemp[large_mask] = (-1 / (np.pi * xl)) * (1 + 1 / xl2 + 15 / (7 * xl2**2))
    # --- Evaluate Regime 3: Normal Domain ---
    normal_mask = ~(singular_mask | large_mask)
    xn = x[normal_mask]
    log_term = np.log(np.abs((sqrt5 - xn) / (sqrt5 + xn)))
    linear_term = (-3 / 10 / np.pi) * xn
    Hftemp[normal_mask] = linear_term + (3 / 4 / sqrt5 / np.pi) * (1 - (xn**2) / 5.0) * log_term

    Hftilde = np.mean(Hftemp / H, 1)
    # % Equation(4.8)

    if p <= n:
        denom1 = (np.pi * (p / n) * lam * ftilde) ** 2
        denom2 = (1 - (p / n) - np.pi * (p / n) * lam * Hftilde) ** 2
        dtilde = lam / (denom1 + denom2)
        # % Equation(4.3)
    else:
        Hftilde0 = (
            (1 / np.pi)
            * (
                3 / 10.0 / h**2
                + 3
                / 4.0
                / sqrt5
                / h
                * (1 - 1 / 5.0 / h**2)
                * np.log((1 + sqrt5 * h) / (1 - sqrt5 * h))
            )
            * np.mean(1 / lam)
        )
        # % Equation(C.8)
        dtilde0 = 1 / (np.pi * (p - n) / n * Hftilde0)
        # % Equation(C.5)
        dtilde1 = lam / (np.pi**2 * lam**2.0 * (ftilde**2 + Hftilde**2))
        # % Eq. (C.4)
        dtilde = np.concatenate([dtilde0 * np.ones(p - n), dtilde1])

    return dtilde

ref_lw_analytical_unstable(lam, n, eps=1e-08)

Ledoit-Wolf Analytical (nonlinear) shrinkage — numerically unstable variant.

Identical in formula to ref_lw_analytical but without the numerical stability fixes for singularities and large tails. Kept for reference only; use ref_lw_analytical for proper reference.

Reference implementation from https://github.com/matzhaugen/analytic_shrinkage.

Parameters:

Name	Type	Description	Default
lam	ndarray	1-D array of empirical eigenvalues.	required
n	int	Effective sample size.	required
eps	float	Threshold below which eigenvalues are treated as numerically zero. Default is 1e-8.	1e-08

Returns:

Type	Description
ndarray	Analytically shrunk eigenvalues.

See Also

shrinkr.functional.lw_analytical Optimized, numerically stable implementation. This is a reference implementation intended for validation.

Source code in shrinkr/reference/_lw_analytical.py

def ref_lw_analytical_unstable(lam: np.ndarray, n: int, eps: float = 1e-8):
    """Ledoit-Wolf Analytical (nonlinear) shrinkage — numerically unstable variant.

    Identical in formula to [`ref_lw_analytical`][shrinkr.reference.ref_lw_analytical] but without the numerical
    stability fixes for singularities and large tails. Kept for reference only;
    use [`ref_lw_analytical`][shrinkr.reference.ref_lw_analytical] for proper reference.

    Reference implementation from <https://github.com/matzhaugen/analytic_shrinkage>.

    Parameters
    ----------
    lam : np.ndarray
        1-D array of empirical eigenvalues.
    n : int
        Effective sample size.
    eps : float, optional
        Threshold below which eigenvalues are treated as numerically zero.
        Default is 1e-8.

    Returns
    -------
    np.ndarray
        Analytically shrunk eigenvalues.

    See Also
    --------
    [`shrinkr.functional.lw_analytical`][]
        Optimized, numerically stable implementation. This is a reference
        implementation intended for validation.

    """
    assert len(lam.shape) == 1

    lam = lam.astype(np.float64)
    p = lam.shape[0]

    # compute analytical nonlinear shrinkage kernel formula
    lam = lam[np.maximum(0, p - n) :]
    if any(lam / sum(lam) < eps):
        raise ValueError("Matrix is singular")

    L = np.tile(lam[:, None], (1, np.minimum(p, n)))
    h = np.power(n, -1 / 3.0)
    # % Equation(4.9)
    H = h * L.T
    x = (L - L.T) / H
    ftilde = (3 / 4.0 / np.sqrt(5)) * np.mean(np.maximum(1 - x**2.0 / 5.0, 0) / H, 1)
    # % Equation(4.7)
    Hftemp = (-3 / 10 / np.pi) * x + (3 / 4.0 / np.sqrt(5) / np.pi) * (1 - x**2.0 / 5.0) * np.log(
        np.abs((np.sqrt(5) - x) / (np.sqrt(5) + x))
    )
    # % Equation(4.8)
    Hftemp[np.abs(x) == np.sqrt(5)] = (-3 / 10 / np.pi) * x[np.abs(x) == np.sqrt(5)]
    Hftilde = np.mean(Hftemp / H, 1)
    if p <= n:
        dtilde = lam / (
            (np.pi * (p / n) * lam * ftilde) ** 2
            + (1 - (p / n) - np.pi * (p / n) * lam * Hftilde) ** 2
        )
    # % Equation(4.3)
    else:
        Hftilde0 = (
            (1 / np.pi)
            * (
                3 / 10.0 / h**2
                + 3
                / 4.0
                / np.sqrt(5)
                / h
                * (1 - 1 / 5.0 / h**2)
                * np.log((1 + np.sqrt(5) * h) / (1 - np.sqrt(5) * h))
            )
            * np.mean(1 / lam)
        )
        # % Equation(C.8)
        dtilde0 = 1 / (np.pi * (p - n) / n * Hftilde0)
        # % Equation(C.5)
        dtilde1 = lam / (np.pi**2 * lam**2.0 * (ftilde**2 + Hftilde**2))
        # % Eq. (C.4)
        dtilde = np.concatenate([dtilde0 * np.ones(p - n), dtilde1])

    return dtilde

ref_lw_linear(X, assume_centered=False, block_size=1000)

Ledoit-Wolf linear shrinkage estimator (reference implementation).

Reference implementation from scikit-learn.

Parameters:

Name	Type	Description	Default
X	ndarray	Data matrix of shape (n_samples, n_features).	required
assume_centered	bool	If True, data is not mean-centered before computing the covariance. Default is False.	False
block_size	int	Size of blocks into which the covariance matrix is split for computation. Default is 1000.	1000

Returns:

Name	Type	Description
sample_cov_star	ndarray	Shrinkage-regularized covariance matrix of shape (n_features, n_features).
shrinkage	float	Optimal shrinkage coefficient.

See Also

shrinkr.functional.lw_linear Optimized implementation of this method. Go there for additional notes and references. Functions ref_* are reference implementations intended for validation.

Source code in shrinkr/reference/_lw_linear.py

def ref_lw_linear(
    X: np.ndarray, assume_centered: bool = False, block_size: int = 1000
) -> tuple[np.ndarray, float]:
    """Ledoit-Wolf linear shrinkage estimator (reference implementation).

    Reference implementation from [scikit-learn](https://github.com/scikit-learn/scikit-learn/blob/fe2edb3cd/sklearn/covariance/_shrunk_covariance.py).

    Parameters
    ----------
    X : np.ndarray
        Data matrix of shape (n_samples, n_features).
    assume_centered : bool, optional
        If True, data is not mean-centered before computing the covariance.
        Default is False.
    block_size : int, optional
        Size of blocks into which the covariance matrix is split for computation.
        Default is 1000.

    Returns
    -------
    sample_cov_star : np.ndarray
        Shrinkage-regularized covariance matrix of shape (n_features, n_features).
    shrinkage : float
        Optimal shrinkage coefficient.

    See Also
    --------
    [`shrinkr.functional.lw_linear`][]
        Optimized implementation of this method.
        Go there for additional notes and references.
        Functions ref_* are reference implementations intended for validation.

    """
    # for only one feature, the result is the same whatever the shrinkage
    if len(X.shape) == 2 and X.shape[1] == 1:
        return np.var(X.reshape(-1)).reshape(1, 1), 0.0
    if X.ndim == 1:
        X = np.reshape(X, (1, -1))

    if X.shape[0] == 1:
        print("Only one sample available. You may want to reshape your data array")
    n_samples, n_features = X.shape

    # optionally center data
    if not assume_centered:
        X = X - X.mean(0)

    sample_cov = (X.T @ X) / X.shape[0]
    # A non-blocked version of the computation is present in the tests
    # in tests/test_covariance.py

    # number of blocks to split the covariance matrix into
    n_splits = int(n_features / block_size)
    X2 = X**2
    emp_cov_trace = np.sum(X2, axis=0) / n_samples
    mu = np.sum(emp_cov_trace) / n_features

    beta_ = 0.0  # sum of the coefficients of <X2.T, X2>
    delta_ = 0.0  # sum of the *squared* coefficients of <X.T, X>
    # starting block computation
    for i in range(n_splits):
        for j in range(n_splits):
            rows = slice(block_size * i, block_size * (i + 1))
            cols = slice(block_size * j, block_size * (j + 1))
            beta_ += np.sum(np.dot(X2.T[rows], X2[:, cols]))
            delta_ += np.sum(np.dot(X.T[rows], X[:, cols]) ** 2)
        rows = slice(block_size * i, block_size * (i + 1))
        beta_ += np.sum(np.dot(X2.T[rows], X2[:, block_size * n_splits :]))
        delta_ += np.sum(np.dot(X.T[rows], X[:, block_size * n_splits :]) ** 2)

    for j in range(n_splits):
        cols = slice(block_size * j, block_size * (j + 1))
        beta_ += np.sum(np.dot(X2.T[block_size * n_splits :], X2[:, cols]))
        delta_ += np.sum(np.dot(X.T[block_size * n_splits :], X[:, cols]) ** 2)

    delta_ += np.sum(np.dot(X.T[block_size * n_splits :], X[:, block_size * n_splits :]) ** 2)
    delta_ /= n_samples**2
    beta_ += np.sum(np.dot(X2.T[block_size * n_splits :], X2[:, block_size * n_splits :]))

    # use delta_ to compute beta
    beta = 1.0 / (n_features * n_samples) * (beta_ / n_samples - delta_)
    # delta is the sum of the squared coefficients of (<X.T,X> - mu*Id) / p
    delta = delta_ - 2.0 * mu * emp_cov_trace.sum() + n_features * mu**2
    delta /= n_features
    # get final beta as the min between beta and delta
    # We do this to prevent shrinking more than "1", which would invert
    # the value of covariances
    beta = min(beta, delta)
    # finally get shrinkage
    shrinkage = 0 if beta == 0 else beta / delta

    sample_cov_star = ((1 - shrinkage) * sample_cov) + (shrinkage * mu * np.eye(n_features))

    return sample_cov_star, shrinkage

ref_oas(sample_cov, n, p=None)

Oracle Approximating Shrinkage (OAS) covariance estimator (reference implementation).

Reference implementation from scikit-learn.

Parameters:

Name	Type	Description	Default
sample_cov	ndarray	Sample covariance matrix of shape (p, p).	required
n	int	Number of observations used to compute the sample covariance.	required
p	int	Number of variables. If None, inferred from sample_cov.	None

Returns:

Name	Type	Description
sample_cov_star	ndarray	Shrinkage-regularized covariance matrix of shape (p, p).
shrinkage	float	Optimal shrinkage coefficient.

See Also

shrinkr.functional.oas Optimized implementation of this method. Go there for additional notes and references. Functions ref_* are reference implementations intended for validation.

Source code in shrinkr/reference/_oas.py

def ref_oas(sample_cov: np.ndarray, n: int, p: int | None = None) -> tuple[np.ndarray, float]:
    """Oracle Approximating Shrinkage (OAS) covariance estimator (reference implementation).

    Reference implementation from [scikit-learn](https://github.com/scikit-learn/scikit-learn/blob/fe2edb3cd/sklearn/covariance/_shrunk_covariance.py).

    Parameters
    ----------
    sample_cov : np.ndarray
        Sample covariance matrix of shape (p, p).
    n : int
        Number of observations used to compute the sample covariance.
    p : int, optional
        Number of variables. If None, inferred from ``sample_cov``.

    Returns
    -------
    sample_cov_star : np.ndarray
        Shrinkage-regularized covariance matrix of shape (p, p).
    shrinkage : float
        Optimal shrinkage coefficient.

    See Also
    --------
    [`shrinkr.functional.oas`][]
        Optimized implementation of this method.
        Go there for additional notes and references.
        Functions ref_* are reference implementations intended for validation.

    """
    if len(sample_cov.shape) != 2:
        raise ValueError("Expected a square numpy matrix")
    if sample_cov.shape[0] != sample_cov.shape[1]:
        raise ValueError("Expected a square numpy matrix")

    if p is None:
        p = sample_cov.shape[0]

    alpha = np.mean(sample_cov**2)
    mu = np.trace(sample_cov) / p
    mu_squared = mu**2

    # The factor 1 / p**2 will cancel out since it is in both the numerator and
    # denominator
    num = alpha + mu_squared
    den = (n + 1) * (alpha - mu_squared / p)
    shrinkage = 1.0 if den == 0 else min(num / den, 1.0)

    # The shrunk covariance is defined as:
    # (1 - shrinkage) * S + shrinkage * F (cf. Eq. 4 in [1])
    # where S is the empirical covariance and F is the shrinkage target defined as
    # F = trace(S) / n_features * np.identity(n_features) (cf. Eq. 3 in [1])
    shrunk_cov = (1.0 - shrinkage) * sample_cov
    shrunk_cov.flat[:: (p + 1)] += shrinkage * mu

    return shrunk_cov, shrinkage

Monte carlo module

shrinkr.monte_carlo

Submodule with simple monte carlo methods for getting testing data.

get_guassian_lda_samples(p=20, n_per_class=100, seed=0)

Generate data for Gaussian LDA with two classes. Balanced dataset.

Parameters:

Name	Type	Description	Default
p	int	Number of features (dimension). Default is 20.	20
n_per_class	int	Number of samples per class. Default is 100.	100
seed	int	Seed for numpy.random.default_rng. Default is 0.	0

Source code in shrinkr/monte_carlo.py

def get_guassian_lda_samples(p: int = 20, n_per_class: int = 100, seed: int = 0):
    """Generate data for Gaussian LDA with two classes. Balanced dataset.

    Parameters
    ----------
    p : int, optional
        Number of features (dimension). Default is 20.
    n_per_class : int, optional
        Number of samples per class. Default is 100.
    seed : int, optional
        Seed for ``numpy.random.default_rng``. Default is 0.
    """
    _, _, cov = get_large_sample_cov(p, n_per_class, seed)

    rng = np.random.default_rng(seed)

    mu = rng.standard_normal(p)  # The difference in means vector
    mu: np.ndarray = mu / np.linalg.vector_norm(mu)
    mu0 = -mu / 2.0
    mu1 = mu / 2.0

    y = np.concatenate([np.zeros(n_per_class), np.ones(n_per_class)])
    X0 = rng.multivariate_normal(mean=mu0, cov=cov, size=n_per_class)
    X1 = rng.multivariate_normal(mean=mu1, cov=cov, size=n_per_class)
    X = np.concatenate([X0, X1], axis=0)

    return X, y

get_large_sample_cov(p=20, n=200, seed=0, add_diagonal=0.1)

Generate Gaussian data with a random positive semi-definite covariance matrix.

The population covariance is constructed as A @ A.T (normalized to unit trace) plus a small ridge add_linear for numerical stability.

Parameters:

Name	Type	Description	Default
p	int	Number of features (dimension). Default is 20.	20
n	int	Number of samples. Default is 200.	200
seed	int	Seed for numpy.random.default_rng. Default is 0.	0
add_diagonal	float	Value added to the diagonal of the true covariance matrix. Default is 1e-1.	0.1

Returns:

Name	Type	Description
X	np.ndarray of shape (n, p)	Simulated data matrix.
sample_cov	np.ndarray of shape (p, p)	Sample covariance matrix computed from X.
real_cov	np.ndarray of shape (p, p)	Population (true) covariance matrix.

Source code in shrinkr/monte_carlo.py

def get_large_sample_cov(
    p: int = 20, n: int = 200, seed: int = 0, add_diagonal: float = 1e-1
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Generate Gaussian data with a random positive semi-definite covariance matrix.

    The population covariance is constructed as ``A @ A.T`` (normalized to unit
    trace) plus a small ridge ``add_linear`` for numerical stability.

    Parameters
    ----------
    p : int, optional
        Number of features (dimension). Default is 20.
    n : int, optional
        Number of samples. Default is 200.
    seed : int, optional
        Seed for ``numpy.random.default_rng``. Default is 0.
    add_diagonal : float, optional
        Value added to the diagonal of the true covariance matrix. Default is 1e-1.

    Returns
    -------
    X : np.ndarray of shape (n, p)
        Simulated data matrix.
    sample_cov : np.ndarray of shape (p, p)
        Sample covariance matrix computed from ``X``.
    real_cov : np.ndarray of shape (p, p)
        Population (true) covariance matrix.
    """
    rng = np.random.default_rng(seed)

    # Create real cov
    A = rng.standard_normal((p, p))
    real_cov = A @ A.T
    real_cov /= np.trace(real_cov)

    real_cov.flat[:: (p + 1)] += add_diagonal

    # Empirical cov
    X = rng.multivariate_normal(mean=np.zeros(p), cov=real_cov, size=n)
    sample_cov = np.cov(X, rowvar=False)

    return X, sample_cov, real_cov

get_small_sample_cov(n=50, seed=0)

Generate a small sample from a fixed 2-D multivariate normal distribution.

The population covariance is [[0.4, 0.2], [0.2, 0.8]].

Parameters:

Name	Type	Description	Default
n	int	Number of samples to generate. Default is 50.	50
seed	int	Seed for numpy.random.RandomState. Default is 0.	0

Returns:

Name	Type	Description
X	np.ndarray of shape (n, 2)	Simulated data matrix.
sample_cov	np.ndarray of shape (2, 2)	Sample covariance matrix computed from X.
real_cov	np.ndarray of shape (2, 2)	Population (true) covariance matrix.

Source code in shrinkr/monte_carlo.py

def get_small_sample_cov(n: int = 50, seed: int = 0) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Generate a small sample from a fixed 2-D multivariate normal distribution.

    The population covariance is ``[[0.4, 0.2], [0.2, 0.8]]``.

    Parameters
    ----------
    n : int, optional
        Number of samples to generate. Default is 50.
    seed : int, optional
        Seed for ``numpy.random.RandomState``. Default is 0.

    Returns
    -------
    X : np.ndarray of shape (n, 2)
        Simulated data matrix.
    sample_cov : np.ndarray of shape (2, 2)
        Sample covariance matrix computed from ``X``.
    real_cov : np.ndarray of shape (2, 2)
        Population (true) covariance matrix.
    """
    real_cov: np.ndarray = np.array([[0.4, 0.2], [0.2, 0.8]])
    rng = np.random.RandomState(seed)
    X = rng.multivariate_normal(mean=[0, 0], cov=real_cov, size=n)
    sample_cov = np.cov(X, rowvar=False)
    return X, sample_cov, real_cov