Covariance estimators

yetanotherspdnet.functions.m_estimators — from samples \(x_1, \dots, x_n \in \mathbb{R}^p\) to a scatter matrix. The sample covariance is \(\frac1n \sum_i x_i x_i^\top\); an M-estimator solves the fixed point

\[ \Sigma = \frac1n \sum_{i=1}^n u\!\left(x_i^\top \Sigma^{-1} x_i\right) x_i x_i^\top, \]

whose weight \(u\) down-weights outlying samples (Tyler: \(u(q) = p/q\); Student-t: \(u(q) = (p+\nu)/(\nu+q)\); Huber). m_estimator differentiates the unrolled iterations; MEstimator differentiates the fixed point implicitly, with a memory cost independent of the number of iterations. The layers SampleCovariance and MEstimation (Layers) wrap them.

autograd path	manual backward
tyler_function()	–	Tyler weight :math:u(q) = p / q.
student_function()	–	Student-t weight :math:u(q) = (p + \nu) / (\nu + q).
huber_function()	–	Huber weight: :math:u(q) = 1/\beta if :math:q \le \delta, else :math:\delta / (\beta q).
normalize_trace()	–	Scale SPD matrices so that :math:\operatorname{tr}(\Sigma) = p.
normalize_determinant()	–	Scale SPD matrices so that :math:\det(\Sigma) = 1.
sample_covariance()	–	Sample covariance matrix (SCM).
m_estimator_step()	–	One fixed-point iteration :math:\Sigma \mapsto F(\Sigma) of an M-estimator.
m_estimator()	MEstimator	M-estimator of scatter by fixed-point iterations (autograd path).
tyler_estimator()	–	Tyler’s M-estimator of scatter, normalized to fix its scale.
student_estimator()	–	Student-t M-estimator of scatter.

Details

Robust covariance estimation: sample covariance and M-estimators.

An M-estimator of scatter is a fixed point of

\[\Sigma = F(\Sigma) = \frac{1}{n} \sum_{i=1}^{n} u\big(x_i^\top \Sigma^{-1} x_i\big)\, x_i x_i^\top\]

for a weight function \(u\) that down-weights samples with a large Mahalanobis distance: Tyler (\(u(q) = p/q\)), Student-t (\(u(q) = (p + \nu)/(\nu + q)\)) or Huber. The sample covariance matrix corresponds to \(u \equiv 1\).

Two gradient paths are available, as elsewhere in the library:

• m_estimator() unrolls the fixed-point iterations and lets autograd differentiate through them (memory grows with the number of iterations);

• MEstimator differentiates the fixed point implicitly: the backward solves the adjoint equation \(w = g + J_\Sigma^\top w\) by iteration and returns \(J_X^\top w\), with memory independent of the number of iterations.

Tyler’s weight is scale-invariant (\(F(c\Sigma) = cF(\Sigma)\)), so its fixed point is only defined up to scale: use it with normalize="trace" or normalize="determinant", which is applied at every iteration and pins the scale down.

tyler_function(quadratic, n_features)

Tyler weight \(u(q) = p / q\).

Parameters:

• quadratic (torch.Tensor of shape (..., n_samples)) – Squared Mahalanobis distances \(q_i = x_i^\top \Sigma^{-1} x_i\)

• n_features (int) – Dimension \(p\) of the samples

Returns:

weights (torch.Tensor of shape (..., n_samples)) – Sample weights

Return type:

Tensor

student_function(quadratic, n_features, nu)

Student-t weight \(u(q) = (p + \nu) / (\nu + q)\).

Maximum-likelihood weight for a multivariate Student-t distribution with \(\nu\) degrees of freedom; tends to the sample covariance when \(\nu \to \infty\) and to Tyler’s weight when \(\nu \to 0\).

Parameters:

• quadratic (torch.Tensor of shape (..., n_samples)) – Squared Mahalanobis distances

• n_features (int) – Dimension \(p\) of the samples

• nu (float) – Degrees of freedom

Returns:

weights (torch.Tensor of shape (..., n_samples)) – Sample weights

Return type:

Tensor

