Numerical Stability

SPD operations involve eigendecomposition, matrix logarithms, and other operations that can be numerically sensitive. SPD Learn provides a comprehensive configuration system for managing numerical stability.

This document describes both the theoretical foundations and practical techniques for working with Symmetric Positive Definite (SPD) matrices in Riemannian geometry-based learning frameworks. The theoretical approaches are based on established methods in the field, including those described in the MENDR framework [Chen et al., 2025].

Overview

Working with SPD matrices presents numerical challenges:

1. Small eigenvalues: Operations like \(\log(\lambda)\) become undefined or unstable when eigenvalues approach zero.

2. Condition number: Ill-conditioned matrices (large ratio of max/min eigenvalues) cause precision loss in matrix operations.

3. Gradient computation: The Loewner matrix formulation requires careful handling of equal or nearly-equal eigenvalues.

4. Mixed precision: Half-precision (float16/bfloat16) training requires larger stability margins.

SPD Learn addresses these challenges with dtype-aware numerical thresholds that automatically adjust based on the precision of your computations.

See also

Geometric Concepts for the mathematical foundations of SPD matrices, eigendecomposition, Riemannian metrics, and geometric operations.

Covariance Matrix Rank Deficiency

The sample covariance matrix (SCM) from a data matrix \(\mathbf{X} \in \reals^{C \times T}\) (C channels, T time samples) is computed as:

\[\mathbf{SCM} = \frac{1}{T-1} \mathbf{X} \mathbf{X}^\top\]

This matrix is guaranteed to be symmetric positive semi-definite, but may have zero eigenvalues if \(T < C\). This rank deficiency is a common source of numerical instability when applying SPD operations that require strictly positive eigenvalues.

Numerical Stability Techniques

Trace Normalization with Epsilon Regularization

To ensure numerical stability during forward and backward passes, a two-stage regularization is applied:

Stage 1 - Pre-normalization regularization:

\[\mathbf{SCM}_{\text{reg}} = \mathbf{SCM} + \epsilon \mathbf{I}\]

Stage 2 - Trace normalization with post-regularization:

\[\mathbf{SCM}_{\text{norm}} = \frac{\mathbf{SCM}_{\text{reg}}}{\text{tr}(\mathbf{SCM}_{\text{reg}})} + \epsilon \mathbf{I}\]

where \(\epsilon = 10^{-5}\) is a typical choice and \(\I\) is the identity matrix.

Rationale:

• Pre-normalization \(\epsilon\) prevents division by near-zero traces

• Trace normalization ensures bounded eigenvalues

• Post-normalization \(\epsilon\) guarantees minimum eigenvalue \(\lambda_{\min} \geq \epsilon\)

Symmetrization

Floating-point operations can introduce small asymmetries. Explicit symmetrization ensures the SPD property:

\[\mathbf{X}_{\text{sym}} = \frac{\mathbf{X} + \mathbf{X}^\top}{2}\]

This should be applied after any operation that might introduce asymmetry.

Cholesky Decomposition for Gradient Flow

Instead of directly optimizing over the SPD manifold, a numerically stable approach is to parameterize SPD matrices via Cholesky decomposition:

\[\mathbf{M} = \mathbf{L} \mathbf{L}^\top\]

where \(\mathbf{L} \in \reals^{n \times n}\) is a learnable lower triangular matrix with positive diagonal elements.

Gradient computation via the product rule:

\[d\mathbf{M} = (d\mathbf{L}) \mathbf{L}^\top + \mathbf{L} (d\mathbf{L})^\top\]

This parameterization:

• Guarantees SPD output without explicit manifold projections

• Enables use of standard Euclidean optimizers (Adam, SGD)

• Provides stable gradient flow through the decomposition

SVD-Based Stable Differentiation

For operations requiring eigendecomposition (log, exp, power), using SVD provides numerical stability:

\[\mathbf{A} = \mathbf{U} \mathbf{S} \mathbf{V}^\top\]

For symmetric matrices, \(\mathbf{U} = \mathbf{V}\), and the singular values equal the absolute eigenvalues.

Gradient decomposition:

The gradient flow through SVD can be decomposed into:

1. Diagonal component: Gradients with respect to singular values

2. Off-diagonal component: Gradients with respect to singular vectors

The orthonormality constraints on \(\mathbf{U}\) and \(\mathbf{V}\) provide natural regularization and prevent gradient explosion.

