---
title: "Lecture 16: Riemannian Manifolds — Fréchet Means"
subtitle: "First-order conditions, existence, and uniqueness 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 and derive the first-order condition $\mathbb{E}\{\log_\mu(X)\} = 0$ for a Fréchet mean on a Riemannian manifold.
- Explain the role of the cut locus in the differentiability of the Fréchet function.
- State the existence theorem of Bhattacharya and Patrangenaru (2003).
- Contrast uniqueness on Hadamard manifolds (global) versus positively curved manifolds (local, small-ball).
- Apply Afsari's theorem to determine conditions for unique Fréchet means.
- Explain why empirical Fréchet means on finite-dimensional Hadamard manifolds attain the rate $O_P(n^{-1/2})$ under a finite-second-moment condition.
- Apply Riemannian gradient descent (RGD) to compute Fréchet means on the SPD manifold and compare with the metric-space approach.
## Gradient of Squared Distance
For a smooth $f : \mathcal{M} \to \mathbb{R}$, the **gradient** $\nabla f(p) \in T_p\mathcal{M}$ satisfies $df_p(v) = g_p(\nabla f(p), v)$ for all $v$.
::: {#prp-gradient-squared .proposition title="Gradient of squared distance"}
Fix $y \in \mathcal{M}$. If $y \notin \operatorname{Cut}(x)$, then $f_y(x) = d^2(x, y)$ is smooth near $x$ and
$$
\nabla f_y(x) = -2\log_x(y).
$$
:::
This formula is important for manifold statistics: it connects the geometry (log map) to optimization (gradient of the Fréchet function).
## First-Order Characterization
On a Riemannian manifold, the Fréchet function gains differentiability away from the cut locus. The explicit gradient formula is the key tool.
::: {#prp-frechet-first-order .proposition title="First-order condition for a population Fréchet mean"}
Let $(\mathcal{M}, g)$ be a complete Riemannian manifold. Let $X$ be a random element with $\mathbb{E}\{d^2(o, X)\} < \infty$, and $F(x) = \mathbb{E}\{d^2(x, X)\}$. Suppose $\mu$ is a population Fréchet mean and $\mathbb{P}\{X \in \operatorname{Cut}(\mu)\} = 0$. Then
$$
\mathbb{E}\{\log_\mu(X)\} = 0.
$$
:::
**Proof sketch.** For $v \in T_\mu\mathcal{M}$, let $\gamma(t) = \exp_\mu(tv)$. By the first-variation formula in the above, for $y \notin \operatorname{Cut}(\mu)$,
$$
\left.\frac{d}{dt}d^2(\gamma(t),y)\right|_{t=0}
=-2\langle \log_\mu(y),v\rangle_\mu.
$$
The local difference quotients are dominated by an integrable multiple of $1+d(\mu,X)$, so differentiation may pass through the expectation. Hence
$$
dF_\mu(v)=-2\left\langle \mathbb{E}\{\log_\mu(X)\},v\right\rangle_\mu.
$$
Because a differentiable function has zero differential at a minimizer, this expression vanishes for every $v$, and nondegeneracy of the metric gives the result.
::: {.corollary title="Sample first-order condition"}
If $\mu$ minimizes $F_n(x) = \frac{1}{n}\sum d^2(x, x_i)$ and $x_i \notin \operatorname{Cut}(\mu)$ for all $i$, then $\sum_{i=1}^n \log_\mu(x_i) = 0$.
:::
::: {.callout-tip title="Generating data with a given Fréchet mean"}
Let $Z$ be a random element of an injectivity domain in $T_\mu\mathcal{M}$, with $\mathbb{E}Z=0$, and set $X=\exp_\mu(Z)$. Then $\log_\mu(X)=Z$, so $\mu$ is a **stationary point** of the Fréchet function. This alone does not prove that $\mu$ is a mean: a stationary point can be a maximum or a nonglobal local minimum. If the support lies in an Afsari ball, or if $\mathcal M$ is Hadamard, the relevant uniqueness and convexity results make the stationary point the unique Fréchet mean.
:::
## Existence and Uniqueness
In [Lecture 3](lecture-03.qmd) we proved the fundamental existence result of @BhattacharyaPatrangenaru2003: on a metric space where every closed bounded set is compact, the Fréchet mean set is nonempty and compact whenever the Fréchet function is finite at one point.
::: {#cor-hopf-rinow-frechet .corollary title="Riemannian manifolds satisfy the BP2003 condition"}
By the Hopf--Rinow theorem, a complete connected finite-dimensional Riemannian manifold has the property that every closed bounded set is compact. Hence the existence theorem of @BhattacharyaPatrangenaru2003 applies: if $\int d^2(p,x)\,dQ(x) < \infty$ for some $p$, then the Fréchet mean set is nonempty and compact.
:::
For uniqueness, recall from [Lecture 2](lecture-02.qmd) that on a Hadamard manifold, the $\operatorname{CAT}(0)$ inequality makes $d^2(\cdot,y)$ strongly geodesically convex [@Sturm2003]. After integration, the Fréchet function remains strongly geodesically convex, so the Fréchet mean is unique whenever the second moment is finite.
On positively curved manifolds like spheres, uniqueness requires the data to be sufficiently concentrated.
::: {#thm-afsari .theorem title="Local uniqueness (Afsari 2011, Theorem 2.1)"}
Let $\mathcal{M}$ be complete with sectional curvatures bounded above by $\kappa$ and injectivity radius $\operatorname{inj}(\mathcal{M})$. Define
$$
\rho_\kappa=\frac12\min\!\left\{\operatorname{inj}(\mathcal M),\frac{\pi}{\sqrt{\kappa}}\right\},
$$
where $\pi/\sqrt{\kappa}=\infty$ if $\kappa\le 0$. If $Q$ is supported in a ball $B(o,\rho)$ with $\rho<\rho_\kappa$, then its squared-distance Fréchet mean exists uniquely, lies in $B(o,\rho)$, and is the only stationary point of $F$ **inside that ball** [@Afsari2011].
:::
::: {.callout-note title="The standard remedy"}
Work inside a ball small relative to the curvature bound and injectivity radius. Afsari's theorem then localizes the unique global minimizer and rules out other stationary points in that ball. It does not assert that every squared-distance term is globally convex throughout an arbitrary positively curved ball.
:::
## Convergence Rate on Hadamard Manifolds
::: {#cor-hadamard-rate .corollary title="Parametric rate in finite dimensions"}
Let $X$ be a random element on an $m$-dimensional Hadamard manifold with $\mathbb{E} d^2(X, x) < \infty$. Then
$$
d(\hat{\mu}_n, \mu) = O_P(n^{-1/2}).
$$
:::
**Proof idea.** The $\operatorname{CAT}(0)$ variance inequality gives quadratic growth,
$$
F(q)-F(\mu)\ge d^2(q,\mu).
$$
Thus the growth exponent is $2$. Locally, an $m$-dimensional Riemannian manifold has polynomial covering numbers, so the square-root log entropy grows more slowly than every power $(\delta/\varepsilon)^\alpha$ with $\alpha>0$. Taking $\alpha<1$ in the general empirical-process rate theorem yields $d(\hat\mu_n,\mu)=O_P(n^{-1/2})$ [@Schotz2019]. The finite-dimensional assumption is doing real work here; the same conclusion is not automatic in an arbitrary infinite-dimensional Hadamard space.
## Interactive Exploration: Fréchet Mean on a Spherical Cap
On $\mathbb{S}^2$, $\kappa=1$ and $\operatorname{inj}(\mathbb S^2)=\pi$, so Afsari's sufficient radius is $\pi/2$. For a finite sample, containment in a north-centered cap of radius strictly below $\pi/2$ guarantees a unique Fréchet mean. The demo can also move beyond that range: the guarantee then disappears, but nonuniqueness does **not** follow automatically. The optimizer uses backtracking and reports whether it actually reached a stationary point.
**Visual guide:**
- **Colored dots** = data points on $\mathbb{S}^2$, color-coded by distance from the cap center
- **Red marker** = estimated sample Fréchet mean $\hat{\mu}$
- **Green marker** = cap center $\mu^*$ (the population mean in the guaranteed range)
- **Dashed circle** = hemisphere boundary ($z = 0$, i.e., distance $\pi/2$ from the north pole)
- **Status indicators** distinguish the sufficient Afsari check from numerical convergence
```{ojs}
//| label: fig-frechet-mean-s2
//| fig-cap: "Interactive: Sample Fréchet mean for spherical-cap data on 𝕊²"
hm_n_control = Inputs.range([5, 150], {step: 1, value: 25, label: "Sample size n"})
hm_spread_control = Inputs.range([0.1, 1.8], {step: 0.02, value: 0.66, label: "Spread (cap radius, rad)"})
hm_seed_control = Inputs.range([1, 30], {step: 1, value: 3, label: "Random seed"})
hm_n = Generators.input(hm_n_control)
hm_spread = Generators.input(hm_spread_control)
hm_seed = Generators.input(hm_seed_control)
hm_controls_view = html`
<style>
#fig-frechet-mean-s2 .quarto-subfloat-caption,
#fig-rgd-spd .quarto-subfloat-caption { display:none; }
.hm-slider-grid { display:flex; flex-wrap:wrap; gap:6px 20px; width:100%; margin:0 0 12px; font:0.85em system-ui,sans-serif; container-type:inline-size; }
.hm-slider-grid > * { flex:1 1 calc((100% - 40px)/3); min-width:0; margin:0; }
.hm-slider-grid input[type="number"] { width:7.5rem !important; }
@container (max-width:700px) { .hm-slider-grid > * { flex-basis:calc((100% - 20px)/2); } }
@container (max-width:480px) { .hm-slider-grid > * { flex-basis:100%; } }
</style>
<div class="hm-slider-grid">
<div>${hm_n_control}</div>
<div>${hm_spread_control}</div>
<div>${hm_seed_control}</div>
</div>`
function runFrechetS2(n, spread, seed) {
var 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);
// 3D helpers
function vdot(a, b) { return a[0]*b[0]+a[1]*b[1]+a[2]*b[2]; }
function vnorm(a) { return Math.sqrt(vdot(a,a)); }
function vscale(a, s) { return [a[0]*s, a[1]*s, a[2]*s]; }
function vadd(a, b) { return [a[0]+b[0], a[1]+b[1], a[2]+b[2]]; }
function vsub(a, b) { return [a[0]-b[0], a[1]-b[1], a[2]-b[2]]; }
function vnormalize(a) { var n = vnorm(a); return n < 1e-14 ? [0,0,1] : [a[0]/n, a[1]/n, a[2]/n]; }
function vcross(a, b) { return [a[1]*b[2]-a[2]*b[1], a[2]*b[0]-a[0]*b[2], a[0]*b[1]-a[1]*b[0]]; }
// North pole: cap center and population mean in the Afsari-guaranteed range.
var muTrue = [0, 0, 1];
// Generate data on S² within a spherical cap of radius `spread` around muTrue
var X = [];
for (var i = 0; i < n; i++) {
// Rejection sampling: generate uniform on sphere, keep if within cap
var pt;
var done = false;
while (!done) {
// Random point on S² (uniform)
var z = 1 - 2 * rng();
var r = Math.sqrt(Math.max(0, 1 - z*z));
var phi = rng() * 2 * Math.PI;
pt = [r*Math.cos(phi), r*Math.sin(phi), z];
// Check distance from north pole
var d = Math.acos(Math.max(-1, Math.min(1, vdot(pt, muTrue))));
if (d <= spread) done = true;
}
X.push(pt);
}
// ---- Fréchet mean via gradient descent ----
// On S²: log_p(q) = θ/sin(θ) * (q - cos(θ)*p) where θ = acos(<p,q>)
function logS2(p, q) {
var dot = vdot(p, q);
dot = Math.max(-1, Math.min(1, dot));
var theta = Math.acos(dot);
if (theta < 1e-12) return [0, 0, 0];
if (Math.PI - theta < 1e-10) return null; // multivalued at the antipode
var coeff = theta / Math.sin(theta);
var proj = vsub(q, vscale(p, dot));
return vscale(proj, coeff);
}
// exp_p(v) = cos(|v|)*p + sin(|v|)*v/|v|
function expS2(p, v) {
var nv = vnorm(v);
if (nv < 1e-12) return [p[0], p[1], p[2]];
var c = Math.cos(nv), s = Math.sin(nv);
return vadd(vscale(p, c), vscale(v, s/nv));
}
function objective(p) {
var total = 0;
for (var j = 0; j < n; j++) {
var dot = Math.max(-1, Math.min(1, vdot(p, X[j])));
var theta = Math.acos(dot);
total += theta * theta;
}
return total / n;
}
// Riemannian gradient descent with Armijo backtracking.
var muHat = [0, 0, 1]; // initialize at north pole
var maxIter = 150;
var gradNorms = [];
var objectiveValues = [];
var backtracks = 0;
var acceptedUpdates = 0;
var logDefined = true;
for (var iter = 0; iter < maxIter; iter++) {
// avgLog is half the negative gradient of the mean squared-distance objective.
var avgLog = [0, 0, 0];
for (var i = 0; i < n; i++) {
var li = logS2(muHat, X[i]);
if (li === null) {
logDefined = false;
break;
}
avgLog = vadd(avgLog, li);
}
if (!logDefined) break;
avgLog = vscale(avgLog, 1/n);
var gn = vnorm(avgLog);
gradNorms.push(gn);
var currentObjective = objective(muHat);
objectiveValues.push(currentObjective);
if (gn < 1e-8) break;
var stepSize = 1;
var accepted = false;
for (var bt = 0; bt < 20; bt++) {
var candidate = vnormalize(expS2(muHat, vscale(avgLog, stepSize)));
var candidateObjective = objective(candidate);
if (candidateObjective <= currentObjective - 2e-4 * stepSize * gn * gn) {
muHat = candidate;
accepted = true;
acceptedUpdates++;
backtracks += bt;
break;
}
stepSize *= 0.5;
}
if (!accepted) break;
}
// First-order residual at the returned iterate.
var sumLogs = [0, 0, 0];
for (var i = 0; i < n; i++) {
var finalLog = logS2(muHat, X[i]);
if (finalLog === null) {
logDefined = false;
break;
}
sumLogs = vadd(sumLogs, finalLog);
}
var sumLogNorm = logDefined ? vnorm(sumLogs) : NaN;
var avgLogNorm = logDefined ? sumLogNorm / n : NaN;
var converged = logDefined && avgLogNorm < 1e-8;
// Geodesic distances from the cap center
var dists = X.map(function(xi) { return Math.acos(Math.max(-1, Math.min(1, vdot(xi, muTrue)))); });
// Distance from muHat to muTrue
var muErr = Math.acos(Math.max(-1, Math.min(1, vdot(muHat, muTrue))));
// Check Afsari condition
var maxDist = Math.max.apply(null, dists);
var afsariOK = maxDist < Math.PI/2;
return {
X: X, muTrue: muTrue, muHat: muHat,
sumLogNorm: sumLogNorm, avgLogNorm: avgLogNorm, dists: dists,
muErr: muErr, maxDist: maxDist, afsariOK: afsariOK,
n: n, spread: spread, gradNorms: gradNorms,
objectiveValues: objectiveValues, backtracks: backtracks,
acceptedUpdates: acceptedUpdates,
logDefined: logDefined, converged: converged
};
}
hmRes = runFrechetS2(hm_n, hm_spread, hm_seed);
viewof hm_view = {
const root = html`
<div style="font-family: system-ui, sans-serif; max-width: 920px;">
<div style="display: flex; gap: 20px; flex-wrap: wrap;">
<div style="flex: 1 1 440px; min-width: 0;">
<h4>Fréchet Mean for a Spherical Cap on 𝕊²</h4>
<svg class="hm-globe" width="440" height="440" viewBox="0 0 440 440"
style="display:block; width:100%; max-width:440px; height:auto; border:1px solid #dee2e6; border-radius:4px; cursor:grab; touch-action:none; user-select:none; background:linear-gradient(180deg,#fbfdff 0%,#f5f8fb 100%);">
<defs>
<radialGradient id="hm-sphere-glow" cx="36%" cy="28%" r="70%">
<stop offset="0%" stop-color="#ffffff"/>
<stop offset="55%" stop-color="#e7f5ff"/>
<stop offset="100%" stop-color="#d8edf7"/>
</radialGradient>
</defs>
<g transform="translate(220, 220)">
<circle r="188" fill="url(#hm-sphere-glow)" stroke="#ced4da" stroke-width="1.2"/>
<g class="hm-scene"></g>
</g>
<g transform="translate(15, 402)">
<circle cx="0" cy="0" r="4.5" fill="#e03131" stroke="#fff" stroke-width="1.5"/><text x="9" y="4" font-size="10" fill="#495057">μ̂ (estimated)</text>
<circle cx="140" cy="0" r="4.5" fill="#2b8a3e" stroke="#fff" stroke-width="1.5"/><text x="149" y="4" font-size="10" fill="#495057">μ* (cap center)</text>
<line x1="280" y1="0" x2="295" y2="0" stroke="#f08c00" stroke-width="1.6" stroke-dasharray="6,3"/>
<text x="300" y="4" font-size="10" fill="#495057">Afsari boundary</text>
</g>
<text x="220" y="426" text-anchor="middle" font-size="11" fill="#6c757d">Drag to rotate the globe</text>
</svg>
</div>
<div style="flex: 1; min-width: 290px;">
<div style="padding: 12px; background: #f8f9fa; border-radius: 6px; margin-bottom: 12px;">
<h4 style="margin-top: 0;">Convergence & Diagnostics</h4>
<table style="width: 100%; border-collapse: collapse; font-size: 0.88em;">
<tr><td style="padding: 3px 8px;">Sample size n</td>
<td style="padding: 3px 8px; text-align: right;">${hmRes.n}</td></tr>
<tr><td style="padding: 3px 8px;">Cap radius (spread)</td>
<td style="padding: 3px 8px; text-align: right;">${(hmRes.spread*180/Math.PI).toFixed(1)}°</td></tr>
<tr><td style="padding: 3px 8px;">Max data distance from μ*</td>
<td style="padding: 3px 8px; text-align: right; font-family: monospace;">${(hmRes.maxDist*180/Math.PI).toFixed(2)}°</td></tr>
<tr><td style="padding: 3px 8px;">d(μ̂, μ*)</td>
<td style="padding: 3px 8px; text-align: right; font-family: monospace; font-weight: bold;">${(hmRes.muErr*180/Math.PI).toFixed(3)}°</td></tr>
<tr><td colspan="2"><hr style="margin: 4px 0;"></td></tr>
<tr><td style="padding: 3px 8px;">‖(1/n) Σ log<sub>μ̂</sub>(X<sub>i</sub>)‖</td>
<td style="padding: 3px 8px; text-align: right; font-family: monospace;">${hmRes.logDefined ? hmRes.avgLogNorm.toExponential(2) : 'undefined (cut locus)'}</td></tr>
<tr><td style="padding: 3px 8px;">‖∇F<sub>n</sub>(μ̂)‖</td>
<td style="padding: 3px 8px; text-align: right; font-family: monospace;">${hmRes.logDefined ? (2*hmRes.avgLogNorm).toExponential(2) : 'undefined'}</td></tr>
<tr><td style="padding: 3px 8px;">Accepted updates</td>
<td style="padding: 3px 8px; text-align: right;">${hmRes.acceptedUpdates}</td></tr>
</table>
</div>
<div style="padding: 12px; border-radius: 6px; ${hmRes.afsariOK ? 'background: #d3f9d8; border-left: 4px solid #2b8a3e;' : 'background: #ffe3e3; border-left: 4px solid #e03131;'}">
<b>Checked Afsari condition:</b> max distance from μ* = ${(hmRes.maxDist*180/Math.PI).toFixed(1)}°
${hmRes.afsariOK
? '<span style="color:#2b8a3e;">< π/2 (90°) → ✓ unique Fréchet mean guaranteed</span>'
: '<span style="color:#e03131;">≥ π/2 (90°) → no conclusion from this sufficient check</span>'}
</div>
<div style="margin-top: 10px; padding: 12px; border-radius: 6px; ${hmRes.converged ? 'background:#d3f9d8; border-left:4px solid #2b8a3e;' : 'background:#fff3cd; border-left:4px solid #f08c00;'}">
<b>Numerical status:</b>
${hmRes.converged
? '<span style="color:#2b8a3e;">✓ stationary residual below 10⁻⁸</span>'
: '<span style="color:#9c6b00;">not certified as stationary; inspect the residual</span>'}
</div>
</div>
</div>
</div>`;
const svg = root.querySelector(".hm-globe");
const scene = root.querySelector(".hm-scene");
const R = 185;
const res = hmRes;
root.value = {rotY: 20, rotX: 15};
function clamp(x, lo, hi) { return Math.max(lo, Math.min(hi, x)); }
function rotate(v) {
const ry = root.value.rotY * Math.PI / 180;
const rx = root.value.rotX * Math.PI / 180;
const x1 = v[0] * Math.cos(ry) + v[2] * Math.sin(ry);
const y1 = v[1];
const z1 = -v[0] * Math.sin(ry) + v[2] * Math.cos(ry);
return [x1, y1 * Math.cos(rx) - z1 * Math.sin(rx), y1 * Math.sin(rx) + z1 * Math.cos(rx)];
}
function proj(v) {
const rv = rotate(v);
return {x: rv[0] * R, y: -rv[1] * R, z: rv[2]};
}
// Draw only front-facing runs. This prevents lines on the far side of the
// sphere from appearing on top of the opaque globe.
function visiblePolyline(points, attrs) {
const runs = [];
let current = [];
for (let i = 0; i < points.length; i++) {
const q = proj(points[i]);
if (q.z >= 0) {
current.push(q);
} else if (current.length > 1) {
runs.push(current);
current = [];
} else {
current = [];
}
}
if (current.length > 1) runs.push(current);
return runs.map(function(run) {
const d = run.map(function(q, i) {
return (i === 0 ? "M" : "L") + " " + q.x.toFixed(1) + " " + q.y.toFixed(1);
}).join(" ");
return '<path d="' + d + '" fill="none" ' + attrs + '/>';
}).join("");
}
function latitude(z, n) {
const r = Math.sqrt(Math.max(0, 1 - z*z));
const pts = [];
for (let j = 0; j <= n; j++) {
const phi = j * 2 * Math.PI / n;
pts.push([r * Math.cos(phi), r * Math.sin(phi), z]);
}
return pts;
}
function meridian(phi, n) {
const pts = [];
for (let j = 0; j <= n; j++) {
const theta = j * Math.PI / 2 / n;
pts.push([Math.sin(theta) * Math.cos(phi), Math.sin(theta) * Math.sin(phi), Math.cos(theta)]);
}
return pts;
}
function render() {
const parts = [];
for (let j = 1; j <= 9; j++) {
const z = Math.cos(j * Math.PI / 20);
parts.push(visiblePolyline(latitude(z, 120), 'stroke="#74c0fc" stroke-width="0.75" opacity="0.42"'));
}
for (let m = 0; m < 24; m++) {
parts.push(visiblePolyline(meridian(m * Math.PI / 12, 80), 'stroke="#4dabf7" stroke-width="0.7" opacity="0.34"'));
}
for (let j = 1; j <= 5; j++) {
const z = j / 6;
parts.push(visiblePolyline(latitude(z, 120), 'stroke="#a5d8ff" stroke-width="2.2" opacity="0.18"'));
}
parts.push(visiblePolyline(latitude(0, 160), 'stroke="#f08c00" stroke-width="2.1" stroke-dasharray="8,4" opacity="0.9"'));
const eqLabelPt = proj([1, 0, 0]);
if (eqLabelPt.z >= 0) {
parts.push('<text x="' + (eqLabelPt.x + 6).toFixed(0) + '" y="' + (eqLabelPt.y - 7).toFixed(0) + '" font-size="10" fill="#f08c00" opacity="0.92">z = 0 (π/2)</text>');
}
const dataPts = [];
for (let i = 0; i < res.n; i++) {
const pt = proj(res.X[i]);
dataPts.push({x: pt.x, y: pt.y, z: pt.z, d: res.dists[i]});
}
dataPts.sort(function(a, b) { return a.z - b.z; });
for (let i = 0; i < dataPts.length; i++) {
const dp = dataPts[i];
if (dp.z < 0) continue;
const t = Math.min(1, dp.d / (Math.PI/2));
const rCol = Math.round(25 + 205*t);
const gCol = Math.round(113 - 60*t);
const bCol = Math.round(255 - 110*t);
const opacity = 0.72 + 0.25*dp.z;
parts.push('<circle cx="' + dp.x.toFixed(1) + '" cy="' + dp.y.toFixed(1) + '" r="4.2" fill="rgb(' + rCol + ',' + gCol + ',' + bCol + ')" stroke="rgba(255,255,255,0.75)" stroke-width="0.8" opacity="' + opacity.toFixed(2) + '"/>');
}
const muTP = proj(res.muTrue);
if (muTP.z >= 0) {
parts.push('<circle cx="' + muTP.x.toFixed(1) + '" cy="' + muTP.y.toFixed(1) + '" r="6.5" fill="#2b8a3e" stroke="#fff" stroke-width="2"/>');
parts.push('<text x="' + (muTP.x + 10).toFixed(0) + '" y="' + (muTP.y - 8).toFixed(0) + '" font-size="11" fill="#2b8a3e" font-weight="bold">μ*</text>');
}
const muHP = proj(res.muHat);
if (muHP.z >= 0) {
parts.push('<circle cx="' + muHP.x.toFixed(1) + '" cy="' + muHP.y.toFixed(1) + '" r="7" fill="#e03131" stroke="#fff" stroke-width="2.5"/>');
parts.push('<text x="' + (muHP.x + 10).toFixed(0) + '" y="' + (muHP.y - 8).toFixed(0) + '" font-size="11" fill="#e03131" font-weight="bold">μ̂</text>');
}
scene.innerHTML = parts.join("\\n");
}
let dragging = false;
let lastX = 0;
let lastY = 0;
svg.addEventListener("pointerdown", function(event) {
dragging = true;
lastX = event.clientX;
lastY = event.clientY;
svg.setPointerCapture(event.pointerId);
svg.style.cursor = "grabbing";
});
svg.addEventListener("pointermove", function(event) {
if (!dragging) return;
const dx = event.clientX - lastX;
const dy = event.clientY - lastY;
lastX = event.clientX;
lastY = event.clientY;
root.value = {
rotY: (root.value.rotY + dx * 0.55 + 360) % 360,
// Keep vertical rotation aligned with the pointer movement.
rotX: clamp(root.value.rotX + dy * 0.55, -80, 80)
};
render();
root.dispatchEvent(new CustomEvent("input", {bubbles: true}));
});
function stopDrag(event) {
if (!dragging) return;
dragging = false;
svg.style.cursor = "grab";
if (event.pointerId !== undefined && svg.hasPointerCapture(event.pointerId)) {
svg.releasePointerCapture(event.pointerId);
}
}
svg.addEventListener("pointerup", stopDrag);
svg.addEventListener("pointercancel", stopDrag);
svg.addEventListener("pointerleave", stopDrag);
render();
return root;
}
```
::: {.callout-note title="How to read this"}
- Data points are **color-coded** by geodesic distance from $\mu^*$ (blue = close, red = far).
- The globe shows the visible half of the upper-hemisphere grid; drag it to change the viewpoint.
- The <b style="color:#f08c00;">dashed orange curve</b> is the visible part of $z = 0$, at distance $\pi/2$ from $\mu^*$.
- The <b style="color:#e03131;">red dot</b> $\hat\mu$ is the returned sample mean estimate.
- The residual tests the necessary first-order condition numerically; it is not displayed as successful unless the tolerance is met.
:::
::: {.callout-tip title="Try these experiments"}
- **Spread < $\pi/2$ (≈ 1.57 rad):** The displayed north-centered cap satisfies Afsari's sufficient condition, so the sample mean is unique. This theorem guarantees the minimizer, not convergence of an arbitrary fixed-step algorithm; the demo therefore uses backtracking.
- **Spread > $\pi/2$:** The checked sufficient condition fails. This means “no guarantee,” not “multiple means.” Compare seeds and inspect the convergence residual rather than inferring nonuniqueness from the warning.
- **Increase $n$:** Across repeated seeds, the typical distance $d(\hat\mu,\mu^*)$ should decrease on the root-$n$ scale. A single realization need not decrease monotonically with $n$.
- **Rotate the view:** Drag the globe left/right to rotate and up/down to tilt. The equator (Afsari boundary) is clearest from a side view.
- **Check the first-order condition:** A small $\|(1/n)\sum_i\log_{\hat\mu}(X_i)\|$ is a numerical stationarity check. Outside a convexity region, stationarity alone would not certify a global minimum.
:::
## Application: Fréchet Means on the SPD Manifold — Optimization via Gradient Descent
### From Metric Space to Manifold
Throughout Lectures 2, 4, and 6–9, we treated the SPD cone $\mathcal{S}_{++}^m$ as a **metric space** — a set equipped with a distance (log-Euclidean or affine-invariant) — and computed Fréchet means, ran ANOVA, and fitted regression curves. The computations relied on geometric "tricks" that were, in hindsight, exploiting manifold structure:
- **Under the log-Euclidean metric:** the Fréchet mean is simply $\exp(\frac{1}{n}\sum_i \log \Sigma_i)$ — a closed form that works because $\log : (\mathcal{S}_{++}^m, d_{\mathrm{LE}}) \to (\operatorname{Sym}(m), \|\cdot\|_F)$ is a **global isometry** to a flat Euclidean space [@ArsignyEtAl2006].
- **Under the affine-invariant metric:** no general closed form exists; the Fréchet mean usually requires **iterative computation**. The Riemannian structure supplies explicit derivatives and geodesic update rules.
Viewing $\mathcal{S}_{++}^m$ as a **Riemannian manifold** changes the picture. We now have:
1. **A tangent space** $T_\Sigma\mathcal{S}_{++}^m \cong \operatorname{Sym}(m)$ at every point — a vector space where we can do linear algebra.
2. **Exponential and logarithmic maps** $\exp_\Sigma$, $\log_\Sigma$ — which move between the manifold and the tangent space.
3. **An intrinsic gradient** $\nabla F$ — which gives the direction of steepest descent *respecting the manifold geometry*.
4. **Geodesic convexity** — which, together with a suitable step-size or line-search rule, supports global convergence guarantees on Hadamard manifolds.
### Riemannian Gradient Descent for the Fréchet Mean
For the function $f_S(\Sigma)=d^2(\Sigma,S)$, differentiated with respect to $\Sigma$,
$$
\operatorname{grad}_\Sigma f_S=-2\log_\Sigma(S).
$$
For the empirical Fréchet function $F_n(\Sigma) = \frac{1}{n}\sum_{i=1}^n d^2(\Sigma, \Sigma_i)$, linearity gives
$$
\nabla F_n(\Sigma) = -\frac{2}{n}\sum_{i=1}^n \log_\Sigma(\Sigma_i).
$$
The **Riemannian gradient descent (RGD)** iteration is then
$$
\Sigma^{(k+1)} = \exp_{\Sigma^{(k)}}\!\big(-\eta_k \,\nabla F_n(\Sigma^{(k)})\big),
$$
where $\eta_k > 0$ is a step size and $\exp_{\Sigma^{(k)}}$ is the Riemannian exponential map. This is the natural generalization of Euclidean gradient descent: move in the direction of steepest descent *along a geodesic*.
### Two Metrics, Two Algorithms
The concrete form of RGD depends on the chosen Riemannian metric.
::: {.callout-note title="RGD under the affine-invariant metric"}
For the affine-invariant metric $g^{\mathrm{AI}}_\Sigma(U, V) = \operatorname{tr}(\Sigma^{-1}U\,\Sigma^{-1}V)$:
- **Logarithmic map:**
$$
\log^{\mathrm{AI}}_\Sigma(S) = \Sigma^{1/2}\,\log\!\big(\Sigma^{-1/2}S\,\Sigma^{-1/2}\big)\,\Sigma^{1/2},
$$
where $\Sigma^{1/2}$ is the symmetric matrix square root and $\log$ is the principal matrix logarithm.
- **Exponential map:**
$$
\exp^{\mathrm{AI}}_\Sigma(U) = \Sigma^{1/2}\,\exp\!\big(\Sigma^{-1/2}U\,\Sigma^{-1/2}\big)\,\Sigma^{1/2}.
$$
- **RGD step:**
$$
\Sigma^{(k+1)} = \Sigma_k^{1/2}\,\exp\!\Big(\frac{2\eta_k}{n}\sum_{i=1}^n \log\big(\Sigma_k^{-1/2}\Sigma_i\,\Sigma_k^{-1/2}\big)\Big)\,\Sigma_k^{1/2}.
$$
The sign is positive because $-\operatorname{grad}F_n=(2/n)\sum_i\log_\Sigma(\Sigma_i)$. Implementations use matrix factorizations or eigendecompositions for the square root, logarithm, and exponential.
**Key property:** The affine-invariant metric makes $(\mathcal{S}_{++}^m, g^{\mathrm{AI}})$ a **Hadamard manifold** — complete, simply connected, with nonpositive sectional curvature [@PennecFillardAyache2006]. Therefore:
- The Fréchet function $F_n$ is **strongly geodesically convex** and has a unique minimizer.
- RGD equipped with a valid line search converges toward that global minimizer. A fixed step such as $\eta_k=1/2$ is not a universal, data-independent convergence theorem.
- The Fréchet mean is globally unique.
:::
::: {.callout-note title="RGD under the log-Euclidean metric — collapses to Euclidean GD"}
Under the log-Euclidean metric $d_{\mathrm{LE}}(\Sigma_1, \Sigma_2) = \|\log\Sigma_1 - \log\Sigma_2\|_F$, the map $\log : \mathcal{S}_{++}^m \to \operatorname{Sym}(m)$ is a **global isometry**. The manifold is flat (sectional curvature $0$ everywhere). In this case, RGD reduces to **ordinary Euclidean gradient descent in the log-domain**:
1. Transform: $L_i = \log \Sigma_i$, $\;L^{(k)} = \log \Sigma^{(k)}$.
2. Compute Euclidean gradient: $G_k = -\frac{2}{n}\sum_i (L_i - L^{(k)})$.
3. Update in log-domain: $L^{(k+1)} = L^{(k)} - \eta_k G_k$.
4. Map back: $\Sigma^{(k+1)} = \exp(L^{(k+1)})$.
This converges in **one step** with $\eta = \frac{1}{2}$ to the closed-form solution $\hat{\Sigma}_{\mathrm{LE}} = \exp(\frac{1}{n}\sum_i \log \Sigma_i)$.
:::
### What the Manifold View Adds
The transition from "metric space with Fréchet means" to "Riemannian manifold with gradient descent" is not merely a change of language. It unlocks concrete capabilities:
<style>
.l16-capability-table { max-width: 100%; overflow-x: auto; }
</style>
::: {.l16-capability-table}
| Capability | Metric-Space View | Manifold View |
|:-----------|:------------------|:--------------|
| **Optimization algorithm** | General proximal, stochastic, or derivative-free methods | Adds systematic differential methods such as Riemannian gradient descent |
| **Gradient** | Not defined (no linear structure) | $\nabla F_n(\Sigma) = -\frac{2}{n}\sum_i \log_\Sigma(\Sigma_i)$ |
| **Step direction** | Depends on the available metric/geodesic structure | Negative Riemannian gradient in a tangent space |
| **Convergence guarantee** | Depends on metric convexity and the algorithm | Geodesic convexity + an admissible step rule can give global convergence |
| **Line search** | Can compare values along known metric curves | Adds differential slopes and backtracking along exponential-map curves |
| **Momentum / acceleration** | No natural analogue | Parallel transport → Riemannian Adam, Riemannian Nesterov |
| **Tangent-space linearization** | Not available | Reduce nonlinear manifold problems to linear tangent-space problems |
| **Second-order methods** | Not defined | Riemannian Newton, trust-region methods |
:::
**Riemannian gradient descent is a recurring computational template.** When the squared-distance objective is differentiable at the iterates and the needed maps are available, one uses
$$
\mu^{(k+1)} = \exp_{\mu^{(k)}}\!\Big(\frac{\eta_k}{n}\sum_{i=1}^n \log_{\mu^{(k)}}(X_i)\Big),
$$
with the factor $2$ absorbed into the step size. The geometry-specific $\exp$ and $\log$ maps change, and a line search or other safeguard may still be required. At cut loci the objective can be nonsmooth, so this formula is not a universal solver for every manifold and every dataset.
::: {.callout-tip title="From trick to principle"}
In earlier lectures, the log-Euclidean "trick" (log → average → exp) appeared to be a convenient algebraic manipulation. The manifold perspective reveals it as a **special case of Riemannian gradient descent** on a flat manifold, where one gradient step with $\eta=1/2$ solves the quadratic exactly. The affine-invariant case is genuinely curved and generally requires iterative, safeguarded updates.
:::
### Interactive Exploration: RGD on $2 \times 2$ SPD Matrices
The demo below constructs $2\times2$ SPD matrices whose **sample affine-invariant Fréchet mean is exactly** a known target $\Sigma^*$. It centers perturbations in the whitened tangent space, then maps them by the affine-invariant exponential. RGD uses Armijo backtracking, so every accepted update decreases the objective.
**Visual guide:**
- **Colored ellipses** = data SPD matrices (visualized as covariance ellipses), color-coded by distance from the exact affine-invariant sample mean
- **Thick red ellipse** = current RGD iterate $\Sigma^{(k)}$
- **Thick green ellipse** = exact sample affine-invariant mean $\Sigma^*$
- **Dashed blue ellipse** = log-Euclidean closed-form solution $\hat{\Sigma}_{\mathrm{LE}}$ (computed in one step)
- **Gradient norm plot** shows the observed convergence of $\|\nabla F_n(\Sigma^{(k)})\|$
- **Distance plot** tracks $d_{\mathrm{AI}}(\Sigma^{(k)}, \Sigma^*)$
```{ojs}
//| label: fig-rgd-spd
//| fig-cap: "Interactive: Riemannian gradient descent for the Fréchet mean on the SPD manifold (affine-invariant metric)"
spd_n_control = Inputs.range([5, 40], {step: 1, value: 15, label: "Sample size n"})
spd_spread_control = Inputs.range([0.1, 1.2], {step: 0.02, value: 0.56, label: "Tangent spread σ"})
spd_seed_control = Inputs.range([1, 30], {step: 1, value: 7, label: "Random seed"})
spd_lr_control = Inputs.range([0.1, 1.5], {step: 0.05, value: 0.6, label: "Initial step size η"})
spd_maxiter_control = Inputs.range([5, 160], {step: 1, value: 50, label: "Max iterations"})
spd_n = Generators.input(spd_n_control)
spd_spread = Generators.input(spd_spread_control)
spd_seed = Generators.input(spd_seed_control)
spd_lr = Generators.input(spd_lr_control)
spd_maxiter = Generators.input(spd_maxiter_control)
rgd_controls_view = html`
<style>
.rgd-slider-grid { display:flex; flex-wrap:wrap; gap:6px 20px; width:100%; margin:0 0 12px; font:0.85em system-ui,sans-serif; container-type:inline-size; }
.rgd-slider-grid > * { flex:1 1 calc((100% - 40px)/3); min-width:0; margin:0; }
.rgd-slider-grid input[type="number"] { width:7.5rem !important; }
@container (max-width:700px) { .rgd-slider-grid > * { flex-basis:calc((100% - 20px)/2); } }
@container (max-width:480px) { .rgd-slider-grid > * { flex-basis:100%; } }
</style>
<div class="rgd-slider-grid">
<div>${spd_n_control}</div>
<div>${spd_spread_control}</div>
<div>${spd_seed_control}</div>
<div>${spd_lr_control}</div>
<div>${spd_maxiter_control}</div>
</div>`
// A symmetric 2x2 matrix is stored as [a,b,c] = [[a,b],[b,c]].
// The angle formula handles diagonal matrices with either eigenvalue ordering.
function symEig2x2(M) {
var a = M[0], b = M[1], c = M[2];
var trace = a + c;
var disc = Math.sqrt(Math.max(0, (a-c)*(a-c) + 4*b*b));
var lam1 = (trace + disc) / 2;
var lam2 = (trace - disc) / 2;
var angle = 0.5 * Math.atan2(2*b, a-c);
var u1 = Math.cos(angle), v1 = Math.sin(angle);
var u2 = -v1, v2 = u1;
return {vals: [lam1, lam2], vecs: [[u1, v1], [u2, v2]]};
}
function symToFull(M) {
return [M[0], M[1], M[1], M[2]];
}
function fullMul(A, B) {
return [
A[0]*B[0] + A[1]*B[2], A[0]*B[1] + A[1]*B[3],
A[2]*B[0] + A[3]*B[2], A[2]*B[1] + A[3]*B[3]
];
}
function fullToSym(A) {
return [A[0], 0.5*(A[1]+A[2]), A[3]];
}
// A and B are symmetric; return A B A. Ordinary products of symmetric
// matrices need not be symmetric, so both full intermediate matrices matter.
function symCongruence(A, B) {
var Af = symToFull(A), Bf = symToFull(B);
return fullToSym(fullMul(fullMul(Af, Bf), Af));
}
function symFunc(M, f) {
var eig = symEig2x2(M);
var lam1 = eig.vals[0], lam2 = eig.vals[1];
var u1 = eig.vecs[0][0], v1 = eig.vecs[0][1];
var u2 = eig.vecs[1][0], v2 = eig.vecs[1][1];
var flam1 = f(lam1), flam2 = f(lam2);
// Q * diag(f(λ)) * Q^T
return [
flam1*u1*u1 + flam2*u2*u2,
flam1*u1*v1 + flam2*u2*v2,
flam1*v1*v1 + flam2*v2*v2
];
}
function symExp(M) { return symFunc(M, Math.exp); }
function symLog(M) { return symFunc(M, function(x) { return Math.log(Math.max(x, 1e-14)); }); }
function symSqrt(M) { return symFunc(M, function(x) { return Math.sqrt(Math.max(x, 0)); }); }
function symInvSqrt(M) { return symFunc(M, function(x) { return 1/Math.sqrt(Math.max(x, 1e-14)); }); }
function symFrob(M) { return Math.sqrt(M[0]*M[0] + 2*M[1]*M[1] + M[2]*M[2]); }
function distAI(S1, S2) {
var S1invSqrt = symInvSqrt(S1);
return symFrob(symLog(symCongruence(S1invSqrt, S2)));
}
function logAI(Sigma, S) {
var SigSqrt = symSqrt(Sigma);
var SigInvSqrt = symInvSqrt(Sigma);
var inner = symCongruence(SigInvSqrt, S);
return symCongruence(SigSqrt, symLog(inner));
}
function expAI(Sigma, U) {
var SigSqrt = symSqrt(Sigma);
var SigInvSqrt = symInvSqrt(Sigma);
var inner = symCongruence(SigInvSqrt, U);
return symCongruence(SigSqrt, symExp(inner));
}
function meanObjectiveAI(Sigma, data) {
var total = 0;
for (var i = 0; i < data.length; i++) {
var d = distAI(Sigma, data[i]);
total += d*d;
}
return total / data.length;
}
function runRGDdemo(n, spread, seed, lr, maxIter) {
var 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);
function normal() {
var u1 = Math.max(rng(), 1e-12);
var u2 = rng();
return Math.sqrt(-2*Math.log(u1)) * Math.cos(2*Math.PI*u2);
}
// The target has distinct eigenvalues and a moderate condition number.
var meanLog = [0.8, 0.35, 0.3];
var SigmaStar = symExp(meanLog);
var starSqrt = symSqrt(SigmaStar);
// Generate whitened tangent perturbations and center their sample average
// exactly at zero. Off-diagonal entries are divided by sqrt(2) so spread
// acts comparably in the Frobenius norm.
var W = [];
var Wbar = [0, 0, 0];
for (var i = 0; i < n; i++) {
var wi = [spread*normal(), spread*normal()/Math.sqrt(2), spread*normal()];
W.push(wi);
Wbar = [Wbar[0]+wi[0]/n, Wbar[1]+wi[1]/n, Wbar[2]+wi[2]/n];
}
var data = [];
for (var i = 0; i < n; i++) {
var centered = [W[i][0]-Wbar[0], W[i][1]-Wbar[1], W[i][2]-Wbar[2]];
data.push(symCongruence(starSqrt, symExp(centered)));
}
var targetAvgLog = [0, 0, 0];
for (var i = 0; i < n; i++) {
var targetLog = logAI(SigmaStar, data[i]);
targetAvgLog = [
targetAvgLog[0] + targetLog[0]/n,
targetAvgLog[1] + targetLog[1]/n,
targetAvgLog[2] + targetLog[2]/n
];
}
var targetResidual = symFrob(symCongruence(symInvSqrt(SigmaStar), targetAvgLog));
var SigmaK = [1, 0, 1];
var sumLog = [0, 0, 0];
for (var i = 0; i < n; i++) {
var Li = symLog(data[i]);
sumLog = [sumLog[0]+Li[0], sumLog[1]+Li[1], sumLog[2]+Li[2]];
}
var SigmaLE = symExp([sumLog[0]/n, sumLog[1]/n, sumLog[2]/n]);
var trajectory = [];
var gradNorms = [];
var residualNorms = [];
var dists = [];
var objectives = [];
var totalBacktracks = 0;
var acceptedStep = 0;
var stoppedByLineSearch = false;
for (var iter = 0; iter <= maxIter; iter++) {
var avgLog = [0, 0, 0];
for (var i = 0; i < n; i++) {
var vi = logAI(SigmaK, data[i]);
avgLog = [avgLog[0]+vi[0]/n, avgLog[1]+vi[1]/n, avgLog[2]+vi[2]/n];
}
// Intrinsic norm: ||avgLog||_Sigma = ||Sigma^{-1/2} avgLog Sigma^{-1/2}||_F.
var whiteAvgLog = symCongruence(symInvSqrt(SigmaK), avgLog);
var residualNorm = symFrob(whiteAvgLog);
var gradNorm = 2*residualNorm;
var distToTarget = distAI(SigmaK, SigmaStar);
var objective = meanObjectiveAI(SigmaK, data);
gradNorms.push(gradNorm);
residualNorms.push(residualNorm);
dists.push(distToTarget);
objectives.push(objective);
trajectory.push({
iter: iter,
Sigma: [SigmaK[0], SigmaK[1], SigmaK[2]],
gradNorm: gradNorm,
residualNorm: residualNorm,
dist: distToTarget,
objective: objective,
acceptedStep: acceptedStep
});
if (residualNorm < 1e-8 || iter === maxIter) break;
var trialEta = lr;
var accepted = false;
for (var bt = 0; bt < 25; bt++) {
// -eta*grad F = 2*eta*avgLog.
var step = [2*trialEta*avgLog[0], 2*trialEta*avgLog[1], 2*trialEta*avgLog[2]];
var candidate = expAI(SigmaK, step);
var candidateObjective = meanObjectiveAI(candidate, data);
if (candidateObjective <= objective - 1e-4*trialEta*gradNorm*gradNorm) {
SigmaK = candidate;
acceptedStep = trialEta;
totalBacktracks += bt;
accepted = true;
break;
}
trialEta *= 0.5;
}
if (!accepted) {
stoppedByLineSearch = true;
break;
}
}
var dataDists = [];
for (var i = 0; i < n; i++) {
dataDists.push(distAI(data[i], SigmaStar));
}
return {
data: data,
dataDists: dataDists,
SigmaStar: SigmaStar,
SigmaLE: SigmaLE,
trajectory: trajectory,
gradNorms: gradNorms,
residualNorms: residualNorms,
dists: dists,
objectives: objectives,
n: n,
spread: spread,
lr: lr,
targetResidual: targetResidual,
iterations: trajectory.length - 1,
totalBacktracks: totalBacktracks,
stoppedByLineSearch: stoppedByLineSearch
};
}
rgdRes = runRGDdemo(spd_n, spd_spread, spd_seed, spd_lr, spd_maxiter);
function renderRGDDemo(res) {
// Helper to represent an SPD matrix as an ellipse
// Σ = R Λ R^T, ellipse axes = sqrt(λ_i) * R columns
function ellipsePath(Sigma, cx, cy, scale) {
var eig = symEig2x2(Sigma);
var lam1 = Math.sqrt(Math.max(0, eig.vals[0]));
var lam2 = Math.sqrt(Math.max(0, eig.vals[1]));
var u1 = eig.vecs[0][0], v1 = eig.vecs[0][1];
// Angle of first eigenvector
var angle = Math.atan2(v1, u1);
// Build SVG ellipse
var rx = lam1 * scale;
var ry = lam2 * scale;
return '<ellipse cx="'+cx.toFixed(1)+'" cy="'+cy.toFixed(1)+'" rx="'+rx.toFixed(2)+'" ry="'+ry.toFixed(2)+'" transform="rotate('+(-angle*180/Math.PI).toFixed(2)+','+cx.toFixed(1)+','+cy.toFixed(1)+')" />';
}
var svgW = 440, svgH = 400;
var cx = 220, cy = 210;
var matricesForScale = res.data
.concat([res.SigmaStar, res.SigmaLE])
.concat(res.trajectory.map(function(x) { return x.Sigma; }));
var maxAxis = 1;
for (var s = 0; s < matricesForScale.length; s++) {
maxAxis = Math.max(maxAxis, Math.sqrt(Math.max(0, symEig2x2(matricesForScale[s]).vals[0])));
}
var ellScale = Math.min(105, 172/maxAxis);
var parts = [];
// Background grid
parts.push('<rect x="0" y="0" width="'+svgW+'" height="'+svgH+'" fill="#fafbfc" rx="4"/>');
// Draw data ellipses (faded)
var maxDist = Math.max.apply(null, res.dataDists.concat([0.01]));
for (var i = 0; i < res.data.length; i++) {
var d = res.dataDists[i];
var t = Math.min(1, d / maxDist);
var r = Math.round(40 + 200 * t);
var g = Math.round(100 + 80 * (1-t));
var b = Math.round(200 - 60 * t);
var opacity = 0.35 + 0.35 * (1-t);
var ep = ellipsePath(res.data[i], cx, cy, ellScale);
parts.push('<g opacity="'+opacity.toFixed(2)+'">');
parts.push(ep.replace('/>', ' fill="rgb('+r+','+g+','+b+')" stroke="rgba(100,100,100,0.35)" stroke-width="0.7"/>'));
parts.push('</g>');
}
// Log-Euclidean solution
var epLE = ellipsePath(res.SigmaLE, cx, cy, ellScale);
parts.push(epLE.replace('/>', ' fill="none" stroke="#1971c2" stroke-width="2.8" stroke-dasharray="8,5" opacity="0.85"/>'));
parts.push('<text x="'+(cx+85)+'" y="26" font-size="11" fill="#1971c2" font-weight="bold">– – Σ̂_LE (closed form)</text>');
// Earlier RGD iterates.
if (res.trajectory.length > 1) {
var stride = Math.max(1, Math.floor((res.trajectory.length-1)/20));
for (var t = 0; t < res.trajectory.length - 1; t += stride) {
var epT = ellipsePath(res.trajectory[t].Sigma, cx, cy, ellScale);
parts.push(epT.replace('/>', ' fill="none" stroke="#e03131" stroke-width="0.8" opacity="0.3"/>'));
}
}
// Draw the exact mean as a wide green halo, then the final iterate on top.
// Both remain visible even when convergence makes the ellipses coincide.
var epStar = ellipsePath(res.SigmaStar, cx, cy, ellScale);
parts.push(epStar.replace('/>', ' fill="none" stroke="#2b8a3e" stroke-width="5.2" opacity="0.9"/>'));
var last = res.trajectory[res.trajectory.length - 1];
var epLast = ellipsePath(last.Sigma, cx, cy, ellScale);
parts.push(epLast.replace('/>', ' fill="none" stroke="#e03131" stroke-width="2.2" opacity="0.95"/>'));
// Legend
parts.push('<circle cx="15" cy="'+ (svgH - 55) +'" r="5" fill="#2b8a3e" stroke="#fff" stroke-width="1.5"/><text x="26" y="'+ (svgH - 51) +'" font-size="11" fill="#2b8a3e" font-weight="bold">Σ* (exact AI mean)</text>');
parts.push('<circle cx="15" cy="'+ (svgH - 35) +'" r="5" fill="#e03131" stroke="#fff" stroke-width="1.5"/><text x="26" y="'+ (svgH - 31) +'" font-size="11" fill="#e03131" font-weight="bold">Σ⁽ᵏ⁾ (RGD iterate)</text>');
parts.push('<line x1="7" y1="'+(svgH-18)+'" x2="23" y2="'+(svgH-18)+'" stroke="#1971c2" stroke-width="2" stroke-dasharray="6,4"/><text x="26" y="'+(svgH-14)+'" font-size="11" fill="#1971c2" font-weight="bold">Σ̂_LE</text>');
return '<svg width="'+svgW+'" height="'+svgH+'" viewBox="0 0 '+svgW+' '+svgH+'" style="display:block; width:100%; max-width:440px; height:auto; border:1px solid #dee2e6; border-radius:4px; background:#fafbfc;">\n'+parts.join('\n')+'\n</svg>';
}
viewof rgd_view = {
const res = rgdRes;
const last = res.trajectory[res.trajectory.length - 1];
const root = html`
<div style="font-family: system-ui, sans-serif; max-width: 920px;">
<div style="display: flex; gap: 20px; flex-wrap: wrap;">
<div style="flex:1 1 440px; min-width:0;">
<h4>SPD Cone: Data and RGD Iterates</h4>
${html([renderRGDDemo(res)])}
<p style="font-size:0.8em; color:#6c757d; margin-top:4px; max-width:440px;">
Each ellipse represents a 2×2 SPD matrix via its covariance ellipse Σ = QΛQᵀ (axes = √λᵢ, directions = eigenvectors).
Faded ellipses = data; <b style="color:#2b8a3e;">green</b> = exact sample AI mean;
<b style="color:#e03131;">red</b> = RGD iterate.
</p>
</div>
<div style="flex: 1; min-width: 300px;">
<div style="padding: 12px; background: #f8f9fa; border-radius: 6px; margin-bottom: 12px;">
<h4 style="margin-top:0;">Convergence Diagnostics</h4>
<table style="width:100%; border-collapse:collapse; font-size:0.88em;">
<tr><td style="padding:3px 8px;">Sample size n</td>
<td style="padding:3px 8px; text-align:right;">${res.n}</td></tr>
<tr><td style="padding:3px 8px;">Spread σ</td>
<td style="padding:3px 8px; text-align:right;">${res.spread.toFixed(2)}</td></tr>
<tr><td style="padding:3px 8px;">Initial step size η</td>
<td style="padding:3px 8px; text-align:right;">${res.lr.toFixed(2)}</td></tr>
<tr><td style="padding:3px 8px;">Last accepted η</td>
<td style="padding:3px 8px; text-align:right;">${last.acceptedStep ? (last.acceptedStep < 1e-3 ? last.acceptedStep.toExponential(2) : last.acceptedStep.toFixed(4)) : '—'}</td></tr>
<tr><td colspan="2"><hr style="margin:4px 0;"></td></tr>
<tr><td style="padding:3px 8px;">Final ‖∇F<sub>n</sub>‖</td>
<td style="padding:3px 8px; text-align:right; font-family:monospace;">${last.gradNorm.toExponential(2)}</td></tr>
<tr><td style="padding:3px 8px;">Accepted updates</td>
<td style="padding:3px 8px; text-align:right;">${res.iterations}</td></tr>
<tr><td style="padding:3px 8px;">Backtracking reductions</td>
<td style="padding:3px 8px; text-align:right;">${res.totalBacktracks}</td></tr>
<tr><td style="padding:3px 8px;">Constructed residual at Σ*</td>
<td style="padding:3px 8px; text-align:right; font-family:monospace;">${res.targetResidual.toExponential(2)}</td></tr>
<tr><td style="padding:3px 8px;">d<sub>AI</sub>(Σ<sup>(k)</sup>, Σ*)</td>
<td style="padding:3px 8px; text-align:right; font-family:monospace; font-weight:bold;">${last.dist.toExponential(2)}</td></tr>
<tr><td style="padding:3px 8px;">d<sub>AI</sub>(Σ̂<sub>LE</sub>, Σ*)</td>
<td style="padding:3px 8px; text-align:right; font-family:monospace;">${distAI(res.SigmaLE, res.SigmaStar).toExponential(2)}</td></tr>
</table>
</div>
<div style="padding: 12px; border-radius: 6px; ${last.residualNorm < 1e-8 ? 'background:#d3f9d8; border-left:4px solid #2b8a3e;' : 'background:#fff3cd; border-left:4px solid #f08c00;'}">
<b>First-order residual:</b> ‖(1/n) Σ log<sub>Σ̂</sub>(Σ<sub>i</sub>)‖<sub>Σ̂</sub> = ${last.residualNorm.toExponential(2)}
${last.residualNorm < 1e-8
? '<span style="color:#2b8a3e;"> ≈ 0 ✓ Converged</span>'
: (res.stoppedByLineSearch
? '<span style="color:#f08c00;"> → line search could not find an accepted step</span>'
: '<span style="color:#f08c00;"> → increase the iteration limit</span>')}
</div>
</div>
</div>
<div style="margin-top:16px;">
<h4>Gradient Norm and Distance to the Exact AI Mean</h4>
<div style="display:flex; flex-wrap:wrap; gap:20px;">
<svg width="430" height="180" viewBox="0 0 430 180" style="flex:1 1 380px; width:100%; max-width:430px; height:auto; border:1px solid #dee2e6; border-radius:4px;">
${(() => {
var w=430, h=180, mx=45, my=150, gw=370, gh=125;
var gn = res.gradNorms;
var parts = [];
// axes
parts.push('<line x1="'+mx+'" y1="'+my+'" x2="'+(mx+gw)+'" y2="'+my+'" stroke="#495057" stroke-width="1"/>');
parts.push('<line x1="'+mx+'" y1="'+my+'" x2="'+mx+'" y2="'+(my-gh)+'" stroke="#495057" stroke-width="1"/>');
parts.push('<text x="'+(mx+gw/2)+'" y="'+(h-5)+'" text-anchor="middle" font-size="10">Iteration</text>');
parts.push('<text x="12" y="'+(my-gh/2)+'" text-anchor="middle" font-size="10" transform="rotate(-90,12,'+(my-gh/2)+')">‖∇F‖</text>');
// log scale y-axis (approximate)
var maxGN = Math.max.apply(null, gn.concat([1e-16]));
var positiveGN = gn.filter(function(x) { return x > 0; });
var minGN = positiveGN.length ? Math.max(1e-16, Math.min.apply(null, positiveGN)) : 1e-16;
var logMax = Math.log10(maxGN), logMin = Math.log10(minGN);
var logRange = logMax - logMin;
if (logRange < 0.5) { logMax = logMin + 0.5; logRange = 0.5; }
for (var lev = Math.ceil(logMin); lev <= Math.floor(logMax); lev++) {
var y = my - gh * (lev - logMin) / logRange;
parts.push('<line x1="'+(mx-3)+'" y1="'+y.toFixed(1)+'" x2="'+mx+'" y2="'+y.toFixed(1)+'" stroke="#adb5bd" stroke-width="0.7"/>');
parts.push('<text x="'+(mx-5)+'" y="'+(y+4).toFixed(0)+'" text-anchor="end" font-size="9" fill="#868e96">10^'+lev+'</text>');
}
// plot
var path = '';
for (var i = 0; i < gn.length; i++) {
var x = mx + gw * i / Math.max(gn.length-1, 1);
var ly = Math.max(logMin, Math.min(logMax, Math.log10(Math.max(gn[i], 1e-16))));
var y = my - gh * (ly - logMin) / logRange;
path += (i===0 ? 'M' : 'L') + x.toFixed(1) + ',' + y.toFixed(1) + ' ';
}
parts.push('<path d="'+path+'" fill="none" stroke="#e03131" stroke-width="2"/>');
return parts.join('\n');
})()}
</svg>
<svg width="430" height="180" viewBox="0 0 430 180" style="flex:1 1 380px; width:100%; max-width:430px; height:auto; border:1px solid #dee2e6; border-radius:4px;">
${(() => {
var w=430, h=180, mx=45, my=150, gw=370, gh=125;
var ds = res.dists;
var parts = [];
parts.push('<line x1="'+mx+'" y1="'+my+'" x2="'+(mx+gw)+'" y2="'+my+'" stroke="#495057" stroke-width="1"/>');
parts.push('<line x1="'+mx+'" y1="'+my+'" x2="'+mx+'" y2="'+(my-gh)+'" stroke="#495057" stroke-width="1"/>');
parts.push('<text x="'+(mx+gw/2)+'" y="'+(h-5)+'" text-anchor="middle" font-size="10">Iteration</text>');
parts.push('<text x="12" y="'+(my-gh/2)+'" text-anchor="middle" font-size="10" transform="rotate(-90,12,'+(my-gh/2)+')">d_AI(Σ⁽ᵏ⁾,Σ*)</text>');
var maxD = Math.max.apply(null, ds.concat([1e-12]));
var dRange = maxD;
var yTicks = 4;
for (var t = 0; t <= yTicks; t++) {
var dv = maxD * (yTicks-t) / yTicks;
var y = my - gh * dv / dRange;
parts.push('<line x1="'+(mx-3)+'" y1="'+y.toFixed(1)+'" x2="'+mx+'" y2="'+y.toFixed(1)+'" stroke="#adb5bd" stroke-width="0.7"/>');
parts.push('<text x="'+(mx-5)+'" y="'+(y+4).toFixed(0)+'" text-anchor="end" font-size="9" fill="#868e96">'+dv.toFixed(2)+'</text>');
}
var path = '';
for (var i = 0; i < ds.length; i++) {
var x = mx + gw * i / Math.max(ds.length-1, 1);
var y = my - gh * ds[i] / dRange;
path += (i===0 ? 'M' : 'L') + x.toFixed(1) + ',' + y.toFixed(1) + ' ';
}
parts.push('<path d="'+path+'" fill="none" stroke="#2b8a3e" stroke-width="2"/>');
return parts.join('\n');
})()}
</svg>
</div>
</div>
</div>`;
return root;
}
```
::: {.callout-note title="How to read this"}
- The <b style="color:#e03131;">red ellipse</b> is the final RGD iterate $\Sigma^{(k)}$. Each accepted update follows one geodesic segment; the full multi-step trajectory need not be a single geodesic.
- The <b style="color:#1971c2;">dashed blue ellipse</b> $\hat{\Sigma}_{\mathrm{LE}}$ is the log-Euclidean solution (one step). Notice it may differ from the affine-invariant optimum — the two metrics define different geometries.
- The **gradient norm** plot reports observed behavior; it is not itself proof of a universal linear rate.
- The slider sets the initial $\eta$. Armijo backtracking halves it until the objective-decrease test succeeds.
:::
::: {.callout-tip title="Try these experiments"}
- **Small initial η (0.1–0.3):** Accepted steps are usually close to the requested value, but convergence may be slow.
- **Moderate initial η (0.5–1.0):** Often efficient for these generated matrices. This is an empirical observation, not a universal optimum.
- **Large initial η (> 1.2):** Backtracking may reject the proposed step and reduce η. Watch the accepted-step and backtracking diagnostics.
- **Increase spread:** Larger spread means data points are farther apart. The Fréchet mean remains unique (Hadamard property), but convergence may require more iterations.
- **Compare LE vs AI:** The log-Euclidean solution $\hat{\Sigma}_{\mathrm{LE}}$ (dashed blue) may differ visibly from the affine-invariant solution. This is because the two metrics define different geodesics and different notions of "straight-line average." In practice, the affine-invariant mean is preferred when invariance under congruence $\Sigma \mapsto G\Sigma G^\top$ matters.
:::
### Connection to Earlier Applications
Each of the applications from Lectures 2–9 can now be understood through the manifold lens:
- **Fréchet ANOVA (Lecture 4):** The group means are Fréchet means on the SPD manifold. RGD is available when the chosen Riemannian metric provides computable derivatives and maps.
- **Global Fréchet regression (Lecture 6):** The regression estimator $\hat{\mu}(x)$ minimizes a weighted squared-distance criterion. Where differentiable, its first-order equation is a weighted sum of log maps; negative regression weights can destroy convexity and require extra care.
- **Kernel and local-linear regression (Lectures 7–8):** The same weighted optimization primitive appears. RGD can be used under the affine-invariant metric when the weighted objective is well posed, although a global convergence guarantee does not follow merely from the unweighted Hadamard result.
- **TV-regularized estimation (Lecture 9):** The proximal operator for TV regularization on a manifold requires solving a geodesically-constrained optimization. The manifold framework provides the geodesic interpolation and proximal map.
The conceptual shift is that many earlier metric-space calculations admit a differential interpretation on smooth manifolds. Log-Euclidean flatness makes the SPD mean Euclidean in log coordinates, whereas the affine-invariant metric leads to genuinely Riemannian computation.
## Key Takeaways
- The **first-order condition** $\mathbb{E}\{\log_\mu(X)\}=0$ is necessary at a differentiable Fréchet minimizer; it is sufficient only with suitable convexity or uniqueness conditions.
- **Existence** holds on any complete connected Riemannian manifold (BP 2003).
- **Uniqueness** is global on Hadamard manifolds, while Afsari gives a small-support guarantee under an upper curvature bound.
- The **parametric rate** $O_P(n^{-1/2})$ follows in finite-dimensional Hadamard manifolds from quadratic growth, finite second moments, consistency, and local entropy.
- The **cut locus** is the main source of nonsmoothness for squared-distance terms and must be handled explicitly in first-order arguments.
## Exercises
1. **First-order condition check.** For the sample mean on $S^1$, verify numerically that $\sum \log_{\hat{\mu}}(x_i) = 0$ for data within a semicircle. <a href="javascript:void(0)" onclick="showSolution('l16-sol-1')" class="solution-link">📝 Show Solution</a>
2. **Afsari radius on $S^2$.** For $S^2$ ($\kappa = 1$, $\operatorname{inj} = \pi$), compute the maximum ball radius $\rho_\kappa$ guaranteeing unique Fréchet means. <a href="javascript:void(0)" onclick="showSolution('l16-sol-2')" class="solution-link">📝 Show Solution</a>
3. **Hadamard vs. sphere.** Give an explicit example of a distribution on $S^2$ with multiple Fréchet means. Why does the Hadamard case avoid this? <a href="javascript:void(0)" onclick="showSolution('l16-sol-3')" class="solution-link">📝 Show Solution</a>
4. **Rate on Hadamard manifolds.** Explain how the $\operatorname{CAT}(0)$ variance inequality and finite-dimensional local entropy together yield the $n^{-1/2}$ rate. <a href="javascript:void(0)" onclick="showSolution('l16-sol-4')" 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="l16-sol-1"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 1</h4><button onclick="closeSolution('l16-sol-1')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Verify the first-order condition numerically.</div><div class="solution-answer"><strong>Solution:</strong> Suppose the observations lie in an open semicircle. Choose an angular coordinate there and unwrap the observations to real numbers $\theta_i$ in an interval of length $<\pi$. Their intrinsic sample mean is represented by $\hat\theta=n^{-1}\sum_i\theta_i$. Because no observation is antipodal to $\hat\mu$, $\log_{\hat\mu}(x_i)$ is the signed angular displacement $\theta_i-\hat\theta$, and therefore $\sum_i\log_{\hat\mu}(x_i)=\sum_i(\theta_i-\hat\theta)=0$. Outside this convex arc, a minimizer still satisfies the equation whenever no datum lies in its cut locus, but unwrapping and uniqueness require separate checks.</div></div></dialog>
<dialog id="l16-sol-2"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 2</h4><button onclick="closeSolution('l16-sol-2')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Compute the Afsari radius for $S^2$.</div><div class="solution-answer"><strong>Solution:</strong> For the unit sphere, $\kappa=1$ and $\operatorname{inj}(S^2)=\pi$. Hence $\rho_\kappa=\frac12\min\{\pi,\pi/\sqrt1\}=\pi/2$. Afsari's theorem applies when the support is contained in some ball of radius strictly below $\pi/2$. For a finite sample, containment in an open hemisphere gives a positive margin from its boundary and hence such a smaller ball. Failure of this sufficient condition does not imply nonuniqueness.</div></div></dialog>
<dialog id="l16-sol-3"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 3</h4><button onclick="closeSolution('l16-sol-3')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Example of multiple Fréchet means on $S^2$.</div><div class="solution-answer"><strong>Solution:</strong> Put mass $1/2$ at the north pole and mass $1/2$ at the south pole. If $x$ has colatitude $\theta$, then $F(x)=\frac12\theta^2+\frac12(\pi-\theta)^2=(\theta-\pi/2)^2+\pi^2/4$. Thus the global minimum occurs exactly when $\theta=\pi/2$: every equatorial point is a Fréchet mean. On a Hadamard manifold, the $\operatorname{CAT}(0)$ inequality makes the integrated squared-distance objective strongly geodesically convex, ruling out two distinct minimizers.</div></div></dialog>
<dialog id="l16-sol-4"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 4</h4><button onclick="closeSolution('l16-sol-4')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Why $n^{-1/2}$ on Hadamard manifolds?</div><div class="solution-answer"><strong>Solution:</strong> The $\operatorname{CAT}(0)$ inequality implies the variance inequality $F(q)-F(\mu)\ge d^2(q,\mu)$, so the population objective has quadratic growth (growth exponent $2$). In a finite-dimensional manifold, a sufficiently small ball has covering number $N(\varepsilon,B_\delta,d)\le C(\delta/\varepsilon)^m$. Consequently $\sqrt{\log N}$ is bounded by $C_\alpha(\delta/\varepsilon)^\alpha$ for any chosen $\alpha>0$, and one may take $\alpha<1$ in the general rate theorem. With a finite second moment and consistency to localize the estimator, these ingredients give $d(\hat\mu_n,\mu)=O_P(n^{-1/2})$. The growth exponent and entropy exponent are different quantities and should not both be denoted by the same symbol.</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
- @BhattacharyaPatrangenaru2003 — Foundational theory of Fréchet means on manifolds.
- @Afsari2011 — Uniqueness and convexity of Fréchet means on Riemannian manifolds.
- @Sturm2003 — Barycenters and variance inequalities in nonpositively curved metric spaces.
- @Schotz2019 — Empirical Fréchet-mean rates from growth, moment, and entropy conditions.
- @PennecFillardAyache2006 — Fréchet means on the SPD manifold.
- @BhattacharyaPatrangenaru2005 — CLT for intrinsic Fréchet means.
## Self-Assessment Quiz
Test your understanding of this lecture with the interactive MCQ quiz:
👉 **[Lecture 16 Quiz — 10 Multiple-Choice Questions](../quizzes/lecture-16-quiz.qmd)**