sts_cov_estimators.cpp

Go to the documentation of this file.

//=============================================================================================================
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "sts_cov_estimators.h"
 
//=============================================================================================================
// STL INCLUDES
//=============================================================================================================
 
#include <algorithm>
 
//=============================================================================================================
// USED NAMESPACES
//=============================================================================================================
 
using namespace STSLIB;
using namespace Eigen;
 
//=============================================================================================================
// DEFINE MEMBER METHODS
//=============================================================================================================
 
std::pair<MatrixXd, double> StsCovEstimators::ledoitWolf(const MatrixXd& matData)
{
    const int p = static_cast<int>(matData.rows());
    const int n = static_cast<int>(matData.cols());
 
    // Sample covariance (using 1/n, matching Ledoit-Wolf and scikit-learn)
    const MatrixXd S = (matData * matData.transpose()) / static_cast<double>(n);
 
    // Shrinkage target: mu * I_p
    const double mu = S.trace() / static_cast<double>(p);
 
    // delta = ||S - mu * I||^2_F / p
    MatrixXd centered = S - mu * MatrixXd::Identity(p, p);
    const double delta = centered.squaredNorm() / static_cast<double>(p);
 
    // If delta is essentially zero, S is already a scaled identity
    if (delta < 1e-30) {
        return {S, 0.0};
    }
 
    // Efficient computation of beta:
    // sum_k ||x_k x_k^T - S||^2_F = sum_k (x_k^T x_k)^2 - n * ||S||^2_F
    double sumSqNorms = 0.0;
    for (int k = 0; k < n; ++k) {
        const double nrm = matData.col(k).squaredNorm();
        sumSqNorms += nrm * nrm;
    }
    const double betaSum = sumSqNorms - static_cast<double>(n) * S.squaredNorm();
    const double beta = betaSum / (static_cast<double>(n) * static_cast<double>(n) * static_cast<double>(p));
 
    // Optimal shrinkage coefficient, clipped to [0, 1]
    const double alpha = std::clamp(beta / delta, 0.0, 1.0);
 
    // Shrunk covariance: alpha * mu * I + (1 - alpha) * S
    MatrixXd covShrunk = (1.0 - alpha) * S;
    covShrunk.diagonal().array() += alpha * mu;
 
    return {covShrunk, alpha};
}