Logarithmic Loss Functions

Operating on log-eigenvalues rather than raw eigenvalues prevents underflow/overflow:

Masked Autoencoding Loss:

\[\mathcal{L}_{\text{MAE}} = \| \log(\boldsymbol{\lambda}_{\text{masked}}) - \log(\hat{\boldsymbol{\lambda}}_{\text{masked}}) \|^2\]

where \(\boldsymbol{\lambda}\) denotes the vector of eigenvalues.

Benefits:

• Logarithmic scaling compresses the dynamic range

• Equal relative errors contribute equally to the loss

• Prevents gradient explosion from large eigenvalue differences

Condition Number Analysis

Definition and Significance

The condition number of an SPD matrix is:

\[\kappa(\mathbf{A}) = \frac{\lambda_{\max}}{\lambda_{\min}}\]

A high condition number indicates ill-conditioning, where:

• Small perturbations in input lead to large perturbations in output

• Numerical errors are amplified during matrix operations

• Gradient-based optimization becomes unstable

Condition Number Bounds

After trace normalization with \(\epsilon\)-regularization:

\[\kappa(\mathbf{A}_{\text{norm}}) \leq \frac{1}{\epsilon}\]

For \(\epsilon = 10^{-5}\), this bounds the condition number at \(10^5\).

Practical recommendation: Choose \(\epsilon\) to balance:

• Larger \(\epsilon\): Better numerical stability, but information loss

• Smaller \(\epsilon\): Preserves information, but risk of instability

Eigenvalue Perturbation Bounds

For a symmetric matrix \(\mathbf{A}\) with perturbation \(\mathbf{E}\):

Weyl’s Theorem:

\[|\lambda_i(\mathbf{A} + \mathbf{E}) - \lambda_i(\mathbf{A})| \leq \|\mathbf{E}\|_2\]

This bounds how much eigenvalues can change due to numerical errors bounded by \(\|\mathbf{E}\|_2\).

Configuration System

The Global Configuration

SPD Learn provides a global numerical_config object that controls all numerical stability thresholds:

from spd_learn.functional import numerical_config

# View current settings
print(numerical_config)

# Modify a threshold
numerical_config.eigval_clamp_scale = 1e5  # More conservative clamping

# Disable warnings
numerical_config.warn_on_clamp = False

Configuration Parameters

Parameter	Default	Description
eigval_clamp_scale	1e4	Scale for ReEig layer eigenvalue clamping
eigval_log_scale	1e2	Scale for matrix logarithm stability
eigval_sqrt_scale	1e2	Scale for matrix square root
eigval_inv_sqrt_scale	1e3	Scale for inverse square root
eigval_power_scale	1e3	Scale for matrix power operations
loewner_equal_scale	1e2	Scale for detecting equal eigenvalues
batchnorm_var_eps	1e-5	Epsilon for batch normalization scalar dispersion (absolute)
dropout_eps	1e-5	Epsilon for dropout diagonal entries (absolute)
warn_on_clamp	True	Emit warnings when eigenvalues are clamped

How thresholds are computed:

threshold = scale * torch.finfo(dtype).eps

For example, with eigval_clamp_scale=1e4 and dtype=torch.float32:

threshold = 1e4 * 1.19e-7  # ≈ 1.19e-3

Getting Epsilon Values

Use get_epsilon() to retrieve the appropriate threshold for an operation:

from spd_learn.functional import get_epsilon
import torch

# Get epsilon for eigenvalue clamping in float32
eps32 = get_epsilon(torch.float32, "eigval_clamp")
print(f"float32 clamp threshold: {eps32:.2e}")  # ~1.19e-3

# Get epsilon for float64 (tighter threshold)
eps64 = get_epsilon(torch.float64, "eigval_clamp")
print(f"float64 clamp threshold: {eps64:.2e}")  # ~2.22e-12

# Get epsilon for float16 (much larger threshold)
eps16 = get_epsilon(torch.float16, "eigval_clamp")
print(f"float16 clamp threshold: {eps16:.2e}")  # ~9.77e0

Temporary Configuration

Use NumericalContext to temporarily modify settings:

from spd_learn.functional import NumericalContext, get_epsilon
import torch

# Default threshold
print(f"Default: {get_epsilon(torch.float32, 'eigval_clamp'):.2e}")

# Temporarily use more conservative threshold
with NumericalContext(eigval_clamp_scale=1e6):
    print(f"Conservative: {get_epsilon(torch.float32, 'eigval_clamp'):.2e}")

