Skip to contents

Calculates the Riemannian distance between two symmetric positive-definite (SPD) matrices S1 and S2, using the Log-Euclidean metric as described by Mitteroecker & Bookstein (2008) and Arsigny et al. (2006). `d(S1, S2) = ||logm(S1) - logm(S2)||_F` which, for the M&B formulation, is effectively `sqrt(sum(log(lambda_j(S2^-1 S1))^2))`. This function is a core component of RGEOM-001.

Usage

riemannian_distance_spd(
  S1,
  S2,
  regularize_epsilon = 1e-06,
  eigenvalue_floor = 1e-09
)

Arguments

S1

A numeric, symmetric matrix (p x p).

S2

A numeric, symmetric matrix (p x p).

regularize_epsilon

A small positive value used for regularization to ensure matrices are SPD before inversion and eigenvalue computation. Default 1e-6.

eigenvalue_floor

A small positive value to threshold relative eigenvalues before taking the logarithm, preventing issues with `log(0)`. Default 1e-9.

Value

The Riemannian distance (a scalar). Returns `NA` if computation fails.

Examples

S1 <- matrix(c(2.3, -0.3, -0.3, 3.6), 2, 2)
S2 <- matrix(c(3.7, 1.9, 1.9, 2.8), 2, 2)
# riemannian_distance_spd(S1, S2) # Expected: ~1.24156 (Need to import MASS for ginv)

# Identical matrices
# riemannian_distance_spd(S1, S1) # Expected: 0