Matrix functions via the eigendecomposition
Everything Riemannian is built from three matrix functions — square root, inverse, and logarithm — and they are all defined the same way. Any SPD matrix has an eigendecomposition C = U \Lambda U^{\top} with orthonormal eigenvectors U and strictly positive eigenvalues on the diagonal of \Lambda. A matrix function is applied by transforming the eigenvalues and leaving the eigenvectors alone.
f(C) = U \, f(\Lambda) \, U^{\top}, \qquad \log C = U (\log \Lambda) U^{\top}, \qquad C^{1/2} = U \Lambda^{1/2} U^{\top}Matrix functions act on eigenvalues only. Positive-definiteness guarantees the log is real and the square root is unique and SPD.
To apply a function like log or square root to a matrix, rotate to its natural axes, apply the function to each eigenvalue, then rotate back. Positive-definiteness guarantees this is well-defined — no log of a negative number.
- U
- The eigenvectors — the matrix's own natural axes.
- \Lambda
- The diagonal matrix of eigenvalues.
- f(\Lambda)
- Apply the function f to each eigenvalue individually.
- \log C,\ C^{1/2}
- The matrix logarithm and the matrix square root.
These matrix log and exp operations are the workhorse of every Riemannian formula that follows.
The affine-invariant Riemannian distance
The canonical distance on the SPD manifold is the affine-invariant Riemannian metric (AIRM). Its geodesic distance between two covariances has a closed form: whiten the second by the first, take the matrix log, and measure its Frobenius norm.
\delta_R(C_1, C_2) = \big\lVert \log\!\big( C_1^{-1/2} C_2 \, C_1^{-1/2} \big) \big\rVert_F = \Big( \sum_{i=1}^{n} \log^2 \lambda_i \Big)^{1/2}The AIRM geodesic distance, where the λᵢ are the eigenvalues of C₁⁻¹C₂ (the generalized eigenvalues of the pair).
The right distance between two covariances first whitens one by the other, then measures how far the result sits from the identity in log-space. Identical matrices give distance 0.
- C_1^{-1/2} C_2 \, C_1^{-1/2}
- C_2 viewed from C_1's frame — the whitening step.
- \lambda_i
- The generalized eigenvalues of the pair (of C_1^{-1}C_2).
- \log \lambda_i
- How far each eigenvalue is from 1 (i.e. from equality).
- \delta_R(C_1, C_2)
- The affine-invariant Riemannian distance.
This affine-invariant metric measures true geodesic distance on the curved space of covariances.
Read the right-hand form physically. The \lambda_i are the eigenvalues of C_1^{-1} C_2 — the ratio of variance between the two trials along each common direction. If the trials are identical every ratio is 1, every log is 0, the distance is 0. A direction where C_2 has twice the variance contributes the same as one where it has half — because \log 2 and \log \tfrac{1}{2} have equal magnitude. The metric measures multiplicative discrepancy, which is exactly right for power features that live on a ratio scale.
Geodesics on the cone
The shortest path between two covariances is not the straight line C_1 + t(C_2 - C_1) — that chord can even leave the cone. The true geodesic curves along the manifold and stays SPD at every t.
\gamma(t) = C_1^{1/2} \big( C_1^{-1/2} C_2 \, C_1^{-1/2} \big)^{t} \, C_1^{1/2}, \qquad t \in [0, 1]The AIRM geodesic from C₁ (t=0) to C₂ (t=1). The midpoint t=½ is the geometric mean of the pair.
The shortest curved path between two covariance matrices. Slide t from 0 to 1 to smoothly morph C_1 into C_2, and the halfway point is exactly their geometric mean.
- t
- Position along the path, from 0 to 1.
- \gamma(t)
- The covariance matrix at position t on the path.
- C_1^{1/2}
- The square root that frames the path around C_1.
- t = \tfrac{1}{2}
- The midpoint — the geometric mean of the pair.
At t=0 you are at C_1, at t=1 at C_2, and at t=\tfrac12 at the swelling-free average.
Set t = \tfrac{1}{2} and \gamma returns the geometric mean of the two matrices — the swelling-free centre from guide 1. This same interpolation is the basis of the iterative mean and the tangent-space projection in guide 3, so the geodesic is not a curiosity: it is the operation every later algorithm is quietly running.
Affine invariance — the property that pays off
Here is the property the metric is named for, and the reason the whole track exists.
\delta_R\big( A C_1 A^{\top}, \; A C_2 A^{\top} \big) = \delta_R(C_1, C_2), \qquad \forall A \in GL(n)Affine (congruence) invariance: any invertible linear transform of the sensors leaves every AIRM distance unchanged.
The Riemannian distance does not care about any invertible re-mixing of your sensors: re-reference, rescale, or linearly recombine the channels and every distance stays exactly the same. That is why decoders transfer across setups.
- A
- Any invertible linear transform of the sensors.
- A C A^{\top}
- How a covariance matrix transforms under A.
- GL(n)
- The group of all invertible n\times n matrices.
- \delta_R
- The Riemannian distance, left unchanged by A.
Swap to a different amplifier gain or reference montage and the classifier's geometry is untouched — a huge win for cross-session robustness.
Congruence transforms C \mapsto A C A^{\top} are precisely what a change of measurement does: re-referencing, per-channel gain differences, a different electrode montage, or the unknown linear mixing by which cortical sources reach the scalp. AIRM is blind to all of them. So a decoder built on AIRM distances is automatically invariant to a large family of nuisances that plague EEG — no explicit correction, no calibration term, it is baked into the geometry.
The cost: computation and estimation
Nothing is free. Each distance needs an eigendecomposition, so it costs O(n^3) per matrix — fine for the tens of channels typical in motor-imagery EEG, a real concern for high-density or ECoG grids (guide 5). More subtly, the metric only makes sense if C is genuinely positive-definite, which requires a good estimate.
When time samples are few relative to channels, the raw sample covariance is ill-conditioned or singular. The standard remedy is Ledoit–Wolf shrinkage: pull the sample covariance toward a scaled identity by a data-driven amount, guaranteeing positive-definiteness and stabilizing every downstream log and inverse. In practice a shrinkage step in covariance estimation is not optional — it is what keeps the whole Riemannian machine numerically alive.