# Back to default
print(f"Restored: {get_epsilon(torch.float32, 'eigval_clamp'):.2e}")

Checking SPD Validity

Use check_spd_eigenvalues() to validate matrices:

from spd_learn.functional import check_spd_eigenvalues
import torch

# Create a matrix and check its eigenvalues
A = torch.randn(3, 3, dtype=torch.float32)
A = A @ A.T + 0.1 * torch.eye(3)
eigvals = torch.linalg.eigvalsh(A)

is_valid, min_val, num_bad = check_spd_eigenvalues(eigvals)
print(f"Valid: {is_valid}, Min eigenvalue: {min_val:.2e}")

# Optionally raise an error
check_spd_eigenvalues(eigvals, raise_on_failure=True)

Safe Eigenvalue Clamping

Use safe_clamp_eigenvalues() for consistent clamping, important for activation functions like ReEig:

from spd_learn.functional import safe_clamp_eigenvalues
import torch

eigvals = torch.tensor([1e-10, 1e-5, 1e-3, 1.0])

# Clamp with dtype-aware threshold
clamped = safe_clamp_eigenvalues(eigvals, "eigval_log")
print(clamped)  # Small values will be clamped

# Also get mask of which values were clamped
clamped, mask = safe_clamp_eigenvalues(eigvals, "eigval_log", return_mask=True)
print(f"Clamped values at indices: {mask.nonzero().squeeze()}")

Implementation Guidelines

Recommended Numerical Thresholds

Parameter	Typical Value	Purpose
\(\epsilon_{\text{reg}}\)	\(10^{-5}\)	Diagonal regularization
\(\epsilon_{\text{log}}\)	\(10^{-7}\)	Logarithm argument floor
\(\epsilon_{\text{div}}\)	\(10^{-8}\)	Division stability
Max condition number	\(10^{6}\)	Ill-conditioning threshold

Stability Checklist

When implementing SPD matrix operations:

1. Always regularize before computing logarithms

2. Symmetrize after any matrix operation that might introduce asymmetry

3. Check eigenvalues in debug mode to detect near-singular matrices

4. Use double precision (float64) when possible for intermediate computations

5. Clip eigenvalues to \([\epsilon, \infty)\) before taking logarithms

6. Monitor condition numbers during training

Code Example

Pseudocode for stable SPD operations:

def stable_log_euclidean_mean(matrices, eps=1e-5):
    """Compute Log-Euclidean mean with numerical stability."""
    log_sum = 0
    for A in matrices:
        # Regularize
        A_reg = A + eps * eye(n)
        # Symmetrize
        A_sym = (A_reg + A_reg.T) / 2
        # Compute stable log
        eigvals, eigvecs = eigh(A_sym)
        eigvals = maximum(eigvals, eps)  # Clip eigenvalues
        log_A = eigvecs @ diag(log(eigvals)) @ eigvecs.T
        log_sum += log_A

    # Compute mean in tangent space
    log_mean = log_sum / len(matrices)

    # Map back to manifold
    eigvals, eigvecs = eigh(log_mean)
    mean = eigvecs @ diag(exp(eigvals)) @ eigvecs.T

    return (mean + mean.T) / 2  # Final symmetrization

Recommendations for Different Scenarios

Standard Training (float32)

The default settings are well-tuned for float32 training:

# No changes needed for most cases
from spd_learn.models import SPDNet

model = SPDNet(n_chans=64, n_outputs=4)

Ill-Conditioned Matrices

For matrices with high condition numbers (common in EEG/fMRI):

from spd_learn.functional import numerical_config

# Use more conservative clamping
numerical_config.eigval_clamp_scale = 1e5
numerical_config.eigval_log_scale = 1e3

# Or consider using float64
model = model.double()

Mixed Precision Training

For float16/bfloat16 training, be more conservative:

from spd_learn.functional import numerical_config, recommend_dtype_for_spd

# Check if float16 is appropriate
condition_number = 1e6  # Estimated from your data
recommended = recommend_dtype_for_spd(condition_number)
print(f"Recommended dtype: {recommended}")

# If using float16, increase scales
numerical_config.eigval_clamp_scale = 1e6
numerical_config.eigval_log_scale = 1e4

High Precision Requirements

For research or when maximum precision is needed:

import torch
from spd_learn.functional import numerical_config