huber_function(quadratic, delta, beta)

Huber weight: \(u(q) = 1/\beta\) if \(q \le \delta\), else \(\delta / (\beta q)\).

Parameters:

• quadratic (torch.Tensor of shape (..., n_samples)) – Squared Mahalanobis distances

• delta (float) – Threshold above which samples are down-weighted

• beta (float) – Scaling factor (makes the estimator consistent for Gaussian data)

Returns:

weights (torch.Tensor of shape (..., n_samples)) – Sample weights

Return type:

Tensor

normalize_trace(data)

Scale SPD matrices so that \(\operatorname{tr}(\Sigma) = p\).

Parameters:

data (torch.Tensor of shape (..., n_features, n_features)) – SPD matrices

Returns:

normalized (torch.Tensor of shape (..., n_features, n_features)) – SPD matrices with trace equal to n_features

Return type:

Tensor

normalize_determinant(data)

Scale SPD matrices so that \(\det(\Sigma) = 1\).

The determinant is computed through its logarithm to avoid overflow.

Parameters:

data (torch.Tensor of shape (..., n_features, n_features)) – SPD matrices

Returns:

normalized (torch.Tensor of shape (..., n_features, n_features)) – SPD matrices with unit determinant

Return type:

Tensor

sample_covariance(data, assume_centered=False)

Sample covariance matrix (SCM).

