---
title: 'Riemstats: Inference for SPD Connectomes'
url: https://www.emergentmind.com/topics/riemstats
type: topic
---

# Riemstats: Inference for SPD Connectomes

Searching arXiv for the specified package paper and closely related context.
{"query":"ti:\"riemtan, riemstats\" OR abs:\"riemstats\"","max_results":5,"sort_by":"submittedDate","sort_order":"descending"}
{"query":"ti:\"geomstats\" AND abs:\"Riemannian geometry\"","max_results":5,"sort_by":"relevance","sort_order":"descending"}
**riemstats** is an R package for statistical inference on manifold-valued data, with a primary emphasis on symmetric positive definite (SPD) matrices arising in functional connectivity analysis and connectomics. It is presented together with the companion package **riemtan**, which supplies the geometric infrastructure, while **riemstats** provides the inferential layer for multiple-sample analysis, including **Fréchet ANOVA**, **Riemannian ANOVA** with classical MANOVA-style statistics, and **harmonization methods** for multi-site studies. The package is motivated by the observation that SPD matrices carry a non-Euclidean structure that standard Euclidean methods can distort, particularly in the computation of averages, distances, and group comparisons [2508.01342].

## 1. Definition and scope

The package is designed for analyses in which each observation is an SPD matrix, especially a brain connectome or functional connectivity matrix. In that setting, the central statistical problem is not merely to store or manipulate matrices, but to perform inference in a way that respects the Riemannian geometry of the SPD cone rather than flattening the data prematurely into Euclidean arrays. The paper therefore presents **riemstats** as the statistical component of a two-package system: **riemtan** handles geometry, representations, and data containers, while **riemstats** handles grouped inference and harmonization [2508.01342].

This division of labor is explicit. **riemtan** provides Riemannian logarithm and exponential maps, vectorization and unvectorization of tangent objects, conversion among manifold, tangent, and Euclidean-vector representations, Fréchet mean computation, metric abstractions, and classes for storing samples and grouped samples. **riemstats** then builds on these objects to implement multiple-sample procedures. In practical terms, the workflow described in the paper is: represent each subject’s connectome as an SPD matrix, choose a Riemannian metric, compute means and tangent-space coordinates using **riemtan**, and perform group comparison or harmonization using **riemstats**.

The implementation is centered on SPD matrices, although the paper also states that the broader conceptual target is manifold-valued data for which a Riemannian tangent-space workflow is meaningful. This suggests a degree of conceptual generality, but the concrete implementation and examples are clearly optimized for SPD connectome analysis.

## 2. Geometric framework on the SPD manifold

The central manifold is the space of \(p \times p\) symmetric positive definite matrices,
\[
SPD_p,
\]
together with the vector space of symmetric matrices,
\[
Symm_p.
\]
The paper states that \(SPD_p\) is an open subset of \(Symm_p\), hence a differentiable manifold, and that each tangent space is canonically isomorphic to \(Symm_p\):
\[
T_P SPD_p \cong Symm_p.
\]
This distinction is fundamental: a point on the manifold is a positive definite matrix, whereas a tangent vector at that point is a symmetric matrix [2508.01342].

The paper’s geometric motivation is that Euclidean operations on SPD matrices can ignore the curved geometry of the SPD cone, distort distances between connectomes, and produce averages or interpolations that are not geometrically meaningful. The Riemannian framework replaces Euclidean straight lines with geodesics, Euclidean means with Fréchet means, and ordinary vector representations with tangent-space coordinates defined by logarithm maps.

Several metrics are named as supported within the package design: the **Euclidean / Frobenius metric**, the **Affine Invariant Riemannian Metric (AIRM)**, the **log-Euclidean metric**, the **log-Cholesky metric**, and the **Bures-Wasserstein metric**. The paper emphasizes metric flexibility, but it also explicitly notes that it does **not** provide closed-form formulas in the text for the geodesic distances or exp/log maps under each of these metrics. That absence is itself part of the package’s characterization: the package is meant to be metric-modular, but the manuscript is more architectural and workflow-oriented than formula-complete.

The core geometric maps are described in standard terms. Given a reference point \(P \in SPD_p\), a matrix \(X \in SPD_p\) is mapped to the tangent space by
\[
\log_P(X),
\]
and a tangent vector is returned to the manifold by the exponential map. Tangent matrices are then converted into Euclidean coordinates through an isometric linear map
\[
\mathrm{vec}_P : T_P SPD_p \to \mathbb{R}^{p(p+1)/2}.
\]
The inverse is conceptually denoted by \(\mathrm{unvec}_P\). The paper does not print the half-vectorization formula explicitly, but it states that the vectorization is chosen to be an isometry, so tangent norms can be identified with Euclidean norms of the vectorized coordinates.