# Use float64
model = model.double()

# Use tighter thresholds
numerical_config.eigval_clamp_scale = 1e2
numerical_config.eigval_log_scale = 1e1

Common Issues and Solutions

NaN Values During Training

Symptom: Loss becomes NaN after some epochs.

Cause: Usually due to eigenvalues becoming too small or negative.

Solution:

from spd_learn.functional import numerical_config

# 1. Enable warnings to see when clamping occurs
numerical_config.warn_on_clamp = True

# 2. Use more conservative thresholds
numerical_config.eigval_clamp_scale = 1e5

# 3. Add regularization to your covariance matrices
from spd_learn.modules import TraceNorm

trace_norm = TraceNorm(eps=1e-4)

Slow Convergence

Symptom: Model trains but converges slowly or gets stuck.

Cause: Overly conservative thresholds may clip important information.

Solution:

# Try tighter thresholds if your data is well-conditioned
numerical_config.eigval_clamp_scale = 1e3

# Check condition numbers of your data
import torch

cond_numbers = []
for batch in dataloader:
    cov = compute_covariance(batch)
    eigvals = torch.linalg.eigvalsh(cov)
    cond = eigvals.max() / eigvals.min()
    cond_numbers.append(cond)
print(f"Median condition number: {torch.median(torch.stack(cond_numbers))}")

Warnings About Eigenvalue Clamping

Symptom: Many warnings about eigenvalue clamping.

Cause: Your data has small eigenvalues being modified.

Options:

# Option 1: Disable warnings if this is expected
numerical_config.warn_on_clamp = False

# Option 2: Preprocess data with regularization
from spd_learn.modules import Shrinkage

shrinkage = Shrinkage(alpha=0.1)  # Ledoit-Wolf shrinkage

# Option 3: Use higher precision
model = model.double()

API Reference

spd_learn.functional.get_epsilon(dtype: dtype, name: Literal['eigval_clamp', 'eigval_log', 'eigval_sqrt', 'eigval_inv_sqrt', 'eigval_power', 'loewner_equal', 'batchnorm_var', 'dropout', 'trace_norm', 'stiefel_init', 'division_safe'] = 'eigval_clamp', *, config: NumericalConfig | None = None) → float

Get a dtype-aware epsilon value for numerical stability.

This function returns an appropriate epsilon value based on the data type and the intended use case. It scales the machine epsilon by a factor that ensures numerical stability for the specific operation.

Parameters:

• dtype (torch.dtype) – The PyTorch dtype to compute epsilon for.

• name (ThresholdName, default="eigval_clamp") –

  The type of threshold to compute. Options are:

  • "eigval_clamp": General eigenvalue clamping (ReEig layer)

  • "eigval_log": Eigenvalue clamping before log operation

  • "eigval_sqrt": Eigenvalue clamping before sqrt operation

  • "eigval_inv_sqrt": Eigenvalue clamping before inverse sqrt

  • "eigval_power": Eigenvalue clamping before power operation

  • "loewner_equal": Detection of equal eigenvalues in Loewner matrix

  • "batchnorm_var": Batch normalization variance epsilon

  • "dropout": Dropout diagonal epsilon

  • "trace_norm": Trace normalization epsilon

  • "stiefel_init": Stiefel manifold initialization

  • "division_safe": Safe division operations

• config (NumericalConfig, optional) – Configuration to use. If None, uses the global numerical_config.

Returns:

The computed epsilon value.

Return type:

float

Examples

>>> import torch
>>> from spd_learn.functional.numerical import get_epsilon
>>> # Get epsilon for float32 eigenvalue clamping
>>> eps32 = get_epsilon(torch.float32, "eigval_clamp")
>>> print(f"float32 eigval_clamp: {eps32:.2e}")
float32 eigval_clamp: 1.19e-03
>>> # Get epsilon for float64 (more precise)
>>> eps64 = get_epsilon(torch.float64, "eigval_clamp")
>>> print(f"float64 eigval_clamp: {eps64:.2e}")
float64 eigval_clamp: 2.22e-12
>>> # float16 needs larger epsilon
>>> eps16 = get_epsilon(torch.float16, "eigval_clamp")
>>> print(f"float16 eigval_clamp: {eps16:.2e}")
float16 eigval_clamp: 9.77e+00

See also

get_epsilon_tensor

Returns epsilon as a tensor on the correct device.

numerical_config

Global configuration for threshold scales.

