---
title: "Lecture 17: Riemannian Manifolds — CLTs for Intrinsic Fréchet Means"
subtitle: "Asymptotic normality on manifolds"
format:
html:
code-fold: true
code-tools: true
code-copy: true
pdf:
documentclass: scrartcl
pdf-engine: xelatex
toc: true
number-sections: true
geometry:
- margin=1in
colorlinks: true
bibliography: source/ref.bib
---
## Learning Goals
By the end of this lecture, learners should be able to:
- State the Bhattacharya–Patrangenaru (2005) CLT for intrinsic Fréchet means on manifolds.
- Interpret the sandwich covariance $\Lambda^{-1}C(\Lambda^\top)^{-1}$ in terms of the Hessian of the Fréchet function.
- Explain why the cut locus poses a fundamental challenge for CLTs on manifolds.
- Describe the Kendall–Le (2011) non-i.i.d. CLT and its Lindeberg condition.
- Construct asymptotic confidence regions for Fréchet means on the SPD manifold using the CLT.
- Compare manifold-based inference (Wald tests, CLT confidence regions) with metric-space inference (permutation tests, bootstrap) and explain when the manifold structure adds power.
## The Classical Manifold CLT
Let $X_1, X_2, \ldots$ be i.i.d. on an $m$-dimensional Riemannian manifold $(\mathcal{M}, g)$. The asymptotic distribution of $\hat{\mu}_n$ is controlled by the local curvature of the Fréchet function $F(x) = \mathbb{E} d^2(x, X)$ at its minimizer $\mu$.
::: {#thm-bp2005-clt .theorem title="Classical intrinsic-manifold CLT (Bhattacharya & Patrangenaru 2005, Theorem 2.2)"}
Suppose $\operatorname{supp}(Q)$ lies in a closed geodesic ball $\overline B(x_0,r)$ that is disjoint from $\operatorname{Cut}(x_0)$, and
$$
r<\frac{\pi}{4K},
$$
where $K^2$ is the supremum of the sectional curvatures on the ball when that supremum is positive (and $K=0$ otherwise, so the radius restriction is vacuous). Let $\mu$ be the resulting unique intrinsic mean, $\hat{\mu}_n$ a measurable sample intrinsic mean, $\varphi_0 = \exp_{x_0}^{-1}$, $\nu_0 = \varphi_0(\mu)$, $Y=\varphi_0(X)$, and $h_0(\nu,y)=d^2\{\varphi_0^{-1}(\nu),\varphi_0^{-1}(y)\}$. Define
$$
\Lambda_0 = \bigl[\mathbb{E} D_{r,r'} h_0(\nu_0, Y)\bigr]_{r,r'=1}^m, \qquad
C_0 = \operatorname{cov}\bigl(D_1 h_0(\nu_0, Y), \ldots, D_m h_0(\nu_0, Y)\bigr).
$$
Then
$$
\sqrt{n}\,\{\varphi_0(\hat{\mu}_n) - \varphi_0(\mu)\} \xrightarrow{\mathcal{D}} N(0, \Lambda_0^{-1} C_0 (\Lambda_0^\top)^{-1}).
$$
:::
::: {#cor-bp2005-centered .corollary title="Mean-centered special case"}
If $x_0 = \mu$ (anchor at the mean), then with $Z = \log_\mu(X) = \exp_\mu^{-1}(X)$,
$$
\sqrt{n}\,\log_\mu(\hat{\mu}_n) \xrightarrow{\mathcal{D}} N(0, 4\Lambda^{-1}\operatorname{cov}(Z)\Lambda^{-1}),
$$
where $\Lambda = [\mathbb{E} D_{r,r'} h(0, X)]$ is the expected Hessian of the squared-distance function at $\mu$. Moreover, $D_1 h(0, X) = -2Z$, so the score vector is proportional to the log map.
:::
The sandwich covariance $\Lambda^{-1}C(\Lambda^\top)^{-1}$ is the usual $M$-estimator sandwich: $\Lambda$ is the expected Hessian (local curvature of the objective), while $C$ is the covariance of its gradient. There is no general information identity relating these two matrices for squared-distance loss.
An important practical consequence of the BP CLT is the construction of asymptotic confidence regions for the intrinsic mean $\mu$. Since $\sqrt{n}\,\log_\mu(\hat{\mu}_n) \xrightarrow{\mathcal{D}} \mathcal{N}(0, \Sigma)$ with $\Sigma = 4\Lambda^{-1}\operatorname{cov}(Z)\Lambda^{-1}$, Slutsky's lemma yields the pivot
$$
n \cdot \log_{\hat{\mu}_n}(\mu)^\top \hat{\Sigma}^{-1} \log_{\hat{\mu}_n}(\mu) \xrightarrow{\mathcal{D}} \chi^2_m,
$$
where $\hat{\Sigma}$ is a consistent estimator of $\Sigma$ (e.g., the empirical sandwich estimator or a bootstrap variance estimate) and $m = \dim\mathcal{M}$. A $100(1-\alpha)\%$ **confidence region in the tangent space** $T_{\hat{\mu}_n}\mathcal{M}$ is therefore the ellipsoid
$$
C_{1-\alpha}^{\text{tan}} = \left\{ v \in T_{\hat{\mu}_n}\mathcal{M} : n \cdot v^\top \hat{\Sigma}^{-1} v \le \chi^2_{m,1-\alpha} \right\}.
$$
To obtain the **confidence region on the manifold itself**, the exponential map at $\hat{\mu}_n$ transports this ellipsoid from the tangent space onto $\mathcal{M}$:
$$
C_{1-\alpha}^{\mathcal{M}} = \exp_{\hat{\mu}_n}\!\left(C_{1-\alpha}^{\text{tan}}\right) = \left\{ \exp_{\hat{\mu}_n}(v) : v \in T_{\hat{\mu}_n}\mathcal{M},\; n \cdot v^\top \hat{\Sigma}^{-1} v \le \chi^2_{m,1-\alpha} \right\}.
$$
Here $\hat\Sigma$ and the log vector must be represented in the same basis of $T_{\hat\mu_n}\mathcal M$. The exponential map preserves radial geodesic distance but distorts transverse lengths and volume according to curvature, so the image need not remain an ordinary Euclidean ellipsoid. For a shrinking $O(n^{-1/2})$ region, this distortion is of smaller order. In practice, one can use the empirical sandwich estimator or a nonparametric bootstrap; this construction is the manifold version of Corollary 2.2 in @BhattacharyaPatrangenaru2005.
## The Kendall–Le Non-i.i.d. CLT
@KendallLe2011 extend the CLT to independent non-identically distributed observations under a Lindeberg condition.
::: {#thm-kendall-le .theorem title="Non-i.i.d. central approximation (Kendall & Le 2011, Theorem 4)"}
Let $X_1,X_2,\ldots$ be independent, not necessarily identically distributed, on a complete connected Riemannian manifold. Let $\mu$ be a common local Fréchet mean, assume $\mathbb P\{X_i\in\operatorname{Cut}(\mu)\}=0$, and let $\hat\mu_n$ be a consistent measurable local empirical mean. Define
$$
Y_i=\log_\mu(X_i),\quad
\phi_n(\mu)=\frac12\sum_{i=1}^n\mathbb E\|Y_i\|^2,\quad
V_n=\operatorname{cov}\!\left\{\frac{\sum_{i=1}^nY_i}{\sqrt{2\phi_n(\mu)}}\right\},
$$
and let $H_i(\mu)$ be the Hessian of $\frac12d^2(\cdot,X_i)$ at $\mu$. In addition to local Hessian regularity and moment bounds, assume $\phi_n(\mu)$ grows at least linearly, the normalized expected Hessian
$$
\widetilde H_n=\frac{\sum_{i=1}^n\mathbb EH_i(\mu)}{2\phi_n(\mu)}
$$
is asymptotically nonsingular, and, for every $\varepsilon>0$,
$$
\frac1{\phi_n(\mu)}
\sum_{i=1}^n
\mathbb E\!\left[
\|Y_i\|^2;\ \|Y_i\|^2>\varepsilon\phi_n(\mu)
\right]\longrightarrow0.
$$
Then the law of $\sqrt{2\phi_n(\mu)}\,\log_\mu(\hat\mu_n)$ is asymptotically close to
$$
N\!\left(0,\widetilde H_n^{-1}V_n\widetilde H_n^{-1}\right).
$$
This is a central *approximation* with an $n$-dependent Gaussian covariance. It becomes an ordinary weak-convergence statement if that covariance converges.
:::
## Interactive Exploration: CLT on $S^2$
The $S^1$ demo above establishes the CLT on the circle, a one-dimensional manifold. We now extend to the **2-sphere** $S^2 = \{x \in \mathbb{R}^3 : \|x\| = 1\}$, a two-dimensional Riemannian manifold of constant curvature $+1$. On $S^2$, the tangent space at any point is a plane ($\cong \mathbb{R}^2$), and the CLT predicts that $\sqrt{n}\,\log_\mu(\hat{\mu}_n)$ converges to a bivariate normal $\mathcal{N}(0, \Sigma)$.
The demo below simulates data from a **von Mises–Fisher (VMF) distribution** on $S^2$ concentrated at the north pole $\mu_0 = (0,0,1)^\top$. The **left panel** shows one representative sample on the sphere together with the estimated Fréchet mean $\hat{\mu}_n$ and the asymptotic confidence region obtained by mapping the CLT ellipsoid from $T_{\hat{\mu}_n}S^2$ back to the sphere via $\exp_{\hat{\mu}_n}$. The **right panel** shows the empirical distribution of $\sqrt{n}\,\log_{\mu_0}(\hat{\mu}_n)$ across Monte Carlo replications — the CLT predicts convergence to a bivariate Gaussian.
```{ojs}
//| label: fig-clt-s2
//| fig-cap: "Interactive: CLT on S² — data, Fréchet mean, and asymptotic confidence regions"
s2n_ctrl = Inputs.range([10, 200], {step: 10, value: 50, label: "Sample size n"})
s2kap_ctrl = Inputs.range([0.5, 10], {step: 0.5, value: 4, label: "Concentration κ"})
s2nmc_ctrl = Inputs.range([50, 300], {step: 50, value: 200, label: "MC replications"})
s2conf_ctrl = Inputs.range([0.80, 0.99], {step: 0.01, value: 0.95, label: "Confidence level 1−α"})
s2n = Generators.input(s2n_ctrl)
s2kap = Generators.input(s2kap_ctrl)
s2nmc = Generators.input(s2nmc_ctrl)
s2conf = Generators.input(s2conf_ctrl)
s2ctrls = html`
<style>
.s2g{display:flex;flex-wrap:wrap;gap:6px 20px;width:100%;margin:0 auto 12px;font:0.85em system-ui,sans-serif;container-type:inline-size}
.s2g>*{flex:1 1 calc((100% - 40px)/3);min-width:0;margin:0}
.s2g input[type="number"]{width:7.5rem !important}
@container(max-width:700px){.s2g>*{flex-basis:calc((100% - 20px)/2)}}
@container(max-width:480px){.s2g>*{flex-basis:100%}}
</style>
<div class="s2g"><div>${s2n_ctrl}</div><div>${s2kap_ctrl}</div><div>${s2nmc_ctrl}</div><div>${s2conf_ctrl}</div></div>`
function runS2CLT(n, kappa, nMC, conf, seed) {
const rng = (function(a) { return function() { a |= 0; a = a + 0x6D2B79F5 | 0; var t = Math.imul(a ^ a >>> 15, 1 | a); t = t + Math.imul(t ^ t >>> 7, 61 | t) ^ t; return ((t ^ t >>> 14) >>> 0) / 4294967296; }; })(seed || 137);
function dot3(a,b) { return a[0]*b[0]+a[1]*b[1]+a[2]*b[2]; }
function norm3(a) { return Math.sqrt(dot3(a,a)); }
function clamp(x) { return Math.max(-1,Math.min(1,x)); }
// VMF on S² with center at the north pole: f(x) ∝ exp(κz).
function rvmf(N, kap) {
var ek = Math.exp(kap), emk = Math.exp(-kap), den = ek - emk, out = [];
for (var i = 0; i < N; i++) {
var z = Math.log(emk + rng() * den) / kap;
var phi = 2 * Math.PI * rng();
var rxy = Math.sqrt(Math.max(0, 1 - z*z));
out.push([rxy * Math.cos(phi), rxy * Math.sin(phi), z]);
}
return out;
}
function logAt(p,q) {
var c=clamp(dot3(p,q)), th=Math.acos(c), sn=Math.sin(th);
if (th<1e-12) return [0,0,0];
var sc=th/Math.max(1e-12,sn);
return [sc*(q[0]-c*p[0]),sc*(q[1]-c*p[1]),sc*(q[2]-c*p[2])];
}
function expVec(p,v) {
var r=norm3(v);
if (r<1e-12) return p.slice();
var cs=Math.cos(r), sn=Math.sin(r)/r;
return [cs*p[0]+sn*v[0],cs*p[1]+sn*v[1],cs*p[2]+sn*v[2]];
}
function objective(p,pts) {
var f=0;
for (var i=0;i<pts.length;i++) {
var th=Math.acos(clamp(dot3(p,pts[i])));
f += 0.5*th*th/pts.length;
}
return f;
}
// Karcher-flow solver with Armijo backtracking; a second start is used
// for the smallest, most diffuse samples where local minima are likelier.
function karcherFrom(pts,start) {
var p=start.slice();
var np=norm3(p); p=p.map(x=>x/np);
for (var it=0;it<80;it++) {
var v=[0,0,0];
for (var i=0;i<pts.length;i++) {
var l=logAt(p,pts[i]);
v[0]+=l[0]/pts.length; v[1]+=l[1]/pts.length; v[2]+=l[2]/pts.length;
}
var nv=norm3(v);
if (nv<1e-9) break;
var oldF=objective(p,pts), step=1, next=p;
while (step>1e-6) {
next=expVec(p,v.map(x=>step*x));
if (objective(next,pts)<=oldF-0.25*step*nv*nv) break;
step*=0.5;
}
p=next;
if (step<=1e-6) break;
}
return p;
}
function intrinsicMeanS2(pts) {
var s=[0,0,0];
for (var i=0;i<pts.length;i++) {
s[0]+=pts[i][0]; s[1]+=pts[i][1]; s[2]+=pts[i][2];
}
var ns=norm3(s), extr=ns>1e-12?s.map(x=>x/ns):[0,0,1];
var p1=karcherFrom(pts,extr);
if (kappa>=2 || pts.length>=40) return p1;
var p2=karcherFrom(pts,[0,0,1]);
return objective(p1,pts)<=objective(p2,pts)?p1:p2;
}
function tbasis(p) {
var px = p[0], py = p[1], pz = p[2], r = Math.sqrt(px*px + py*py);
if (r > 1e-12) return [[py/r, -px/r, 0], [-px*pz/r, -py*pz/r, r]];
var sg = pz > 0 ? 1 : -1;
return [[1, 0, 0], [0, -sg, 0]];
}
function expS2(p, tb, a, b) {
return expVec(p,[
a*tb[0][0]+b*tb[1][0],
a*tb[0][1]+b*tb[1][1],
a*tb[0][2]+b*tb[1][2]
]);
}
function eig2(a,b,d) {
var tr=a+d, disc=Math.sqrt(Math.max(0,(a-d)*(a-d)+4*b*b));
var l1=Math.max(1e-10,(tr+disc)/2), l2=Math.max(1e-10,(tr-disc)/2);
var e1=[1,0];
if (Math.abs(b)>1e-12 || Math.abs(l1-a)>1e-12) {
e1=[b,l1-a];
var ne=Math.hypot(e1[0],e1[1]); e1=[e1[0]/ne,e1[1]/ne];
} else if (d>a) e1=[0,1];
return {l1:l1,l2:l2,e1:e1,e2:[-e1[1],e1[0]]};
}
function sandwichAtMean(p,pts) {
var tb=tbasis(p), logs=[], h11=0,h12=0,h22=0;
for (var i=0;i<pts.length;i++) {
var lv=logAt(p,pts[i]), x=dot3(lv,tb[0]), y=dot3(lv,tb[1]);
logs.push([x,y]);
var r=Math.hypot(x,y);
if (r<1e-10) { h11+=1; h22+=1; continue; }
var rcot=r/Math.tan(r), ux=x/r, uy=y/r;
h11+=rcot+(1-rcot)*ux*ux;
h12+=(1-rcot)*ux*uy;
h22+=rcot+(1-rcot)*uy*uy;
}
h11/=pts.length; h12/=pts.length; h22/=pts.length;
var mx=logs.reduce((s,z)=>s+z[0],0)/logs.length;
var my=logs.reduce((s,z)=>s+z[1],0)/logs.length;
var c11=0,c12=0,c22=0;
for (var i=0;i<logs.length;i++) {
var x=logs[i][0]-mx,y=logs[i][1]-my;
c11+=x*x/(logs.length-1); c12+=x*y/(logs.length-1); c22+=y*y/(logs.length-1);
}
var detH=h11*h22-h12*h12;
if (detH<1e-8) { h11+=1e-4; h22+=1e-4; detH=h11*h22-h12*h12; }
var i11=h22/detH,i12=-h12/detH,i22=h11/detH;
var a11=i11*c11+i12*c12, a12=i11*c12+i12*c22;
var a21=i12*c11+i22*c12, a22=i12*c12+i22*c22;
var s11=a11*i11+a12*i12;
var s12=a11*i12+a12*i22;
var s22=a21*i12+a22*i22;
return {tb:tb,s11:s11,s12:s12,s22:s22};
}
var mu0 = [0, 0, 1], scaled = [], lastPts = null, lastMuH = null;
for (var mc = 0; mc < nMC; mc++) {
var pts = rvmf(n, kappa);
var mh = intrinsicMeanS2(pts);
var lg = logAt(mu0,mh);
scaled.push([Math.sqrt(n)*lg[0],Math.sqrt(n)*lg[1]]);
if (mc === nMC - 1) { lastPts = pts; lastMuH = mh; }
}
var mu1 = 0, mu2 = 0;
for (var i = 0; i < scaled.length; i++) { mu1 += scaled[i][0]/nMC; mu2 += scaled[i][1]/nMC; }
var s11 = 0, s12 = 0, s22 = 0;
for (var i = 0; i < scaled.length; i++) {
var d1 = scaled[i][0] - mu1, d2 = scaled[i][1] - mu2;
s11 += d1*d1/(nMC-1); s12 += d1*d2/(nMC-1); s22 += d2*d2/(nMC-1);
}
var mcEig=eig2(s11,s12,s22);
var lam1=mcEig.l1,lam2=mcEig.l2,e1=mcEig.e1,e2=mcEig.e2;
var chi2c = -2 * Math.log(1 - conf);
var nEll = 120, ellU = [];
for (var i = 0; i < nEll; i++) {
var t = 2 * Math.PI * i / nEll;
ellU.push([
Math.sqrt(lam1*chi2c) * e1[0]*Math.cos(t) + Math.sqrt(lam2*chi2c) * e2[0]*Math.sin(t),
Math.sqrt(lam1*chi2c) * e1[1]*Math.cos(t) + Math.sqrt(lam2*chi2c) * e2[1]*Math.sin(t)
]);
}
// A genuine one-sample plug-in sandwich estimate for the displayed sample.
var plug=sandwichAtMean(lastMuH,lastPts), plugEig=eig2(plug.s11,plug.s12,plug.s22);
var tb=plug.tb,scCr=Math.sqrt(chi2c/n),confS2=[];
for (var i = 0; i <= nEll; i++) {
var t = 2 * Math.PI * i / nEll;
var a=scCr*(Math.sqrt(plugEig.l1)*plugEig.e1[0]*Math.cos(t)+Math.sqrt(plugEig.l2)*plugEig.e2[0]*Math.sin(t));
var b=scCr*(Math.sqrt(plugEig.l1)*plugEig.e1[1]*Math.cos(t)+Math.sqrt(plugEig.l2)*plugEig.e2[1]*Math.sin(t));
confS2.push(expS2(lastMuH, tb, a, b));
}
var lv0=logAt(lastMuH,mu0),aC=dot3(lv0,tb[0]),bC=dot3(lv0,tb[1]);
var detS=plug.s11*plug.s22-plug.s12*plug.s12;
var qi=(aC*aC*plug.s22-2*aC*bC*plug.s12+bC*bC*plug.s11)/detS;
var covered=n*qi<=chi2c;
return {
n:n, kappa:kappa, nMC:nMC, conf:conf, scaled:scaled,
s11:s11, s12:s12, s22:s22, lam1:lam1, lam2:lam2, e1:e1, e2:e2, chi2c:chi2c,
lastPts:lastPts, lastMuH:lastMuH, mu0:mu0, ellU:ellU, confS2:confS2, covered:covered
};
}
s2R = runS2CLT(s2n, s2kap, s2nmc, s2conf, 137);
html`
<div style="font-family:system-ui,sans-serif;max-width:960px;">
${s2ctrls}
<div style="display:flex;gap:18px;flex-wrap:wrap;margin-bottom:10px;">
${(() => {
var R = s2R;
var ca = Math.cos(22*Math.PI/180), sa = Math.sin(22*Math.PI/180);
function rot(p) { return [p[0], ca*p[1]-sa*p[2], sa*p[1]+ca*p[2]]; }
function proj(p) { return [p[0], p[1]]; }
var scx = 210, scy = 230, sRad = 170;
// Latitude circles + meridians
var latAng = [75,60,45,30,15], gridP = '';
for (var li = 0; li < latAng.length; li++) {
var phi = latAng[li]*Math.PI/180, rho = Math.sin(phi), z0 = Math.cos(phi), pts = '';
for (var i = 0; i <= 80; i++) {
var th = 2*Math.PI*i/80, p = rot([rho*Math.cos(th), rho*Math.sin(th), z0]), pp = proj(p);
pts += (i===0?'M':'L')+' '+(scx+sRad*pp[0]).toFixed(1)+' '+(scy-sRad*pp[1]).toFixed(1);
}
gridP += '<path d="'+pts+'" fill="none" stroke="#dee2e6" stroke-width="0.5"/>';
}
for (var j = 0; j < 8; j++) {
var th0 = j*Math.PI/4, pts = '';
for (var i = 0; i <= 80; i++) {
var phi2 = Math.PI*i/80;
var p = rot([Math.sin(phi2)*Math.cos(th0), Math.sin(phi2)*Math.sin(th0), Math.cos(phi2)]), pp = proj(p);
pts += (i===0?'M':'L')+' '+(scx+sRad*pp[0]).toFixed(1)+' '+(scy-sRad*pp[1]).toFixed(1);
}
gridP += '<path d="'+pts+'" fill="none" stroke="#dee2e6" stroke-width="0.5"/>';
}
// Data points
var dots = '';
for (var i = 0; i < R.lastPts.length; i++) {
var rp = rot(R.lastPts[i]);
if (rp[2] < -0.05) continue;
var pp = proj(rp), op = Math.max(0.12, Math.min(1, rp[2]+0.5));
dots += '<circle cx="'+(scx+sRad*pp[0]).toFixed(1)+'" cy="'+(scy-sRad*pp[1]).toFixed(1)+'" r="1.8" fill="#868e96" opacity="'+op.toFixed(2)+'"/>';
}
// Confidence region on sphere
var cPath = '';
for (var i = 0; i < R.confS2.length; i++) {
var rp = rot(R.confS2[i]), pp = proj(rp);
cPath += (i===0?'M':'L')+' '+(scx+sRad*pp[0]).toFixed(1)+' '+(scy-sRad*pp[1]).toFixed(1);
}
var mu0r = rot(R.mu0), mu0p = proj(mu0r);
var muHr = rot(R.lastMuH), muHp = proj(muHr);
// Tangent space layout
var tml=58, tmr=45, tmt=42, tmb=38, tpw=420-tml-tmr, tph=400-tmt-tmb;
var u1a=[], u2a=[];
for (var i=0;i<R.scaled.length;i++){u1a.push(R.scaled[i][0]);u2a.push(R.scaled[i][1]);}
var elim = Math.max(
Math.abs(Math.min.apply(null,u1a)), Math.abs(Math.max.apply(null,u1a)),
Math.abs(Math.min.apply(null,u2a)), Math.abs(Math.max.apply(null,u2a)),
Math.sqrt(R.lam1*R.chi2c)*1.2, Math.sqrt(R.lam2*R.chi2c)*1.2
);
var u1min=-elim, u1max=elim, u2min=-elim, u2max=elim;
function ux(v){return tml+(v-u1min)/(u1max-u1min)*tpw;}
function uy(v){return tmt+tph-(v-u2min)/(u2max-u2min)*tph;}
// Scatter
var td = '';
for (var i=0;i<R.scaled.length;i++){
td += '<circle cx="'+ux(R.scaled[i][0]).toFixed(1)+'" cy="'+uy(R.scaled[i][1]).toFixed(1)+'" r="1.5" fill="#74c0fc" opacity="0.55"/>';
}
// Confidence ellipse
var ep = '';
for (var i=0;i<R.ellU.length;i++){
ep += (i===0?'M':'L')+' '+ux(R.ellU[i][0]).toFixed(1)+' '+uy(R.ellU[i][1]).toFixed(1);
}
// Marginal histograms
var nB=18, h1=Array(nB).fill(0), h2=Array(nB).fill(0);
for (var i=0;i<R.scaled.length;i++){
var bi1=Math.min(nB-1,Math.max(0,Math.floor((R.scaled[i][0]-u1min)/(u1max-u1min)*nB)));
var bi2=Math.min(nB-1,Math.max(0,Math.floor((R.scaled[i][1]-u2min)/(u2max-u2min)*nB)));
if(bi1>=0&&bi1<nB)h1[bi1]++; if(bi2>=0&&bi2<nB)h2[bi2]++;
}
var hm1=Math.max.apply(null,h1.concat([1])), hm2=Math.max.apply(null,h2.concat([1]));
var hh=40, hw=40, h1r='', h2r='';
for (var i=0;i<nB;i++){
var bh1=h1[i]/hm1*hh, bh2=h2[i]/hm2*hw;
h1r+='<rect x="'+ux(u1min+(i/nB)*(u1max-u1min)).toFixed(1)+'" y="'+(tmt-bh1).toFixed(1)+'" width="'+(tpw/nB-0.5).toFixed(1)+'" height="'+bh1.toFixed(1)+'" fill="#74c0fc" opacity="0.45"/>';
h2r+='<rect x="'+(tml+tpw).toFixed(1)+'" y="'+uy(u2min+((i+1)/nB)*(u2max-u2min)).toFixed(1)+'" width="'+bh2.toFixed(1)+'" height="'+(tph/nB-0.5).toFixed(1)+'" fill="#74c0fc" opacity="0.45"/>';
}
return `
<div style="flex:1;min-width:400px;padding:10px;background:#f8f9fa;border-radius:6px;">
<h4 style="margin:0 0 4px 0;text-align:center;">S²: Data & Confidence Region</h4>
<svg width="420" height="440" viewBox="0 0 420 440" style="display:block;margin:0 auto;">
<defs><clipPath id="s2c"><circle cx="${scx}" cy="${scy}" r="${sRad}"/></clipPath></defs>
<circle cx="${scx}" cy="${scy}" r="${sRad}" fill="#f1f3f5" stroke="#adb5bd" stroke-width="1.2"/>
<g clip-path="url(#s2c)">
${gridP}
<path d="${cPath} Z" fill="#fa5252" opacity="0.22" stroke="#e03131" stroke-width="1.6"/>
${dots}
</g>
<circle cx="${(scx+sRad*mu0p[0]).toFixed(1)}" cy="${(scy-sRad*mu0p[1]).toFixed(1)}" r="5" fill="#1971c2" stroke="white" stroke-width="1.5"/>
<text x="${(scx+sRad*mu0p[0]+8).toFixed(1)}" y="${(scy-sRad*mu0p[1]+4).toFixed(1)}" font-size="11" fill="#1971c2" font-weight="bold">μ₀</text>
<circle cx="${(scx+sRad*muHp[0]).toFixed(1)}" cy="${(scy-sRad*muHp[1]).toFixed(1)}" r="5" fill="#e03131" stroke="white" stroke-width="1.5"/>
<text x="${(scx+sRad*muHp[0]+8).toFixed(1)}" y="${(scy-sRad*muHp[1]+4).toFixed(1)}" font-size="11" fill="#e03131" font-weight="bold">μ̂</text>
</svg>
<p style="font-size:0.78em;color:#495057;text-align:center;margin:5px 0 0 0;">
Gray dots: one sample (n=${R.n}). Blue: true mean μ₀. Red: estimated mean μ̂.<br>
Shaded: ${(R.conf*100).toFixed(0)}% confidence region on S² (exp<sub>μ̂</sub> of CLT ellipsoid).
<b>${R.covered ? '✓ μ₀ inside region' : '✗ μ₀ outside region'}</b>
${R.covered ? '' : '(expected ≈'+(R.conf*100).toFixed(0)+'% of trials)'}
</p>
</div>
<div style="flex:1;min-width:400px;padding:10px;background:#f8f9fa;border-radius:6px;">
<h4 style="margin:0 0 4px 0;text-align:center;">Tangent Space T<sub>μ₀</sub>S² ≅ ℝ²</h4>
<svg width="420" height="400" viewBox="0 0 420 400" style="display:block;margin:0 auto;">
<line x1="${tml}" y1="${tmt+tph}" x2="${tml+tpw}" y2="${tmt+tph}" stroke="#adb5bd" stroke-width="1"/>
<line x1="${tml}" y1="${tmt}" x2="${tml}" y2="${tmt+tph}" stroke="#adb5bd" stroke-width="1"/>
<line x1="${ux(0).toFixed(1)}" y1="${tmt}" x2="${ux(0).toFixed(1)}" y2="${tmt+tph}" stroke="#dee2e6" stroke-width="0.8" stroke-dasharray="4,3"/>
<line x1="${tml}" y1="${uy(0).toFixed(1)}" x2="${tml+tpw}" y2="${uy(0).toFixed(1)}" stroke="#dee2e6" stroke-width="0.8" stroke-dasharray="4,3"/>
${h1r}${h2r}
${td}
<path d="${ep} Z" fill="none" stroke="#e03131" stroke-width="2.2" opacity="0.85"/>
<circle cx="${ux(0).toFixed(1)}" cy="${uy(0).toFixed(1)}" r="3.5" fill="#1971c2"/>
<text x="${(tml+tpw/2).toFixed(1)}" y="${tmt+tph+26}" text-anchor="middle" font-size="10">u₁ = √n · log_μ₀(μ̂)₁</text>
<text x="14" y="${(tmt+tph/2).toFixed(1)}" text-anchor="middle" font-size="10" transform="rotate(-90,14,${(tmt+tph/2).toFixed(1)})">u₂</text>
<text x="${(tml+tpw-5).toFixed(1)}" y="${(tmt-10).toFixed(1)}" text-anchor="end" font-size="9" fill="#868e96">marginal (u₁)</text>
<text x="${(tml+tpw+18).toFixed(1)}" y="${(tmt+12).toFixed(1)}" text-anchor="start" font-size="9" fill="#868e96">marg. (u₂)</text>
</svg>
<p style="font-size:0.78em;color:#495057;text-align:center;margin:5px 0 0 0;">
nMC = ${R.nMC} replications. Blue dot: origin (true μ₀ in tangent coords).<br>
Red ellipse: fitted Gaussian ${(R.conf*100).toFixed(0)}% contour {u : uᵀΣ̂<sub>MC</sub>⁻¹u ≤ χ²₂,₁₋α = ${R.chi2c.toFixed(2)}}.
Marginal histograms help assess the Gaussian approximation coordinate by coordinate.
</p>
</div>`;
})()}
</div>
<div style="padding:10px 14px;background:#e7f5ff;border-radius:6px;border-left:4px solid #1971c2;font-size:0.84em;color:#212529;">
<b>How to read this:</b> The <b>left panel</b> shows one sample (gray dots) drawn from a von Mises–Fisher distribution on S² concentrated at the north pole μ₀ = (0,0,1)ᵀ. The concentration κ controls spread: larger κ → tighter clustering. The red point is the <b>intrinsic</b> sample mean, computed by Karcher flow. Its shaded confidence patch uses a plug-in estimate of the spherical Hessian sandwich from that one sample, then maps the $n^{-1/2}$-scaled tangent ellipse through exp<sub>μ̂</sub>. The <b>right panel</b> shows √n·log<sub>μ₀</sub>(μ̂) across ${s2R.nMC} Monte Carlo replications; its red contour uses the empirical Monte Carlo covariance only as a diagnostic fitted-Gaussian contour. As $n$ grows, the scaled distribution should stabilize toward its Gaussian limit while the unscaled confidence patch shrinks at rate $n^{-1/2}$.
</div>
</div>`
```
## Application: From Metric Space to Manifold — Inference on $\mathcal{S}_{++}^m$ {#sec-l17-application}
### Motivation: What the Manifold Adds
Lectures 4–10 developed a rich set of tools for data in metric spaces: Fréchet ANOVA (@DubeyMuller2019), global and local-linear Fréchet regression (@PetersenMueller2019), and metric-space logistic classification. In those lectures, the SPD cone $\mathcal{S}_{++}^m$ was treated as a **metric space** — we could compute distances, Fréchet means, and perform permutation tests, but we could not exploit differential structure.
Lectures 15–17 have now equipped us with the **Riemannian manifold** structure: tangent spaces $T_\Sigma\mathcal{S}_{++}^m \cong \operatorname{Sym}(m)$, exponential and logarithmic maps, parallel transport, and — crucially — the **Bhattacharya–Patrangenaru CLT** for intrinsic Fréchet means. This section revisits the SPD inference problems from earlier lectures and asks: **what additional capabilities does the manifold structure unlock?**
::: {.callout-note title="Recall: the portfolio risk context"}
The running example through Lectures 4–10 was **portfolio risk management with SPD covariance matrices**. Concretely:
- Each observation $\Sigma_i \in \mathcal{S}_{++}^m$ is a covariance matrix estimated from a rolling window of $m$ asset returns over (say) 60 trading days.
- The **Fréchet mean** $\hat{\mu}_n$ represents the "average" risk structure over the observation period.
- In [Lecture 4](lecture-04.qmd), we tested whether the risk structure shifted between calm and stress market regimes using **permutation-based Fréchet ANOVA**—a metric-space method requiring only distances.
- In [Lectures 6–9](lecture-06.qmd), we modeled how the risk profile evolves as a function of macroeconomic covariates via **global and local-linear Fréchet regression**—again, purely metric-space tools.
The same SPD framework underpins **diffusion tensor imaging (DTI)** (each voxel is a $3 \times 3$ SPD tensor describing water diffusion anisotropy), **brain functional connectivity** (correlation matrices across $m$ brain regions), and **computer vision** (region covariance descriptors).
In all of these applications, we treated $\mathcal{S}_{++}^m$ as a metric space. Now, with the manifold structure in hand, we can ask sharper questions: *How precise is our estimate of the average risk structure? In which directions of the tangent space is the uncertainty largest? Can we construct a confidence ellipsoid whose axes correspond to interpretable patterns of covariance change (e.g., overall scale, relative variances, or cross-asset correlations)?*
:::
The manifold CLT answers these questions. Concretely, it unlocks the following capabilities unavailable in the pure metric-space framework:
1. **Asymptotic confidence regions** — the CLT gives explicit Gaussian limits in the tangent space, enabling elliptical confidence sets for Fréchet means.
2. **Covariance-standardized tests** — Wald-type tests can gain efficiency when the tangent-space CLT and covariance model are accurate.
3. **Tangent-space PCA** — dimension reduction that respects the manifold geometry at the mean.
### The SPD Manifold
Recall from Lecture 15 that $\mathcal{S}_{++}^m$ admits two principal Riemannian metrics:
::: {.callout-tip title="Two metrics on the SPD manifold"}
**Affine-invariant (AI) metric:**
$$
g^{\mathrm{AI}}_\Sigma(U, V) = \operatorname{tr}(\Sigma^{-1}U\,\Sigma^{-1}V), \qquad
d_{\mathrm{AI}}(\Sigma_0, \Sigma_1) = \|\log(\Sigma_0^{-1/2}\Sigma_1\Sigma_0^{-1/2})\|_F.
$$
**Log-Euclidean (LE) metric:** pullback of the Frobenius inner product via the matrix logarithm,
$$
d_{\mathrm{LE}}(\Sigma_0, \Sigma_1) = \|\log\Sigma_0 - \log\Sigma_1\|_F.
$$
Under the LE metric, $(\mathcal{S}_{++}^m, g^{\mathrm{LE}})$ is **globally isometric** to the Euclidean space $(\operatorname{Sym}(m), \langle\cdot,\cdot\rangle_F)$—the manifold is flat [@ArsignyEtAl2007]. In contrast, the AI metric has nonpositive sectional curvature (zero in commuting directions and negative in noncommuting directions), not positive curvature [@PennecFillardAyache2006].
:::
**Why this matters for the CLT:** Under the log-Euclidean metric, the Fréchet mean is
$$
\hat{\mu}_n^{\mathrm{LE}}=\exp\!\left(\frac1n\sum_{i=1}^n\log\Sigma_i\right).
$$
Use the isometric vectorization
$$
\operatorname{svec}(A)=
\bigl(A_{ii}\ (1\le i\le m),\ \sqrt2A_{ij}\ (1\le i<j\le m)\bigr)^\top
\quad\text{in any fixed order},
\qquad
\|\operatorname{svec}(A)\|_2=\|A\|_F.
$$
For the half-squared-distance loss, the Hessian is the identity and the score covariance is $\operatorname{cov}\{\operatorname{svec}(\log X)\}$. Equivalently, for the squared-distance convention used in @thm-bp2005-clt, $\Lambda=2I$ and $C=4\operatorname{cov}\{\operatorname{svec}(\log X)\}$; the factors cancel. Hence
$$
\sqrt{n}\,\operatorname{svec}(\log\hat{\mu}_n^{\mathrm{LE}}-\log\mu)
\xrightarrow{\mathcal D}
N\!\left(0,\operatorname{cov}\{\operatorname{svec}(\log X)\}\right).
$$
In other words, the LE metric makes SPD inference **exactly Euclidean** in the log-domain. This is the bridge between the metric-space methods of Lectures 4–10 and the manifold CLT: under the LE metric, the manifold CLT is ordinary multivariate normality of log-covariance matrices.
### Confidence Regions for Fréchet Means on $\mathcal{S}_{++}^m$
Lecture 4 used permutation calibration for group comparisons, but a hypothesis test by itself does not provide a joint confidence region for one Fréchet mean. The CLT supplies such a region: an **ellipsoid in tangent coordinates**, mapped back to the manifold.
::: {#prp-spd-confidence .proposition title="Asymptotic confidence region on $\mathcal{S}_{++}^m$"}
Let $\Sigma_1, \ldots, \Sigma_n$ be i.i.d. random SPD matrices with Fréchet mean $\mu$ under either the AI or LE metric, and assume the corresponding CLT and nonsingularity conditions. Let $\hat{\mu}_n$ be the sample Fréchet mean and $\hat{\Lambda}_n$, $\hat{C}_n$ be consistent estimators represented in one orthonormal tangent basis. Let $\varphi_{\hat\mu_n}(\mu)$ denote the coordinate vector of $\log_{\hat\mu_n}(\mu)$ in that basis. An asymptotic $100(1-\alpha)\%$ confidence region is
$$
C_{1-\alpha} = \left\{ \mu : n \cdot \varphi_{\hat\mu_n}(\mu)^\top \hat{\Lambda}_n^\top \hat{C}_n^{-1} \hat{\Lambda}_n \varphi_{\hat\mu_n}(\mu) \le \chi^2_{d, 1-\alpha} \right\},
$$
where $d = m(m+1)/2$ is the manifold dimension. Equivalently, the region is the image under $\exp_{\hat{\mu}_n}$ of the ellipsoid $\{v \in T_{\hat{\mu}_n}\mathcal{S}_{++}^m : n \cdot v^\top \hat{\Sigma}^{-1} v \le \chi^2_{d,1-\alpha}\}$, where $\hat{\Sigma} = \hat{\Lambda}_n^{-1}\hat{C}_n(\hat{\Lambda}_n^\top)^{-1}$.
Under the LE metric, use $\varphi_{\hat\mu_n}(\mu)=\operatorname{svec}(\log\mu-\log\hat\mu_n)$. The region simplifies to the standard Gaussian ellipsoid in $\mathbb R^d$:
$$
C_{1-\alpha}^{\mathrm{LE}} = \left\{ \mu : n \cdot \|\operatorname{svec}(\log\hat{\mu}_n - \log\mu)\|^2_{\hat{\Sigma}^{-1}} \le \chi^2_{d, 1-\alpha} \right\},
$$
where $\hat{\Sigma} = \widehat{\operatorname{cov}}\{\operatorname{svec}(\log\Sigma_i)\}$.
:::
**What this buys us over the metric-space approach:**
- **Continuous confidence levels** (not limited by Monte Carlo resolution of randomization tests).
- **Explicit ellipsoid shape** — the eigenvectors of $\hat{\Sigma}$ reveal which directions in the tangent space are most uncertain.
- **Computationally cheap under LE geometry** — the mean and covariance are computed in the log domain, with no randomization loop.
### Two-Sample Tests: Metric vs. Manifold
Lecture 4 introduced the permutation-based Fréchet ANOVA for testing $H_0: \mu_1 = \mu_2$ against $H_1: \mu_1 \ne \mu_2$ in a general metric space. With the manifold CLT, we now have a **parametric alternative**: the Wald test in the tangent space.
::: {#def-wald-spd .definition title="Wald test for equality of two Fréchet means on a manifold"}
Let $\Sigma_1^{(1)}, \ldots, \Sigma_{n_1}^{(1)}$ and $\Sigma_1^{(2)}, \ldots, \Sigma_{n_2}^{(2)}$ be independent samples from two populations on $\mathcal{S}_{++}^m$ with Fréchet means $\mu_1, \mu_2$. Under the LE metric, let $Y_i^{(g)} = \operatorname{svec}(\log\Sigma_i^{(g)})$. If the two groups have a common nonsingular covariance, let $\hat{\Sigma}_p$ be the pooled estimator. Then
$$
W_n = \frac{n_1 n_2}{n_1 + n_2} (\bar{Y}^{(1)} - \bar{Y}^{(2)})^\top \hat{\Sigma}_p^{-1} (\bar{Y}^{(1)} - \bar{Y}^{(2)}),
$$
and under $H_0: \mu_1 = \mu_2$, we have $W_n \xrightarrow{\mathcal{D}} \chi^2_d$.
:::
::: {.callout-important title="Metric vs. manifold: a comparison"}
| Aspect | Permutation test (metric) | Wald test (manifold) |
|---|---|---|
| Validity | Exact under exchangeability | Asymptotic (CLT) |
| Computational cost | $B$ refits of Fréchet means | One eigendecomposition |
| Statistic | Unstandardized distance in the demo | Covariance-standardized distance |
| Confidence region | Not supplied by the test alone | Elliptical via the CLT |
| Requires manifold? | No (metric only) | Yes (tangent space + CLT) |
| Sensitivity to cut locus | No differentiability assumption | Requires CLT-specific cut-locus control |
:::
The key tradeoff: a label-permutation test is exact only under exchangeability under the null (and Monte Carlo permutation adds simulation error). A Wald test is asymptotic but can gain power when it appropriately standardizes an anisotropic covariance—for example, when the mean shift lies in a low-variance direction. Neither test uniformly dominates the other.
### Interactive Demo: Confidence Regions and Power Comparison on $\mathcal{S}_{++}^2$
The following demo simulates $2 \times 2$ SPD matrices under the log-Euclidean metric and compares:
1. **Coverage** of the CLT-based confidence ellipsoid for the first group mean.
2. **Power** of the Wald test vs. the permutation test for two-sample comparisons.
3. **Shape** of the confidence ellipsoid — which directions in $\operatorname{Sym}(2) \cong \mathbb{R}^3$ are most uncertain?
The log coordinates are exactly Gaussian in this simulation but anisotropic. The mean shift is placed in the relatively low-variance off-diagonal direction, a setting in which Mahalanobis standardization can help the Wald test. This is an illustration, not a claim that Wald tests uniformly dominate permutation tests.
```{ojs}
//| label: fig-spd-inference
//| fig-cap: "Interactive: Manifold inference on the SPD cone — confidence regions and power comparison"
spd_n1_control = Inputs.range([10, 100], {step: 5, value: 30, label: "Sample size (group 1)"})
spd_n2_control = Inputs.range([10, 100], {step: 5, value: 30, label: "Sample size (group 2)"})
spd_shift_control = Inputs.range([0, 1.5], {step: 0.05, value: 0.4, label: "Mean shift ‖log μ₁ − log μ₂‖_F"})
spd_concentration_control = Inputs.range([0.3, 3.0], {step: 0.1, value: 1.0, label: "Concentration (1/σ)"})
spd_nMC_control = Inputs.range([50, 300], {step: 50, value: 200, label: "MC replications"})
spd_n1 = Generators.input(spd_n1_control)
spd_n2 = Generators.input(spd_n2_control)
spd_shift = Generators.input(spd_shift_control)
spd_concentration = Generators.input(spd_concentration_control)
spd_nMC = Generators.input(spd_nMC_control)
spdi_controls_view = html`
<style>
.spdi-slider-grid { display:flex; flex-wrap:wrap; gap:6px 20px; width:100%; margin:0 auto 12px; font:0.85em system-ui,sans-serif; container-type:inline-size; }
.spdi-slider-grid > * { flex:1 1 calc((100% - 40px)/3); min-width:0; margin:0; }
.spdi-slider-grid input[type="number"] { width:7.5rem !important; }
@container (max-width:700px) { .spdi-slider-grid > * { flex-basis:calc((100% - 20px)/2); } }
@container (max-width:480px) { .spdi-slider-grid > * { flex-basis:100%; } }
</style>
<div class="spdi-slider-grid">
<div>${spd_n1_control}</div>
<div>${spd_n2_control}</div>
<div>${spd_shift_control}</div>
<div>${spd_concentration_control}</div>
<div>${spd_nMC_control}</div>
</div>`
function runSPDInference(n1, n2, shift, conc, nMC, seed) {
const rng = (function(a) { return function() { a |= 0; a = a + 0x6D2B79F5 | 0; var t = Math.imul(a ^ a >>> 15, 1 | a); t = t + Math.imul(t ^ t >>> 7, 61 | t) ^ t; return ((t ^ t >>> 14) >>> 0) / 4294967296; }; })(seed || 137);
// Standard normal via Box-Muller
function rnorm() {
const u1 = rng(), u2 = rng();
return Math.sqrt(-2 * Math.log(Math.max(u1, 1e-10))) * Math.cos(2 * Math.PI * u2);
}
// Isometric coordinates svec([[a,b],[b,c]]) = [a, √2 b, c].
const d = 3;
const sigma = 1.0 / conc;
// True log-mean of group 1: identity → log = 0
const mu1_log = [0, 0, 0];
// Unit-Frobenius shift in the relatively low-variance off-diagonal direction.
const shift_dir = [0, 1, 0];
const mu2_log = shift_dir.map(v => v * shift);
const waldCover = {count: 0, total: 0};
const waldReject = {count: 0, total: 0};
const permReject = {count: 0, total: 0};
let lastDetail = null;
function drawLogNoise() {
const z1=rnorm(), z2=rnorm(), z3=rnorm();
return [
sigma*(1.10*z1),
sigma*(0.10*z1+0.35*z2),
sigma*(0.25*z1+0.15*z2+0.80*z3)
];
}
for (let mc = 0; mc < nMC; mc++) {
// Generate group 1
const Y1 = [];
for (let i = 0; i < n1; i++) {
const eps=drawLogNoise();
const y = mu1_log.map((m, j) => m + eps[j]);
Y1.push(y);
}
// Generate group 2
const Y2 = [];
for (let i = 0; i < n2; i++) {
const eps=drawLogNoise();
const y = mu2_log.map((m, j) => m + eps[j]);
Y2.push(y);
}
// Sample means
const ybar1 = Array(d).fill(0);
const ybar2 = Array(d).fill(0);
for (let j = 0; j < d; j++) {
for (let i = 0; i < n1; i++) ybar1[j] += Y1[i][j] / n1;
for (let i = 0; i < n2; i++) ybar2[j] += Y2[i][j] / n2;
}
function sampleCov(Y,bar) {
const S=Array(d*d).fill(0);
for (let j = 0; j < d; j++) {
for (let k = 0; k < d; k++) {
let s=0;
for (let i=0;i<Y.length;i++) s+=(Y[i][j]-bar[j])*(Y[i][k]-bar[k]);
S[j*d+k]=s/(Y.length-1);
}
}
return S;
}
const S1=sampleCov(Y1,ybar1),S2=sampleCov(Y2,ybar2);
const Shat=Array(d*d).fill(0);
for (let j=0;j<d*d;j++) {
Shat[j]=((n1-1)*S1[j]+(n2-1)*S2[j])/(n1+n2-2);
}
// ---- Wald test ----
const diff = ybar1.map((v, j) => v - ybar2[j]);
// Pivoted Gaussian elimination for the 3×3 systems.
function solve3(S, b) {
const n = 3;
const A = S.slice();
const x = b.slice();
for (let col = 0; col < n; col++) {
let maxRow = col, maxVal = Math.abs(A[col*n + col]);
for (let row = col+1; row < n; row++) {
if (Math.abs(A[row*n + col]) > maxVal) {
maxVal = Math.abs(A[row*n + col]);
maxRow = row;
}
}
if (maxVal < 1e-12) return null;
if (maxRow !== col) {
for (let j = 0; j < n; j++) {
const tmp = A[col*n + j]; A[col*n + j] = A[maxRow*n + j]; A[maxRow*n + j] = tmp;
}
const tmp = x[col]; x[col] = x[maxRow]; x[maxRow] = tmp;
}
const piv = A[col*n + col];
for (let j = col; j < n; j++) A[col*n + j] /= piv;
x[col] /= piv;
for (let row = 0; row < n; row++) {
if (row === col) continue;
const f = A[row*n + col];
for (let j = col; j < n; j++) A[row*n + j] -= f * A[col*n + j];
x[row] -= f * x[col];
}
}
return x;
}
const Sinv_diff = solve3(Shat, diff);
let W = 0;
for (let j = 0; j < d; j++) W += diff[j] * Sinv_diff[j];
W *= (n1 * n2) / (n1 + n2);
const waldCrit = 7.8147; // chi2(3, 0.95)
// Coverage of the group-1 mean uses group 1's own covariance estimate.
waldCover.total++;
const err1=mu1_log.map((v,j)=>v-ybar1[j]);
const S1invErr=solve3(S1,err1);
let W_conf=0;
for (let j=0;j<d;j++) W_conf+=err1[j]*S1invErr[j];
W_conf*=n1;
if (W_conf<=waldCrit) waldCover.count++;
if (W > waldCrit) waldReject.count++;
waldReject.total++;
// ---- Monte Carlo permutation test (B=99; plus-one p-value) ----
const B = 99;
const allY = Y1.concat(Y2);
const nTotal = n1 + n2;
let Fn_obs = 0;
// Observed: distance between group means
for (let j = 0; j < d; j++) Fn_obs += (ybar1[j] - ybar2[j]) * (ybar1[j] - ybar2[j]);
Fn_obs = Math.sqrt(Fn_obs); // Euclidean distance in log-domain
let n_ge = 0;
for (let b = 0; b < B; b++) {
// Shuffle indices
const perm = [];
for (let i = 0; i < nTotal; i++) perm.push(i);
for (let i = nTotal - 1; i > 0; i--) {
const jj = Math.floor(rng() * (i + 1));
const tmp = perm[i]; perm[i] = perm[jj]; perm[jj] = tmp;
}
const m1 = Array(d).fill(0), m2 = Array(d).fill(0);
for (let i = 0; i < n1; i++) {
const pt = allY[perm[i]];
for (let j = 0; j < d; j++) m1[j] += pt[j] / n1;
}
for (let i = n1; i < nTotal; i++) {
const pt = allY[perm[i]];
for (let j = 0; j < d; j++) m2[j] += pt[j] / n2;
}
let Fn_b = 0;
for (let j = 0; j < d; j++) Fn_b += (m1[j] - m2[j]) * (m1[j] - m2[j]);
Fn_b = Math.sqrt(Fn_b);
if (Fn_b >= Fn_obs) n_ge++;
}
const p_perm = (1 + n_ge) / (B + 1);
if (p_perm < 0.05) permReject.count++;
permReject.total++;
// Save last MC for display
if (mc === nMC - 1) {
function eig3(S) {
const A=[[S[0],S[1],S[2]],[S[3],S[4],S[5]],[S[6],S[7],S[8]]];
for (let it=0;it<30;it++) {
let p=0,q=1,max=Math.abs(A[0][1]);
for (const pair of [[0,2],[1,2]]) {
const v=Math.abs(A[pair[0]][pair[1]]);
if (v>max) { max=v; p=pair[0]; q=pair[1]; }
}
if (max<1e-12) break;
const ang=0.5*Math.atan2(2*A[p][q],A[q][q]-A[p][p]);
const c=Math.cos(ang),s=Math.sin(ang);
const app=c*c*A[p][p]-2*s*c*A[p][q]+s*s*A[q][q];
const aqq=s*s*A[p][p]+2*s*c*A[p][q]+c*c*A[q][q];
for (let k=0;k<3;k++) if (k!==p&&k!==q) {
const akp=A[k][p],akq=A[k][q];
A[k][p]=A[p][k]=c*akp-s*akq;
A[k][q]=A[q][k]=s*akp+c*akq;
}
A[p][p]=app; A[q][q]=aqq; A[p][q]=A[q][p]=0;
}
return [A[0][0],A[1][1],A[2][2]].sort((x,y)=>y-x);
}
const eigVals=eig3(S1);
lastDetail = {
ybar1: ybar1, ybar2: ybar2, diff: diff,
W: W, p_perm: p_perm, waldCrit: waldCrit,
n1: n1, n2: n2, shift: shift, sigma: sigma,
Shat: Shat, eigVals: eigVals
};
}
}
return {
waldCoverage: waldCover.count / Math.max(1, waldCover.total),
waldPower: waldReject.count / Math.max(1, waldReject.total),
permPower: permReject.count / Math.max(1, permReject.total),
lastDetail: lastDetail,
nMC: nMC
};
}
spdRes = runSPDInference(spd_n1, spd_n2, spd_shift, spd_concentration, spd_nMC, 42);
html`
<div style="font-family: system-ui, sans-serif; max-width: 900px;">
${spdi_controls_view}
<div style="padding: 14px; background: #f8f9fa; border-radius: 8px; margin-bottom: 16px;">
<h4 style="margin-top: 0;">Inference Comparison: Wald (Manifold CLT) vs. Permutation (Metric Space)</h4>
<table style="width: 100%; border-collapse: collapse; font-size: 0.9em;">
<tr style="border-bottom: 2px solid #dee2e6;">
<th style="padding: 6px 10px; text-align: left;">Metric</th>
<th style="padding: 6px 10px; text-align: center;">Wald (manifold CLT)</th>
<th style="padding: 6px 10px; text-align: center;">Permutation (metric)</th>
<th style="padding: 6px 10px; text-align: left;">Notes</th>
</tr>
<tr style="border-bottom: 1px solid #e9ecef;">
<td style="padding: 6px 10px;"><b>Group-1 mean coverage (95% target)</b></td>
<td style="padding: 6px 10px; text-align: center; font-weight: bold; color: ${Math.abs((spdRes.waldCoverage||0) - 0.95) < 1.96*Math.sqrt(0.95*0.05/spd_nMC) ? '#2b8a3e' : '#c92a2a'};">${((spdRes.waldCoverage||0)*100).toFixed(1)}%</td>
<td style="padding: 6px 10px; text-align: center;">—</td>
<td style="padding: 6px 10px; font-size: 0.83em; color: #868e96;">Asymptotic χ²₃ ellipsoid</td>
</tr>
<tr>
<td style="padding: 6px 10px;"><b>Power (α = 0.05)</b></td>
<td style="padding: 6px 10px; text-align: center; font-weight: bold; font-size: 1.1em; color: ${spdRes.waldPower > spdRes.permPower ? '#2b8a3e' : '#c92a2a'};">${(spdRes.waldPower*100).toFixed(1)}%</td>
<td style="padding: 6px 10px; text-align: center; font-size: 1.1em; color: ${spdRes.permPower > spdRes.waldPower ? '#2b8a3e' : '#c92a2a'};">${(spdRes.permPower*100).toFixed(1)}%</td>
<td style="padding: 6px 10px; font-size: 0.83em; color: #495057;">Observed Monte Carlo difference = ${(spdRes.waldPower-spdRes.permPower).toFixed(3)}</td>
</tr>
</table>
</div>
${spdRes.lastDetail ? (() => {
const ld = spdRes.lastDetail;
const labels = ['(log Σ)₁₁', '√2 (log Σ)₁₂', '(log Σ)₂₂'];
return `
<div style="display: flex; gap: 16px; flex-wrap: wrap; margin-bottom: 16px;">
<div style="flex: 1; min-width: 280px; padding: 12px; background: #f1f3f5; border-radius: 6px;">
<h4 style="margin-top: 0;">Estimated Log-Means (last replication)</h4>
<table style="width: 100%; border-collapse: collapse; font-size: 0.85em;">
<tr style="border-bottom: 1px solid #dee2e6;">
<th style="padding: 3px 6px;">Dim</th>
<th style="padding: 3px 6px; text-align: right;">μ̂₁ (log)</th>
<th style="padding: 3px 6px; text-align: right;">μ̂₂ (log)</th>
<th style="padding: 3px 6px; text-align: right;">Diff</th>
</tr>
${labels.map((l, j) => `
<tr>
<td style="padding: 3px 6px; font-family: monospace;">${l}</td>
<td style="padding: 3px 6px; text-align: right;">${ld.ybar1[j].toFixed(3)}</td>
<td style="padding: 3px 6px; text-align: right;">${ld.ybar2[j].toFixed(3)}</td>
<td style="padding: 3px 6px; text-align: right; font-weight: bold;">${ld.diff[j].toFixed(3)}</td>
</tr>`).join('')}
</table>
<p style="margin: 8px 0 4px 0; font-size: 0.85em;">
<b>Wald statistic:</b> W = ${ld.W.toFixed(2)}
${ld.W > ld.waldCrit ? '<span style="color:#c92a2a;"> > χ²₃,₀.₉₅ = 7.81 → reject H₀</span>' : '<span style="color:#2b8a3e;"> < 7.81 → do not reject H₀</span>'}
<br>
<b>Permutation p-value (B=99):</b> ${ld.p_perm.toFixed(2)}
${ld.p_perm < 0.05 ? '<span style="color:#c92a2a;"> → reject H₀</span>' : '<span style="color:#2b8a3e;"> → do not reject H₀</span>'}
</p>
</div>
<div style="flex: 1; min-width: 280px; padding: 12px; background: #f1f3f5; border-radius: 6px;">
<h4 style="margin-top: 0;">Confidence Ellipsoid Shape</h4>
<p style="font-size: 0.85em; color: #495057;">
Eigenvalues of Σ̂ (variance along each tangent direction):
</p>
<table style="width: 100%; border-collapse: collapse; font-size: 0.85em;">
<tr style="border-bottom: 1px solid #dee2e6;">
<th style="padding: 3px 6px;">Axis</th>
<th style="padding: 3px 6px; text-align: right;">Variance</th>
<th style="padding: 3px 6px; text-align: right;">Std Dev</th>
</tr>
${ld.eigVals.map((v, j) => `
<tr>
<td style="padding: 3px 6px;">Axis ${j+1}</td>
<td style="padding: 3px 6px; text-align: right; font-family: monospace;">${v.toFixed(4)}</td>
<td style="padding: 3px 6px; text-align: right; font-family: monospace;">${Math.sqrt(Math.abs(v)).toFixed(4)}</td>
</tr>`).join('')}
</table>
<p style="font-size: 0.82em; color: #868e96; margin-top: 8px;">
The ellipsoid is widest along Axis 1 (largest eigenvalue) — this direction has the most uncertainty.
For 2×2 SPD matrices, these axes correspond to linear combinations of the log-variance and log-covariance entries.
</p>
</div>
</div>`;
})() : ''}
<div style="padding: 12px; background: #e7f5ff; border-radius: 6px; border-left: 4px solid #1971c2;">
<b>How to read this:</b>
<ul style="margin: 4px 0 0 0; padding-left: 18px; font-size: 0.85em; color: #495057;">
<li>The coverage row is computed for group 1 at every shift value; it should approach 95% as n₁ grows.</li>
<li>Increase <b>shift</b> to compare empirical rejection probabilities. Here the shift lies in a low-variance direction, so covariance standardization favors the Wald statistic.</li>
<li>Decrease <b>concentration</b> (increase σ) to reduce the signal-to-noise ratio. The log data remain exactly Gaussian in this particular simulation.</li>
<li>The permutation test uses the unstandardized log-Euclidean distance and B=99 random permutations per replication; its p-values therefore have resolution 0.01.</li>
<li>The ellipsoid eigenvalues quantify anisotropic uncertainty in isometric svec coordinates.</li>
</ul>
</div>
</div>
`
```
::: {.callout-tip title="Key observations from the demo"}
- The Wald test weights a direction by its estimated inverse variance; the distance-based permutation statistic weights all isometric coordinates equally. Their power can therefore differ under anisotropy.
- At shift zero, the permutation test has finite-sample validity from exchangeability, whereas the displayed Wald cutoff is asymptotic.
- Larger samples improve covariance estimation and make the $\chi^2_3$ approximation more accurate.
- The ellipsoid axes reveal which linear combinations of log-variance and log-covariance entries are estimated most precisely. A scalar distance alone does not encode this directional uncertainty.
:::
### Discussion: When Does the Manifold Structure Help?
The manifold perspective on $\mathcal{S}_{++}^m$ — and on Riemannian manifolds more generally — adds value over the pure metric-space perspective in several ways:
**1. Tangent-space coordinates enable classical multivariate statistics.** Under the log-Euclidean metric, the SPD manifold is isometric to $\mathbb{R}^{m(m+1)/2}$. This means every Euclidean multivariate tool — PCA, linear discriminant analysis, MANOVA, Gaussian graphical models — can be applied **rigorously** in the log-domain, with the assurance that distances and means are preserved by the isometry.
**2. The CLT provides explicit covariance estimates.** The sandwich formula $\hat{\Lambda}^{-1}\hat{C}(\hat{\Lambda}^\top)^{-1}$ yields confidence regions and standardized effect sizes in tangent coordinates. A permutation test can be accompanied by an effect estimate, but its scalar test statistic alone does not provide a covariance ellipsoid.
**3. Geodesic regression and PGA.** Beyond the Fréchet mean, smooth Riemannian structure supplies tangent vectors, exponential maps, and parallel transport for **geodesic regression** and **principal geodesic analysis** (PGA; @FletcherJoshi2004). General geodesic metric spaces may support related constructions, but not this differential toolkit.
**4. Dimension reduction with geometric meaning.** Tangent-space PCA at the Fréchet mean identifies the directions of maximum variance **on the manifold**. For SPD matrices, these directions correspond to interpretable patterns of covariance change — e.g., an eigenvector aligned with the off-diagonal entry indicates that the correlation between two assets is the primary source of variability.
**5. When the manifold structure does NOT help.** With small samples or a poor Gaussian approximation, an exchangeability-based permutation test may be more reliable than an asymptotic Wald cutoff. Cut-locus complications matter on manifolds such as spheres; the AI and LE SPD geometries themselves are complete and have no cut locus. If only a group-comparison p-value is needed, a distance-based test may also be simpler.
::: {.callout-important title="Take-home message"}
The manifold CLT is **complementary** to metric-space inference, not a replacement. Use the manifold approach when:
- The CLT approximation and covariance estimate are empirically credible.
- You need **continuous confidence regions** or **effect sizes with directions**.
- You want to do **dimension reduction** (PGA, geodesic PCA).
Stick with metric-space permutation/bootstrap methods when:
- Sample sizes are small.
- Exchangeability is plausible but the tangent-space approximation is suspect.
- Only a hypothesis test decision (reject/retain) is needed.
- Computational budget is generous and asymptotic approximations are suspect.
:::
## Key Takeaways
- The **BP 2005 CLT** gives $\sqrt{n}$-asymptotic normality with sandwich covariance $\Lambda^{-1}C(\Lambda^\top)^{-1}$ under small-ball and cut-locus-avoidance conditions.
- In the **mean-centered chart**, the limiting covariance is $4\Lambda^{-1}\operatorname{cov}(Z)\Lambda^{-1}$, where $Z = \log_\mu(X)$.
- **Kendall–Le (2011)** extend the CLT to non-i.i.d. data with a Lindeberg condition, using parallel transport to linearize around $\mu$.
- The **manifold CLT enables** asymptotic confidence regions, Wald tests, and tangent-space PCA that are unavailable in a pure metric-space framework.
- The **Wald test** and **permutation test** are complementary: covariance standardization can help Wald in anisotropic settings, while permutation calibration is exact under exchangeability.
## Exercises
1. **Euclidean cancellation.** For squared Euclidean distance, compute $\Lambda$ and $C$ in the BP convention and show that the sandwich covariance reduces to $\operatorname{cov}(X)$. <a href="javascript:void(0)" onclick="showSolution('l17-sol-1')" class="solution-link">📝 Show Solution</a>
2. **BP conditions on $S^2$.** Suppose a distribution on the unit sphere is supported in a closed geodesic ball of radius $r$ about $x_0$. What radius bound in @thm-bp2005-clt guarantees the classical BP CLT? Why does a von Mises–Fisher law with finite $\kappa$ not satisfy this theorem? <a href="javascript:void(0)" onclick="showSolution('l17-sol-2')" class="solution-link">📝 Show Solution</a>
3. **Confidence region on $\mathcal{S}_{++}^2$.** Given $n = 50$ i.i.d. $2 \times 2$ SPD matrices with sample Fréchet mean $\hat{\mu}_n$ and estimated log-covariance $\hat{\Sigma}$, construct the 95% asymptotic confidence ellipsoid for $\mu$ under the LE metric. Interpret the axes of the ellipsoid. <a href="javascript:void(0)" onclick="showSolution('l17-sol-5')" class="solution-link">📝 Show Solution</a>
4. **Wald vs. permutation power.** From the interactive demo in Section @sec-l17-application, explain why covariance standardization can make the Wald test more sensitive to a shift in a low-variance direction. When would you prefer permutation calibration? <a href="javascript:void(0)" onclick="showSolution('l17-sol-6')" class="solution-link">📝 Show Solution</a>
<style>
.solution-link { font-size: 0.9em; text-decoration: none; white-space: nowrap; margin-left: 0.3em; }
.solution-link:hover { text-decoration: underline; }
.solution-dialog { padding: 0; max-width: 720px; }
.solution-dialog-header { display: flex; justify-content: space-between; align-items: flex-start; border-bottom: 1px solid #dee2e6; padding: 1.25rem 1.5rem 1rem; background: #f8f9fa; border-radius: 8px 8px 0 0; }
.solution-dialog-header h4 { margin: 0; font-size: 1.15rem; }
.solution-dialog-close { background: none; border: 1px solid #adb5bd; border-radius: 4px; padding: 0.2rem 0.75rem; cursor: pointer; font-size: 0.9rem; color: #495057; white-space: nowrap; flex-shrink: 0; }
.solution-dialog-close:hover { background: #e9ecef; }
.solution-original { padding: 1rem 1.5rem; background: #f1f3f5; border-left: 4px solid #868e96; margin: 1rem 1.5rem; border-radius: 4px; font-size: 0.95rem; }
.solution-answer { padding: 0.5rem 1.5rem 1.5rem; }
.solution-answer strong { color: #2b8a3e; }
dialog { border: none; border-radius: 8px; box-shadow: 0 8px 32px rgba(0,0,0,0.22); padding: 0; max-width: 750px; width: 90vw; }
dialog::backdrop { background: rgba(0,0,0,0.45); }
</style>
<dialog id="l17-sol-1"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 1</h4><button onclick="closeSolution('l17-sol-1')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Verify the sandwich covariance for squared Euclidean distance.</div><div class="solution-answer"><strong>Solution:</strong> In Euclidean coordinates, $h(\theta,X)=\|\theta-X\|^2$, so $\nabla_\theta h(\mu,X)=2(\mu-X)$ and $\nabla_\theta^2h(\mu,X)=2I$. Thus $\Lambda=2I$ and $C=4\operatorname{cov}(X)$. Therefore $\Lambda^{-1}C(\Lambda^\top)^{-1}=(\frac12I)\,4\operatorname{cov}(X)\,(\frac12I)=\operatorname{cov}(X)$, exactly the covariance in the ordinary sample-mean CLT.</div></div></dialog>
<dialog id="l17-sol-2"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 2</h4><button onclick="closeSolution('l17-sol-2')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Apply the BP support condition on $S^2$.</div><div class="solution-answer"><strong>Solution:</strong> The unit $2$-sphere has sectional curvature $1$, so $K=1$ and the theorem requires $r<\pi/4$. Such a ball is also disjoint from the antipodal cut locus of its center. A von Mises–Fisher density with finite concentration is strictly positive on all of $S^2$, so its support is the entire sphere and cannot fit inside any such ball. High concentration makes violations rare but does not change the support; the classical small-ball theorem still does not apply.</div></div></dialog>
<dialog id="l17-sol-3"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 3</h4><button onclick="closeSolution('l17-sol-3')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Why can a degenerate quadratic term produce a slower rate?</div><div class="solution-answer"><strong>Solution:</strong> In one dimension, suppose the population criterion is locally $F(t)-F(0)\asymp t^r$ for an even $r>2$, while the empirical gradient contributes noise of order $n^{-1/2}$. Balancing $t^{r-1}$ against $n^{-1/2}$ gives $t=O_p\{n^{-1/(2r-2)}\}$, slower than $n^{-1/2}$. For example, a quartic leading term ($r=4$) gives rate $n^{-1/6}$. The limit is obtained by a nonlinear transformation of Gaussian score noise and is generally non-Gaussian. A singular Hessian signals that the quadratic argument fails, but higher-order growth and regularity are also needed; singularity by itself does not determine the rate or uniqueness.</div></div></dialog>
<dialog id="l17-sol-4"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 4</h4><button onclick="closeSolution('l17-sol-4')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Role of parallel transport in the Kendall-Le proof.</div><div class="solution-answer"><strong>Solution:</strong> The gradient condition at $\hat{\mu}_n$ is $\sum \log_{\hat{\mu}_n}(X_i) = 0$ — these are vectors in $T_{\hat{\mu}_n}\mathcal{M}$. To linearize around $\mu$, we need to compare vectors in $T_{\hat{\mu}_n}\mathcal{M}$ to vectors in $T_\mu\mathcal{M}$. **Parallel transport** $\Pi_{\hat{\mu}_n, \mu}$ moves the log vectors to the common tangent space $T_\mu\mathcal{M}$ while preserving inner products. This gives $\sum \Pi_{\hat{\mu}_n, \mu} \log_{\hat{\mu}_n}(X_i) = 0$ in $T_\mu\mathcal{M}$. Then we expand $\Pi_{\gamma(t), \mu} \log_{\gamma(t)}(X_i)$ along the geodesic from $\mu$ to $\hat{\mu}_n$, using the Hessian identity. Without parallel transport, there would be no way to add vectors from different tangent spaces.</div></div></dialog>
<dialog id="l17-sol-5"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 5</h4><button onclick="closeSolution('l17-sol-5')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Construct a 95% confidence ellipsoid on $\mathcal{S}_{++}^2$.</div><div class="solution-answer"><strong>Solution:</strong> **Construction:** Under the LE metric, let $Y_i = \operatorname{svec}(\log\Sigma_i) \in \mathbb{R}^3$, where $\operatorname{svec}([[a,b],[b,c]])=(a,\sqrt2b,c)^\top$. The sample Fréchet mean is $\hat{\mu}_n^{\mathrm{LE}} = \exp(\frac{1}{n}\sum_i \log\Sigma_i)$ with log-coordinates $\bar{Y} = \frac{1}{n}\sum_i Y_i$. Estimate $\hat{\Sigma} = \frac{1}{n-1}\sum_i (Y_i - \bar{Y})(Y_i - \bar{Y})^\top$. The 95% confidence ellipsoid is $\{y \in \mathbb{R}^3 : n(y - \bar{Y})^\top \hat{\Sigma}^{-1} (y - \bar{Y}) \le \chi^2_{3,0.95} = 7.815\}$. Map it back via $\mu = \exp\{\operatorname{svec}^{-1}(y)\}$.
**Interpretation of axes:** In the eigendecomposition $\hat{\Sigma}=V\Lambda V^\top$, eigenvector $v_j$ is a linear combination of the two log-diagonal directions and the scaled log-off-diagonal direction. The corresponding semi-axis length is $\sqrt{\chi^2_{3,0.95}\lambda_j/n}$. Thus the largest eigenvalue identifies the least precisely estimated pattern of log-covariance change; no particular entry is necessarily the largest-variance direction.</div></div></dialog>
<dialog id="l17-sol-6"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 6</h4><button onclick="closeSolution('l17-sol-6')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Compare Wald and permutation calibration.</div><div class="solution-answer"><strong>Solution:</strong> Write the mean difference in covariance-eigenvector coordinates. The Wald statistic weights its component in direction $j$ by $1/\hat\lambda_j$, whereas the unstandardized log-Euclidean distance gives every isometric coordinate the same weight. A shift along a reliably estimated low-variance direction can therefore be much more visible to Wald. This is not uniform dominance: estimating a nearly singular covariance can be unstable, and the $\chi^2$ cutoff is only asymptotic. Prefer permutation calibration when exchangeability holds and finite-sample type-I error is the priority, especially when the covariance estimate or Gaussian approximation is doubtful. Report an effect estimate and uncertainty separately regardless of which calibration is used.</div></div></dialog>
<script>
function showSolution(id) { const d = document.getElementById(id); if (d) { d.showModal(); if (window.MathJax && MathJax.typesetPromise) MathJax.typesetPromise([d]).catch(function(e) { console.log('MathJax error:', e); }); } }
function closeSolution(id) { const d = document.getElementById(id); if (d) d.close(); }
document.addEventListener('click', function(e) { if (e.target.tagName === 'DIALOG') e.target.close(); });
</script>
## Further Reading
- @BhattacharyaPatrangenaru2005 — Classical CLT on manifolds.
- @KendallLe2011 — Non-i.i.d. CLT with Lindeberg condition.
- @HotzHuckemann2015Circle — Uniqueness, computation, and antipode-sensitive CLTs for intrinsic means on $S^1$.
- @HotzLeWood2024 — CLT with cut locus in the support.
- @DubeyMuller2019 — Fréchet ANOVA and two-sample permutation tests in metric spaces (Lecture 4).
- @PetersenMueller2019 — Fréchet regression for metric-space valued responses (Lectures 6–8).
- @PennecFillardAyache2006 — Fréchet means and tangent-space PCA on the SPD manifold.
- @FletcherJoshi2004 — Principal geodesic analysis (PGA) on manifolds.
- @ArsignyEtAl2007 — Log-Euclidean metrics for SPD matrices and the LE isometry to $\mathbb{R}^d$.
## Self-Assessment Quiz
Test your understanding of this lecture with the interactive MCQ quiz:
👉 **[Lecture 17 Quiz — 10 Multiple-Choice Questions](../quizzes/lecture-17-quiz.qmd)**