The Fréchet mean is described as the Riemannian analogue of the arithmetic mean. The intended definition is
\[
\mu_F = \arg\min_{P \in SPD_p} \frac{1}{n}\sum_{i=1}^n d^2(P,X_i),
\]
where \(d\) is the chosen geodesic distance. The associated sample variation is
\[
V = \frac{1}{n}\sum_{i=1}^n d^2(X_i,\hat\mu_F).
\]
These quantities underlie essentially all subsequent procedures in **riemstats**.

## 3. Statistical procedures for multiple samples

The package implements three principal inferential families: **Fréchet ANOVA**, **Riemannian ANOVA with classical test statistics**, and **harmonization methods for multi-site studies** [2508.01342].

Fréchet ANOVA is described as a metric-space analogue of ANOVA for testing equality of group Fréchet means. The null hypothesis is
\[
H_0:\ \mu_{F,1}=\mu_{F,2}=\cdots=\mu_{F,G}.
\]
The method is attributed to Dubey and Müller (2019). The paper states that the procedure uses sample Fréchet means and sample variations, constructs an \(F\)-ratio analogue, adds a second statistic \(U_n\) to account for scale differences, and produces a combined statistic \(T\) whose asymptotic null law satisfies
\[
T \overset{d}{\longrightarrow} \chi^2_1.
\]
In the package interface, `frechet_anova(super_sample)` returns a p-value. Because the method operates at the metric-space level, it is presented as relatively general, but the paper also notes that it may miss directional covariance information that a tangent-space method can exploit.

Riemannian ANOVA is the more differential-structure-aware alternative. Here the data are mapped into a tangent space, vectorized, and analyzed with classical MANOVA-style statistics adapted to the Riemannian setting. The paper explicitly states that **Wilks’ Lambda** and **Pillai’s trace** are implemented. In the examples, `log_wilks_lambda` appears as the default statistic and `pillais_trace` is also available. The manuscript does not print the determinant or trace formulas for the Riemannian versions, nor the full between-group and within-group covariance decompositions, but it characterizes the procedure as a tangent-space analogue of MANOVA that can exploit covariance anisotropy and may have better power than distance-based methods when covariance is non-spherical.

The third major component is harmonization for multi-site neuroimaging studies. The paper states that **ComBat harmonization** and **rigid harmonization** are implemented for SPD matrices. The functions shown are
```r
combat_harmonization(batch_sample)
```
and
```r
rigid_harmonization(batch_sample)
```
The conceptual basis is that connectomes from each site can be mapped to tangent spaces around site-specific Fréchet means and then adjusted either by manifold-aware transport or by Euclidean harmonization on tangent-space features. The manuscript explicitly discusses both direct manifold harmonization and feature harmonization as broader strategies, but it does **not** provide the exact correction equations for the package implementations. It also notes clustering-based evaluation criteria such as the **Calinski–Harabasz score**, **Davies–Bouldin index**, and **silhouette score**, though it does not present a full harmonization benchmark in the displayed text.

## 4. Software architecture and representation system

The software architecture is explicitly described as **modular**, **composition-based**, **high-level**, **extensible**, and **parallel-aware**. A central design principle is seamless movement among three linked representations of the same data: manifold points in \(SPD_p\), tangent images in \(T_PSPD_p\), and Euclidean vector coordinates [2508.01342].

The key abstraction is the S3 metric class `riem_met`, whose elements are lists containing four functions: `log`, `exp`, `vec`, and `unvec`. This design makes the package metric-modular: adding a new metric means supplying these four operations rather than rewriting downstream statistical procedures.

At the sample level, the paper describes the R6 class `CSample`, which stores three linked forms of a sample: manifold representation (`conns`), tangent representation (`tan_imgs`), and vectorized representation (`vec_imgs`). It also stores the sample size, matrix size, chosen metric, Fréchet mean, centering status, variation, and sample covariance. For grouped analysis, the package uses `CSuperSample`, another R6 class storing a list of `CSample` objects together with total sample size, matrix size, manifold dimension, common metric, pooled Fréchet mean, pooled covariance, and within-group and total covariance summaries.

The class `TangentImageHandler` manages conversion among manifold, tangent, and vectorized forms. This layered architecture is meant to permit geometry-aware internal computation while still exposing Euclidean features when standard statistical machinery is needed.

The paper also emphasizes numerical stability and memory efficiency. Inputs are expected to use packed SPD or symmetric matrix classes from the **Matrix** package, including `dpoMatrix`, `dspMatrix`, and `dppMatrix`. This is presented as a response to limitations in earlier tools; the paper specifically remarks that `Riemann` may sometimes output Fréchet means with negative eigenvalues, whereas the proposed implementation addresses such issues partly through its matrix-class design.