spd_learn.functional.get_epsilon_tensor(dtype: dtype, name: Literal['eigval_clamp', 'eigval_log', 'eigval_sqrt', 'eigval_inv_sqrt', 'eigval_power', 'loewner_equal', 'batchnorm_var', 'dropout', 'trace_norm', 'stiefel_init', 'division_safe'] = 'eigval_clamp', *, device: str | device | None = None, config: NumericalConfig | None = None) → Tensor

Get a dtype-aware epsilon value as a tensor.

Similar to get_epsilon(), but returns a tensor on the specified device. This is useful when the epsilon needs to be used in tensor operations that require matching devices.

Parameters:

• dtype (torch.dtype) – The PyTorch dtype to compute epsilon for.

• name (ThresholdName, default="eigval_clamp") – The type of threshold to compute.

• device (str or torch.device, optional) – The device to place the tensor on. If None, uses CPU.

• config (NumericalConfig, optional) – Configuration to use. If None, uses the global numerical_config.

Returns:

A scalar tensor containing the epsilon value.

Return type:

torch.Tensor

Examples

>>> import torch
>>> from spd_learn.functional.numerical import get_epsilon_tensor
>>> eps = get_epsilon_tensor(torch.float32, "eigval_clamp", device="cpu")
>>> print(eps)
tensor(0.0012)

spd_learn.functional.safe_clamp_eigenvalues(eigenvalues: Tensor, name: Literal['eigval_clamp', 'eigval_log', 'eigval_sqrt', 'eigval_inv_sqrt', 'eigval_power', 'loewner_equal', 'batchnorm_var', 'dropout', 'trace_norm', 'stiefel_init', 'division_safe'] = 'eigval_clamp', *, config: NumericalConfig | None = None, return_mask: bool = False) → Tensor | tuple

Safely clamp eigenvalues with dtype-aware threshold.

This function clamps eigenvalues to ensure they are positive and numerically stable. It uses a dtype-aware threshold to balance stability and precision.

Parameters:

• eigenvalues (torch.Tensor) – The eigenvalues to clamp.

• name (ThresholdName, default="eigval_clamp") – The type of threshold to use.

• config (NumericalConfig, optional) – Configuration to use. If None, uses the global numerical_config.

• return_mask (bool, default=False) – If True, also return a boolean mask indicating which eigenvalues were clamped.

Returns:

The clamped eigenvalues. If return_mask=True, returns a tuple of (clamped_eigenvalues, clamped_mask).

Return type:

torch.Tensor or tuple

Examples

>>> import torch
>>> from spd_learn.functional.numerical import safe_clamp_eigenvalues
>>> eigvals = torch.tensor([1e-10, 1e-5, 1e-3, 1.0])
>>> clamped = safe_clamp_eigenvalues(eigvals, "eigval_log")
>>> print(clamped)
tensor([1.1921e-05, 1.1921e-05, 1.0000e-03, 1.0000e+00])

spd_learn.functional.check_spd_eigenvalues(eigenvalues: Tensor, name: Literal['eigval_clamp', 'eigval_log', 'eigval_sqrt', 'eigval_inv_sqrt', 'eigval_power', 'loewner_equal', 'batchnorm_var', 'dropout', 'trace_norm', 'stiefel_init', 'division_safe'] = 'eigval_clamp', *, config: NumericalConfig | None = None, raise_on_failure: bool = False) → tuple

Check if eigenvalues satisfy SPD requirements.

Parameters:

• eigenvalues (torch.Tensor) – The eigenvalues to check.

• name (ThresholdName, default="eigval_clamp") – The threshold to use for the positivity check.

• config (NumericalConfig, optional) – Configuration to use. If None, uses the global numerical_config.

• raise_on_failure (bool, default=False) – If True, raise an error when eigenvalues fail the check.

Returns:

A tuple of (is_valid, min_eigenvalue, num_below_threshold).

Return type:

tuple

Raises:

ValueError – If raise_on_failure=True and eigenvalues are not valid.

Examples

>>> import torch
>>> from spd_learn.functional.numerical import check_spd_eigenvalues
>>> eigvals = torch.tensor([1e-10, 0.1, 1.0])
>>> is_valid, min_val, num_bad = check_spd_eigenvalues(eigvals)
>>> print(f"Valid: {is_valid}, Min: {min_val:.2e}, Bad count: {num_bad}")
Valid: False, Min: 1.00e-10, Bad count: 1