\[\hat\Sigma = \frac{1}{n'} \sum_{i=1}^{n} (x_i - \bar x)(x_i - \bar x)^\top\]

with \(n' = n - 1\) when the data are centered here, \(n'= n\) (and \(\bar x = 0\)) when assume_centered.

Parameters:

• data (torch.Tensor of shape (..., n_samples, n_features)) – Samples

• assume_centered (bool, optional) – Whether the data are already centered. Default is False

Returns:

covariance (torch.Tensor of shape (..., n_features, n_features)) – Sample covariance matrices

Return type:

Tensor

m_estimator_step(covariance, data, weight_function, shrinkage=None, normalize=None)

One fixed-point iteration \(\Sigma \mapsto F(\Sigma)\) of an M-estimator.

\[F(\Sigma) = \frac{1}{n} \sum_{i=1}^{n} u\big(x_i^\top \Sigma^{-1} x_i\big)\, x_i x_i^\top\]

followed, if requested, by the shrinkage \(\beta F(\Sigma) + (1 - \beta) I\) and by a normalization. The Mahalanobis distances are computed with a Cholesky factorization.

Parameters:

• covariance (torch.Tensor of shape (..., n_features, n_features)) – Current estimate

• data (torch.Tensor of shape (..., n_samples, n_features)) – Centered samples

• weight_function (Callable) – Weight \(u\), called on the tensor of squared distances of shape (..., n_samples)

• shrinkage (float | None, optional) – Shrinkage coefficient \(\beta \in (0, 1]\) towards the identity. Default is None (no shrinkage)

• normalize (Callable | None, optional) – Normalization applied to the result (e.g. normalize_trace()). Default is None

Returns:

covariance (torch.Tensor of shape (..., n_features, n_features)) – Updated estimate

Return type:

Tensor

m_estimator(

data,

weight_function,

n_iterations=30,

tol=1e-06,

assume_centered=False,

init=None,

shrinkage=None,

normalize=None,

)

M-estimator of scatter by fixed-point iterations (autograd path).

Iterates m_estimator_step() from init until the largest relative Frobenius change over the batch falls below tol or n_iterations is reached. Gradients flow through the unrolled iterations.

Parameters:

• data (torch.Tensor of shape (..., n_samples, n_features)) – Samples

• weight_function (Callable) – Weight \(u\) of squared Mahalanobis distances, e.g. functools.partial(student_function, n_features=p, nu=3.0)

• n_iterations (int, optional) – Maximum number of iterations. Default is 30

• tol (float, optional) – Stopping tolerance on the relative change. Default is 1e-6

• assume_centered (bool, optional) – Whether the data are already centered. Default is False

• init (torch.Tensor of shape (n_features, n_features) or (..., n_features, n_features), optional) – Initial estimate. Default is the identity

• shrinkage (float | None, optional) – Shrinkage coefficient towards the identity, applied at each iteration. Default is None

• normalize (str | None, optional) – "trace", "determinant" or None, applied at each iteration. Required for scale-invariant weights such as Tyler’s. Default is None

Returns:

covariance (torch.Tensor of shape (..., n_features, n_features)) – Estimated scatter matrices

Return type:

Tensor

class MEstimator(*args, **kwargs)

M-estimator of scatter with an implicit (fixed-point) backward.

The forward pass iterates without building a graph. At the fixed point \(\Sigma^\star = F(\Sigma^\star, X)\), the implicit function theorem gives \(\partial L / \partial X = J_X^\top w\) where \(w\) solves \(w = g + J_\Sigma^\top w\) (\(g\) the incoming gradient). The adjoint equation is solved by fixed-point iteration, which converges at the rate of the forward iteration since \(J_\Sigma\) has spectral radius below one at an attracting fixed point. The vector-Jacobian products of one step are obtained with autograd.

Use as MEstimator.apply(data, weight_function, n_iterations, tol, assume_centered, init, shrinkage, normalize) with the arguments of m_estimator().

static forward(

ctx,

data,

weight_function,

n_iterations=30,

tol=1e-06,

assume_centered=False,

init=None,

shrinkage=None,

normalize=None,

)

Forward pass: fixed-point iterations without graph

Parameters:

• ctx (torch.autograd.function._ContextMethodMixin) – Context object to save tensors for the backward pass

• data (Tensor) – See m_estimator()

• weight_function (Callable) – See m_estimator()

• n_iterations (int) – See m_estimator()

• tol (float) – See m_estimator()

• assume_centered (bool) – See m_estimator()

• init (Tensor | None) – See m_estimator()

• shrinkage (float | None) – See m_estimator()

• normalize (str | None) – See m_estimator()

Returns:

covariance (torch.Tensor of shape (..., n_features, n_features)) – Estimated scatter matrices

Return type:

Tensor

static backward(ctx, grad_output)

Backward pass: implicit differentiation at the fixed point

Parameters:

• ctx (torch.autograd.function._ContextMethodMixin) – Context object with the saved tensors

• grad_output (torch.Tensor of shape (..., n_features, n_features)) – Gradient of the loss with respect to the estimate

Returns:

grad_data (torch.Tensor of shape (..., n_samples, n_features)) – Gradient of the loss with respect to the samples; None for the other arguments

Return type:

tuple

tyler_estimator(

data,

n_iterations=30,

tol=1e-06,

assume_centered=False,

normalize='trace',

use_autograd=False,

)

Tyler’s M-estimator of scatter, normalized to fix its scale.

Parameters:

• data (torch.Tensor of shape (..., n_samples, n_features)) – Samples (at least n_features + 1 of them)

• n_iterations (int, optional) – Maximum number of iterations. Default is 30

• tol (float, optional) – Stopping tolerance. Default is 1e-6

• assume_centered (bool, optional) – Whether the data are already centered. Default is False

• normalize (str, optional) – "trace" or "determinant". Default is "trace"

• use_autograd (bool, optional) – Unrolled autograd path (True) or implicit backward (False, default)

Returns:

covariance (torch.Tensor of shape (..., n_features, n_features)) – Tyler estimates

Return type:

Tensor

student_estimator(

data,

nu,

n_iterations=30,

tol=1e-06,

assume_centered=False,

use_autograd=False,

)

Student-t M-estimator of scatter.

Parameters:

• data (torch.Tensor of shape (..., n_samples, n_features)) – Samples

• nu (float) – Degrees of freedom (\(\nu > 0\))

• n_iterations (int, optional) – Maximum number of iterations. Default is 30

• tol (float, optional) – Stopping tolerance. Default is 1e-6

• assume_centered (bool, optional) – Whether the data are already centered. Default is False

• use_autograd (bool, optional) – Unrolled autograd path (True) or implicit backward (False, default)

Returns:

covariance (torch.Tensor of shape (..., n_features, n_features)) – Student-t estimates

Return type:

Tensor