## 5. Algorithms and practical workflow

The typical workflow begins with SPD matrices derived from functional connectivity data, which are parsed into **Matrix** SPD classes and assembled into `CSample` objects. Several such samples are then combined into a `CSuperSample`, after which one may compute Fréchet means, center the data, obtain tangent-space vectors, and run either group comparison or harmonization [2508.01342].

The Fréchet mean estimator, although implemented in **riemtan**, is central to **riemstats**. The paper describes it as a mini-batch gradient descent method with configurable batch size, configurable learning rate with default \(0.2\), relative-change convergence criterion, default tolerance \(0.05\), default maximum iterations \(20\), and parallel relocation of tangent vectors when the reference point changes. The update rule is given in pseudocode form:
\[
\text{gradient} = \text{mean}(\text{tangent\_vectors\_in\_batch}),
\]
followed by
\[
\text{new\_ref\_point} = \text{Exp}(\text{ref\_point}, \text{lr} \times \text{gradient}).
\]
The procedure initializes at \(I\), computes tangent vectors, updates the reference point batch by batch, relocates tangent vectors into the new tangent space, and stops when the relative change falls below the tolerance.

The representative statistical functions are:
```r
frechet_anova(super_sample)
```
for Fréchet ANOVA,
```r
riem_anova(super_sample, stat_fun = log_wilks_lambda, den = 100)
```
for tangent-space Riemannian ANOVA, and
```r
combat_harmonization(batch_sample)
rigid_harmonization(batch_sample)
```
for harmonization. The paper strongly recommends centering around the **Fréchet mean** rather than the identity matrix \(I\), since tangent coordinates at \(I\) may distort results. The typical sequence is:
```r
sample$compute_f_mean()
sample$center()
sample$compute_vecs()
```
after which Euclidean analyses can be applied to the vectorized tangent images if desired.

## 6. Empirical evaluation, applications, and limitations

The empirical evaluation in the paper is concentrated more heavily on **riemtan** than on **riemstats**, but these results are directly relevant because **riemstats** depends on the same geometric infrastructure. The benchmark focuses on Fréchet mean runtime under **AIRM**, chosen because it is computationally demanding. For \(20 \times 20\) matrices with \(n=100\), `Riemann` is reported as slightly faster, with 3.245s versus 3.534s for `riemtan`. For larger problems, however, the new implementation becomes substantially faster; for \(40 \times 40\) matrices with \(n=400\), the reported runtimes are 78.828s for `Riemann` and 12.562s for `riemtan`, approximately a **6.3× speedup** [2508.01342].

The paper also reports that performance is not monotone in the number of cores: intermediate parallelization performs better than maximal parallelization, because excessive core allocation introduces overhead. Similarly, intermediate mini-batch sizes are described as best, with batch size 64 often outperforming smaller values such as 32 and sometimes competing well with full-batch 128. These findings are relevant because they condition the practical scalability of group-level inference in **riemstats**.

Application examples are drawn mainly from connectomics, including a Human Connectome Project workflow associated with Abbas et al. (2021). The paper states that preliminary experiments on HCP Young Adult data suggest that Riemannian ANOVA is appropriately calibrated and may outperform distance-based methods when covariance is non-spherical, although detailed tables are not shown in the text provided.

Several limitations are stated explicitly. Many formulas are not given in closed form, including metric formulas, exact test-statistic formulas, and harmonization equations. Direct timing and accuracy studies for **riemstats** methods are limited relative to the benchmarking of **riemtan**. Harmonization details are underspecified, and the manuscript does not provide a full API reference, package class diagram, or explicit availability information such as CRAN or GitHub URLs. The paper also makes clear that, despite conceptual breadth, the implementation is centered on SPD matrices and connectomics rather than being a fully general manifold-statistics platform.

Within the broader software landscape, this positioning is significant. Packages such as **geomstats** provide reusable manifold classes, Riemannian metrics, exponential and logarithm maps, geodesics, weighted Fréchet means, variance, and tangent PCA for machine learning workflows [1805.08308]. **riemstats**, by contrast, is specialized toward grouped SPD/connectome inference in R, with an emphasis on ANOVA-style testing and multi-site harmonization built atop the representation system supplied by **riemtan** [2508.01342]. This suggests that its distinctive contribution lies less in generic manifold computation than in the statistical analysis of multiple samples of connectomes under metric-flexible SPD geometry.

Source: https://www.emergentmind.com/topics/riemstats