spd_learn.functional.get_loewner_threshold(eigenvalues: Tensor, *, config: NumericalConfig | None = None) → float

Get threshold for detecting equal eigenvalues in Loewner matrix.

The Loewner matrix computation requires special handling when eigenvalues are equal or nearly equal. This function returns an appropriate threshold for detecting such cases.

Parameters:

• eigenvalues (torch.Tensor) – The eigenvalues (used to determine dtype).

• config (NumericalConfig, optional) – Configuration to use. If None, uses the global numerical_config.

Returns:

The threshold for eigenvalue equality detection.

Return type:

float

Notes

The threshold is computed as:

threshold = scale * max(1, |eigenvalues|.max()) * eps

This adaptive threshold accounts for the magnitude of eigenvalues, providing better numerical stability for matrices with large eigenvalues.

spd_learn.functional.is_half_precision(dtype: dtype) → bool

Check if dtype is half precision (float16 or bfloat16).

Parameters:

dtype (torch.dtype) – The dtype to check.

Returns:

True if the dtype is float16 or bfloat16.

Return type:

bool

spd_learn.functional.recommend_dtype_for_spd(condition_number: float, *, prefer_speed: bool = False) → dtype

Recommend a dtype based on expected matrix condition number.

Parameters:

• condition_number (float) – The expected condition number of the SPD matrices.

• prefer_speed (bool, default=False) – If True, prefer faster dtypes when possible.

Returns:

The recommended dtype.

Return type:

torch.dtype

Examples

>>> from spd_learn.functional.numerical import recommend_dtype_for_spd
>>> # Well-conditioned matrices can use float32
>>> print(recommend_dtype_for_spd(1e3))
torch.float32
>>> # Ill-conditioned matrices need float64
>>> print(recommend_dtype_for_spd(1e10))
torch.float64

class spd_learn.functional.NumericalConfig(eigval_clamp_scale: float = 10000.0, eigval_log_scale: float = 100.0, eigval_sqrt_scale: float = 100.0, eigval_inv_sqrt_scale: float = 1000.0, eigval_power_scale: float = 1000.0, loewner_equal_scale: float = 100.0, stiefel_init_scale: float = 1000.0, division_safe_scale: float = 100000.0, batchnorm_var_eps: float = 1e-05, dropout_eps: float = 1e-05, trace_norm_eps: float = 1e-06, warn_on_clamp: bool = True, strict_spd_check: bool = False, _threshold_cache: Dict[tuple, float]=<factory>)

Bases: object

Global configuration for numerical stability thresholds.

This class provides centralized control over numerical stability parameters used throughout the spd_learn library. All thresholds are specified as multipliers of the machine epsilon for the given dtype.

The actual threshold for a given dtype is computed as:

threshold = scale * torch.finfo(dtype).eps

For example, with eigval_clamp_scale=1e4 and dtype=torch.float32:

threshold = 1e4 * 1.19e-7 ≈ 1.19e-3

Parameters:

• eigval_clamp_scale (float) – Scale factor for general eigenvalue clamping (ReEig layer). Default: 1e4 (yields ~1e-3 for float32).

• eigval_log_scale (float) – Scale factor for eigenvalue clamping before log operation. Default: 1e2 (yields ~1e-5 for float32).

• eigval_sqrt_scale (float) – Scale factor for eigenvalue clamping before sqrt operation. Default: 1e2 (yields ~1e-5 for float32).

• eigval_inv_sqrt_scale (float) – Scale factor for eigenvalue clamping before inverse sqrt. Default: 1e3 (yields ~1e-4 for float32).

• eigval_power_scale (float) – Scale factor for eigenvalue clamping before power operation. Default: 1e3 (yields ~1e-4 for float32).

• loewner_equal_scale (float) – Scale factor for detecting equal eigenvalues in Loewner matrix. Default: 1e2 (yields ~1e-5 for float32).

• batchnorm_var_eps (float) – Absolute epsilon for batch normalization scalar dispersion. This is a scalar value (mean squared Frobenius norm in tangent space), not a variance matrix. Default: 1e-5.

• dropout_eps (float) – Absolute epsilon for dropout diagonal entries. Default: 1e-5.

• trace_norm_eps (float) – Absolute epsilon for trace normalization. Default: 1e-6.

• stiefel_init_scale (float) – Scale factor for Stiefel manifold initialization. Default: 1e3 (yields ~1e-4 for float32).

• division_safe_scale (float) – Scale factor for safe division operations. Default: 1e5 (yields ~1e-2 for float32).

• warn_on_clamp (bool) – Whether to emit warnings when eigenvalues are clamped. Default: True.

• strict_spd_check (bool) – Whether to perform strict SPD checks (slower but safer). Default: False.

Notes

The default scale factors are chosen to balance numerical stability with accuracy [Higham, 2002]. More conservative (larger) values provide better stability but may reduce precision. Less conservative (smaller) values preserve more information but risk numerical issues.

For mixed-precision training (fp16), consider using larger scale factors as the machine epsilon for fp16 is much larger (~9.77e-4).

batchnorm_var_eps: float = 1e-05

clear_cache() → None

Clear the threshold cache after configuration changes.

division_safe_scale: float = 100000.0

dropout_eps: float = 1e-05

eigval_clamp_scale: float = 10000.0

eigval_inv_sqrt_scale: float = 1000.0

eigval_log_scale: float = 100.0

eigval_power_scale: float = 1000.0

eigval_sqrt_scale: float = 100.0

get_scale(name: Literal['eigval_clamp', 'eigval_log', 'eigval_sqrt', 'eigval_inv_sqrt', 'eigval_power', 'loewner_equal', 'batchnorm_var', 'dropout', 'trace_norm', 'stiefel_init', 'division_safe']) → float

Get the scale factor for a given threshold name.

Parameters:

name (ThresholdName) – The name of the threshold.

Returns:

The scale factor for the threshold.

Return type:

float

is_absolute(name: Literal['eigval_clamp', 'eigval_log', 'eigval_sqrt', 'eigval_inv_sqrt', 'eigval_power', 'loewner_equal', 'batchnorm_var', 'dropout', 'trace_norm', 'stiefel_init', 'division_safe']) → bool

Check if a threshold uses absolute values (not scaled by eps).

Parameters:

name (ThresholdName) – The name of the threshold.

Returns:

True if the threshold is absolute, False if scaled.

Return type:

bool

loewner_equal_scale: float = 100.0

stiefel_init_scale: float = 1000.0

strict_spd_check: bool = False

summary(dtype: dtype = torch.float32) → str

Return formatted string showing all thresholds for a given dtype.

Parameters:

dtype (torch.dtype, default=torch.float32) – The dtype to compute thresholds for.

Returns:

Formatted summary of all threshold values.

Return type:

str

Examples

>>> from spd_learn.functional.numerical import numerical_config
>>> print(numerical_config.summary(torch.float32))
Numerical Configuration Summary (dtype=torch.float32)
==================================================
...

trace_norm_eps: float = 1e-06

warn_on_clamp: bool = True

class spd_learn.functional.NumericalContext(**kwargs)

Bases: object

Context manager for temporarily modifying numerical configuration.

This context manager allows temporary modification of the global numerical configuration. The original configuration is restored when exiting the context.

Parameters:

**kwargs – Configuration parameters to temporarily override.

Examples

>>> from spd_learn.functional.numerical import (
...     numerical_config, NumericalContext, get_epsilon
... )
>>> import torch
>>> # Default epsilon
>>> print(f"Default: {get_epsilon(torch.float32, 'eigval_clamp'):.2e}")
Default: 1.19e-03
>>> # Temporarily use more conservative threshold
>>> with NumericalContext(eigval_clamp_scale=1e6):
...     print(f"Conservative: {get_epsilon(torch.float32, 'eigval_clamp'):.2e}")
Conservative: 1.19e-01
>>> # Back to default
>>> print(f"Restored: {get_epsilon(torch.float32, 'eigval_clamp'):.2e}")
Restored: 1.19e-03

References

[1]

Nicholas J Higham. Accuracy and Stability of Numerical Algorithms. SIAM, 2nd edition, 2002. doi:10.1137/1.9780898718027.

[2]

Matthew Chen, Micky Nnamdi, Justin Shao, Andrew Hornback, Hongyun Huang, Ben Tamo, Yishan Zhong, Benoit Marteau, Wenqi Shi, and May Dongmei Wang. Mendr: manifold explainable neural data representations. 2025. arXiv:2508.04956.

See also

• Frequently Asked Questions — Troubleshooting common issues

• User Guide — Getting started with SPD Learn

• Theory — Theory section overview

• Geometric Concepts — Understanding SPD geometry