---
title: "Lecture 18: Riemannian Manifolds — Geodesic Regression"
subtitle: "Fitting geodesic curves to manifold-valued responses"
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:
- Formulate the geodesic regression model as the manifold analogue of simple linear regression.
- Write the geodesic least-squares objective and describe how it is minimized.
- Define symmetric spaces and give key examples (Euclidean, sphere, SPD, $\mathrm{SO}(m)$).
- Compute and interpret the intrinsic $R^2$ and its permutation test.
## Geodesic Regression Model
When the response space is a Riemannian manifold, the most direct parametric analogue of linear regression is to assume the regression mean follows a geodesic. This was developed by @Fletcher2013.
For a scalar predictor $x \in \mathbb{R}$ and manifold-valued response $Y \in \mathcal{M}$, fix $p \in \mathcal{M}$ (intercept) and $v \in T_p\mathcal{M}$ (slope). The geodesic is
$$
\gamma_{p,v}(x) = \exp_p(x v).
$$
The stochastic model is
$$
Y = \exp_{\gamma_{p,v}(X)}(\varepsilon), \qquad \varepsilon \in T_{\gamma_{p,v}(X)}\mathcal{M},
$$
where the error is represented in the tangent space at the fitted point. When the logarithm is uniquely defined, a natural centering condition is
$$
\mathbb{E}\!\left\{\log_{\gamma_{p,v}(X)}(Y)\mid X\right\}=0.
$$
This is a tangent-space centering condition; by itself it does not guarantee that $\gamma_{p,v}(x)$ is the unique conditional Fréchet mean unless the conditional Fréchet function is suitably convex. In Euclidean space, $\exp_p(u) = p + u$, so the model reduces to $Y = p + Xv + \varepsilon$.
## Geodesic Least-Squares
Given observations $(x_i, y_i) \in \mathbb{R} \times \mathcal{M}$, the objective is
$$
E_n(p, v) = \frac{1}{2}\sum_{i=1}^n d^2(y_i, \exp_p(x_i v)), \qquad (p, v) \in T\mathcal{M}.
$$
Any minimizer $(\hat{p}, \hat{v})$ gives the fitted curve $\hat{\gamma}(x) = \exp_{\hat{p}}(x\hat{v})$.
**Computation.** Writing $\mu_i = \exp_p(x_i v)$ and $\epsilon_i = \log_{\mu_i}(y_i) \in T_{\mu_i}\mathcal{M}$, Fletcher derived:
$$
\nabla_p E_n = -\sum_{i=1}^n D_1(\exp)_{(p, x_i v)}^* \epsilon_i, \qquad
\nabla_v E_n = -\sum_{i=1}^n x_i\, D_2(\exp)_{(p, x_i v)}^* \epsilon_i,
$$
where $D_1, D_2$ are differentials of the exponential map with respect to base point and tangent argument, and $*$ denotes the adjoint. Estimation uses **intrinsic gradient descent** on the tangent bundle.
For $D_1$, the tangent vector $v$ is extended by parallel transport while the base point varies; this specifies the horizontal variation needed because $(p,v)$ belongs to the tangent bundle rather than a Cartesian product.
**Intrinsic $R^2$.**
Fletcher's intrinsic coefficient of determination:
$$
R^2 = 1 - \frac{n^{-1}\sum_i d^2(y_i, \hat{\gamma}(x_i))}{\inf_\omega n^{-1}\sum_i d^2(y_i, \omega)},
$$
where the denominator is the sample Fréchet variance of the responses. This reduces to the usual $R^2$ in Euclidean space. A **permutation test** for $H_0: v = 0$ (no trend) permutes the $x_i$ relative to $y_i$, recomputes $R^2$, and compares the observed value to the permutation distribution.
For a global least-squares solution, the constant curve through a sample Fréchet mean is always an admissible competitor, so $0\le R^2\le 1$. A merely local or approximate numerical solution need not preserve the lower bound.
## Existence and Uniqueness on Symmetric Spaces
::: {#def-symmetric-space .definition title="Symmetric space"}
A connected Riemannian manifold $\mathcal{M}$ is a **symmetric space** if for every $x \in \mathcal{M}$, there exists an involutive isometry $\sigma_x$ (a "geodesic symmetry") such that $\sigma_x(\gamma(t)) = \gamma(-t)$ for every geodesic $\gamma$ with $\gamma(0) = x$.
:::
::: {.callout-note collapse="true" title="Involutive Isometry"}
An isometry $f : \mathcal{M} \to \mathcal{M}$ is **involutive** if applying it twice yields the identity: $f \circ f = \mathrm{id}_\mathcal{M}$. In other words, $f$ is its own inverse — $f^{-1} = f$. The geodesic symmetry $\sigma_x$ at a point $x$ satisfies $\sigma_x \circ \sigma_x = \mathrm{id}$, reflecting the manifold locally.
:::
Key properties: symmetric spaces are **homogeneous** (look the same everywhere) and **complete** [@Lee2018].
**Examples.**
| Space | Geodesic symmetry at $p$ |
|---|---|
| $\mathbb{R}^m$ | $\sigma_x(y) = 2x - y$ (reflection) |
| $\mathbb{S}^m$ | $\sigma_x(y) = 2\langle x,y\rangle x - y$ |
| $\mathrm{SO}(m)$ | $\sigma_I(g) = g^{-1}$ (inversion) |
| $\mathcal{S}_{++}^m$ (affine-invariant) | $\sigma_P(\Sigma) = P\Sigma^{-1}P$ |
| Grassmann, projective spaces | Natural involutive isometries |
Fix data $\mathcal D=\{(x_i,y_i)\}_{i=1}^n\subset\mathbb R\times\mathcal M$, with $n\ge2$ and the $x_i$ not all equal. Translating every predictor by the same constant only reparametrizes the same fitted geodesics, so choose an index $i_0$ and take $x_{i_0}=0$ without loss of generality.
- For $r>0$, define the admissible class of **$r$-close geodesic parameters**
$$
\mathcal G_r(\mathcal D)
=
\left\{(p,v)\in T\mathcal M:
\begin{array}{l}
d\!\left(\exp_p(x_i v),y_i\right)\le r,\\[-2pt]
y_i\notin\operatorname{Cut}\!\left(\exp_p(x_i v)\right)
\end{array}
\text{ for every }i
\right\}.
$$
Thus “unique among $r$-close geodesics” means uniqueness of $(p,v)$ when $E_n$ is restricted to $\mathcal G_r(\mathcal D)$; it does not rule out a lower value outside this class.
- Let $\mathrm{vol}$ be Riemannian volume and define
$$
Z(q,\sigma)=\int_{\mathcal M}\exp\!\left\{-\frac{d^2(q,y)}{2\sigma^2}\right\}\,d\mathrm{vol}(y).
$$
Homogeneity of a symmetric space makes this integral independent of $q$; write its common value as $Z(\sigma)$. Whenever $Z(\sigma)<\infty$, the **intrinsic isotropic Gaussian** with center $q$ and scale $\sigma>0$ is
$$
f(y;q,\sigma)=Z(\sigma)^{-1}
\exp\!\left\{-\frac{d^2(q,y)}{2\sigma^2}\right\}.
$$
::: {.callout-note collapse="true" title="Riemannian volume"}
Let $(\mathcal M,g)$ be a $d$-dimensional Riemannian manifold. In local coordinates $z=(z^1,\ldots,z^d)$, write
$$
g_{ab}(z)=g_z\!\left(\frac{\partial}{\partial z^a},
\frac{\partial}{\partial z^b}\right)
$$
for the metric matrix. The **Riemannian volume measure** is the measure whose coordinate expression is
$$
d\mathrm{vol}_g(z)
=
\sqrt{\det\{g_{ab}(z)\}}\,
dz^1\cdots dz^d.
$$
Thus, for a nonnegative or integrable function $h$ supported in a coordinate neighborhood $U$,
$$
\int_U h(y)\,d\mathrm{vol}_g(y)
=
\int_{z(U)}
h\!\left(z^{-1}(u)\right)
\sqrt{\det\{g_{ab}(u)\}}\,du.
$$
The factor $\sqrt{\det(g_{ab})}$ transforms with the coordinate Jacobian, so these local expressions define a chart-independent measure.
:::
::: {#thm-fletcher .theorem title="Geodesic regression on symmetric spaces (Fletcher 2013, Theorems 2–4)"}
Let $(\mathcal M,g)$ be a connected Riemannian symmetric space, and let $\mathcal D=\{(x_i,y_i)\}_{i=1}^n$ satisfy the nondegenerate-design and predictor-origin conditions above. Suppose that the responses $y_1,\ldots,y_n$ lie in a geodesically convex set $A\subset\mathcal M$. Assume one of the following two curvature-and-concentration conditions:
1. **Positive-curvature case:** the sectional curvature is positive and upper bounded by $K_M>0$, and
$$
\operatorname{diam}(A)<r,\qquad r=\frac{\pi}{4\sqrt{K_M}}.
$$
2. **Negative-curvature case:** the sectional curvature is negative and lower bounded by $K_M<0$, and
$$
\operatorname{diam}(A)<r,\qquad r=\frac{1}{\sqrt{|K_M|}}.
$$
Then $\mathcal G_r(\mathcal D)$ is nonempty, and the restricted least-squares problem
$$
(\hat p,\hat v)
=
\underset{(p,v)\in\mathcal G_r(\mathcal D)}{\arg\min}\;
E_n(p,v),
\qquad
E_n(p,v)=\frac12\sum_{i=1}^n
d^2\!\left(y_i,\exp_p(x_i v)\right),
$$
has a unique solution.
If, in addition, $Z(\sigma)<\infty$ and the observations are conditionally independent with
$$
Y_i\mid X_i=x_i\sim
f\!\left(\,\cdot\,;\exp_p(x_i v),\sigma\right),
$$
then, for each fixed $\sigma>0$, the conditional log-likelihood on $\mathcal G_r(\mathcal D)$ is
$$
\ell(p,v;\sigma)
=-n\log Z(\sigma)-\frac{1}{\sigma^2}E_n(p,v).
$$
Consequently, its unique maximizer in $(p,v)$ is the same pair $(\hat p,\hat v)$ as the restricted least-squares minimizer.
:::
The result is deliberately local and qualified. It does not assert global uniqueness outside $\mathcal G_r(\mathcal D)$, and its two curvature cases do not cover symmetric spaces with zero or mixed sectional curvature. In the flat Euclidean case, ordinary least squares gives uniqueness when the predictors are not all equal.
## Interactive Exploration: Geodesic Regression on $S^2$
```{ojs}
//| label: fig-geodesic-reg
//| fig-cap: "Interactive: least-squares geodesic regression on the sphere — choose the true regression curve, drag to rotate the 3D globe"
n18_control = Inputs.range([10, 100], {step: 10, value: 40, label: "Sample size n"})
noise18_control = Inputs.range([0.05, 0.5], {step: 0.05, value: 0.15, label: "Noise σ (radians)"})
n18 = Generators.input(n18_control)
noise18 = Generators.input(noise18_control)
// Candidate true regression curves mu(x) on S², x ∈ [-1, 1].
// Geodesic curves carry (p, v); arbitrary curves carry a point-valued mu(x).
curve18_defs = [
{label: "Geodesic: gentle great circle", geodesic: true, p: [0,0,1], v: [0.6,0.3,0]},
{label: "Geodesic: steep great circle", geodesic: true, p: [0,0,1], v: [-0.25,0.85,0]},
{label: "Non-geodesic: small circle (parallel)", geodesic: false,
mu: x => { const la = 1.05, lo = 2.3*x, c = Math.cos(la); return [c*Math.cos(lo), Math.sin(la), c*Math.sin(lo)]; }},
{label: "Non-geodesic: S-bend around the equator", geodesic: false,
mu: x => { const la = 0.9*Math.sin(3*x), lo = 0.8*x, c = Math.cos(la); return [c*Math.cos(lo), Math.sin(la), c*Math.sin(lo)]; }},
{label: "Non-geodesic: Lissajous wiggle", geodesic: false,
mu: x => { const la = 0.8*Math.sin(2*x), lo = 1.2*Math.sin(3*x), c = Math.cos(la); return [c*Math.cos(lo), Math.sin(la), c*Math.sin(lo)]; }}
]
curve18_control = Inputs.select(curve18_defs.map(d => d.label), {value: curve18_defs[0].label, label: "True regression curve"})
curve18 = Generators.input(curve18_control)
curve18_def = curve18_defs.find(d => d.label === curve18) || curve18_defs[0]
geod_controls_view = html`
<style>
#fig-geodesic-reg .quarto-subfloat-caption,
#fig-geod-vs-frechet .quarto-subfloat-caption { display:none; }
.geod-slider-grid { display:flex; flex-wrap:wrap; gap:8px 20px; width:100%; max-width:840px; margin:0 auto 12px; font:0.85em system-ui,sans-serif; container-type:inline-size; }
.geod-slider-grid > * { min-width:0; margin:0; }
.geod-curve-control { flex:1 1 100%; }
.geod-numeric-control { flex:1 1 calc((100% - 20px)/2); }
.geod-slider-grid input[type="number"] { width:10rem !important; }
.geod-slider-grid select { width:100%; max-width:100%; }
@container (max-width:480px) { .geod-numeric-control { flex-basis:100%; } }
</style>
<div class="geod-slider-grid">
<div class="geod-curve-control">${curve18_control}</div>
<div class="geod-numeric-control">${n18_control}</div>
<div class="geod-numeric-control">${noise18_control}</div>
</div>`
function nelderMead18(f, start, steps, maxIter) {
const d = start.length;
let simplex = [start.slice()];
for (let j = 0; j < d; j++) {
const q = start.slice();
q[j] += steps[j];
simplex.push(q);
}
let vals = simplex.map(f);
for (let iter = 0; iter < maxIter; iter++) {
const ord = Array.from({length:d+1}, (_,i) => i).sort((i,j) => vals[i]-vals[j]);
simplex = ord.map(i => simplex[i]);
vals = ord.map(i => vals[i]);
const mean = vals.reduce((a,b) => a+b, 0)/(d+1);
const spread = Math.sqrt(vals.reduce((s,z) => s+(z-mean)**2, 0)/(d+1));
if (spread < 1e-10 * (1 + Math.abs(mean))) break;
const centroid = Array(d).fill(0);
for (let i = 0; i < d; i++) for (let j = 0; j < d; j++) centroid[j] += simplex[i][j]/d;
const trial = factor => centroid.map((z,j) => z + factor*(centroid[j]-simplex[d][j]));
const reflected = trial(1), fr = f(reflected);
if (fr < vals[0]) {
const expanded = trial(2), fe = f(expanded);
if (fe < fr) { simplex[d] = expanded; vals[d] = fe; }
else { simplex[d] = reflected; vals[d] = fr; }
} else if (fr < vals[d-1]) {
simplex[d] = reflected; vals[d] = fr;
} else {
const contracted = fr < vals[d] ? trial(0.5) : centroid.map((z,j) => z + 0.5*(simplex[d][j]-z));
const fc = f(contracted);
if (fc < Math.min(fr, vals[d])) {
simplex[d] = contracted; vals[d] = fc;
} else {
for (let i = 1; i <= d; i++) {
simplex[i] = simplex[0].map((z,j) => z + 0.5*(simplex[i][j]-z));
vals[i] = f(simplex[i]);
}
}
}
}
const best = vals.indexOf(Math.min(...vals));
return {q:simplex[best], value:vals[best]};
}
function runGeodReg(n, sigma, def, 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 || 42);
function rnorm() {
return Math.sqrt(-2*Math.log(Math.max(rng(),1e-15)))*Math.cos(2*Math.PI*rng());
}
function dot(a,b) { return a[0]*b[0]+a[1]*b[1]+a[2]*b[2]; }
function norm(a) { return Math.sqrt(dot(a,a)); }
function cross(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]]; }
function unit(a) { const z = norm(a); return a.map(x => x/z); }
// True regression curve: a geodesic exp_p(xv) or an arbitrary smooth curve mu(x)
const trueMu = def.geodesic ? (x => expS2(def.p, def.v.map(z => z*x))) : def.mu;
function expS2(p, v) {
const r = norm(v);
if (r < 1e-12) return p.slice();
const a = Math.cos(r), b = Math.sin(r)/r;
return [a*p[0]+b*v[0], a*p[1]+b*v[1], a*p[2]+b*v[2]];
}
function sphereDist(p, q) {
return Math.acos(Math.max(-1, Math.min(1, dot(p,q))));
}
// Isotropic Gaussian coordinates in the tangent plane, mapped by exp_mu.
const X = Array.from({length: n}, () => (rng() - 0.5) * 2);
const Y = X.map(xi => {
const mu = trueMu(xi);
const ref = Math.abs(mu[2]) < 0.9 ? [0,0,1] : [1,0,0];
const e1 = unit(cross(ref,mu)), e2 = cross(mu,e1);
const z1 = rnorm(), z2 = rnorm();
const eps = e1.map((z,j) => sigma*(z1*z + z2*e2[j]));
return expS2(mu, eps);
});
// Four data-driven parameters: two for p, one tangent direction, and speed.
function decode(q) {
const th=q[0], ph=q[1], al=q[2], speed=q[3];
const p=[Math.sin(th)*Math.cos(ph),Math.sin(th)*Math.sin(ph),Math.cos(th)];
const et=[Math.cos(th)*Math.cos(ph),Math.cos(th)*Math.sin(ph),-Math.sin(th)];
const ep=[-Math.sin(ph),Math.cos(ph),0];
return {p, v:et.map((z,j) => speed*(Math.cos(al)*z+Math.sin(al)*ep[j]))};
}
const mx=X.reduce((a,b)=>a+b,0)/n, vx=X.reduce((s,x)=>s+(x-mx)**2,0);
const my=[0,1,2].map(j=>Y.reduce((s,y)=>s+y[j],0)/n);
const beta=[0,1,2].map(j=>Y.reduce((s,y,i)=>s+(X[i]-mx)*(y[j]-my[j]),0)/vx);
const intercept=my.map((z,j)=>z-beta[j]*mx);
const p0=unit(intercept);
const vt=beta.map((z,j)=>z-dot(beta,p0)*p0[j]), speed0=Math.max(0.05,norm(vt));
const th0=Math.acos(Math.max(-1,Math.min(1,p0[2]))), ph0=Math.atan2(p0[1],p0[0]);
const et0=[Math.cos(th0)*Math.cos(ph0),Math.cos(th0)*Math.sin(ph0),-Math.sin(th0)];
const ep0=[-Math.sin(ph0),Math.cos(ph0),0];
const al0=Math.atan2(dot(vt,ep0),dot(vt,et0));
const objective=q => {
const z=decode(q);
return X.reduce((s,x,i)=>s+sphereDist(Y[i],expS2(z.p,z.v.map(a=>a*x)))**2,0);
};
const fit=nelderMead18(objective,[th0,ph0,al0,speed0],[0.08,0.08,0.12,0.08],220);
const fitted=decode(fit.q);
// Intrinsic R² relative to the sample Fréchet mean of the Y_i (gradient iteration on S²)
let fm = unit(my);
for (let it = 0; it < 60; it++) {
const w = [0, 0, 0];
for (const y of Y) {
const d = sphereDist(fm, y);
if (d < 1e-10) continue;
const f = d / Math.sin(d), c = dot(fm, y); // log_fm(y) = (d/sin d)(y − cos d · fm)
for (let j = 0; j < 3; j++) w[j] += (y[j] - c*fm[j]) * f / n;
}
if (norm(w) < 1e-10) break;
fm = unit(expS2(fm, w));
}
const sst = Y.reduce((s, y) => s + sphereDist(y, fm)**2, 0);
const r2 = sst > 0 ? 1 - fit.value/sst : 0;
const geodPts = Array.from({length:81},(_,j) => {
const t=-1+2*j/80;
return {x:t,true:trueMu(t),fitted:expS2(fitted.p,fitted.v.map(z=>z*t))};
});
return {X,Y,geodPts,bestV:fitted.v,vTrue:def.geodesic?def.v:null,sigma,n,sse:fit.value,r2,curveLabel:def.label,isGeodesic:def.geodesic};
}
grRes = runGeodReg(n18, noise18, curve18_def, 42);
{
// ---- 3D Globe Renderer for Geodesic Regression on S² ----
function initGlobe(cv, data, LW, LH) {
const ctx = cv.getContext('2d');
const dpr = cv.width / LW; // backing-store scale for HiDPI displays
// --- State ---
let rotY = -0.55, rotX = 0.40;
let velY = 0, velX = 0;
let dragging = false;
let dragStartX, dragStartY, dragStartRotY, dragStartRotX;
let lastDragX, lastDragY, lastDragTime;
let lastTime = performance.now();
let animId = null;
// --- Layout ---
const cx = 218, cy = 218, R = 168;
const camDist = 2.8;
// --- 3D rotation (Y then X) ---
function rotate(p, rx, ry) {
let [x, z] = [p[0], p[2]];
const c1 = Math.cos(ry), s1 = Math.sin(ry);
[x, z] = [c1 * x + s1 * z, -s1 * x + c1 * z];
let y = p[1];
const c2 = Math.cos(rx), s2 = Math.sin(rx);
[y, z] = [c2 * y - s2 * z, s2 * y + c2 * z];
return [x, y, z];
}
// --- Perspective projection: camera at (0,0,+camDist) looking at origin ---
// (front = +z: larger scale and full opacity; back = -z: smaller and faded)
function project(p) {
const [rx, ry, rz] = rotate(p, rotX, rotY);
const s = camDist / (camDist - rz);
return { x: cx + s * rx * R, y: cy - s * ry * R, z: rz, s };
}
// --- Sphere grid generators ---
function spherePt(lat, lon) {
const cl = Math.cos(lat);
return [cl * Math.cos(lon), Math.sin(lat), cl * Math.sin(lon)];
}
function makeParallel(lat, n) {
const pts = [];
for (let i = 0; i <= n; i++) pts.push(spherePt(lat, 2 * Math.PI * i / n));
return pts;
}
function makeMeridian(lon, n) {
const pts = [];
for (let i = 0; i <= n; i++) pts.push(spherePt(-Math.PI / 2 + Math.PI * i / n, lon));
return pts;
}
const N = 128;
const parallels = [-3, -2, -1, 0, 1, 2, 3].map(i => makeParallel(i * Math.PI / 8, N));
const meridians = Array.from({length: 16}, (_, i) => makeMeridian(i * Math.PI / 8, N));
// --- Draw a 3D curve with front/back depth cueing ---
function drawCurve(pts3D, color, width, dashPat) {
const proj = pts3D.map(project);
ctx.lineWidth = width;
ctx.lineCap = ctx.lineJoin = 'round';
for (let pass = 0; pass < 2; pass++) {
// pass 0 = back (faded), pass 1 = front (full color)
ctx.beginPath();
let started = false;
ctx.setLineDash(pass === 1 && dashPat ? dashPat : []);
ctx.strokeStyle = pass === 0
? color.replace('rgb', 'rgba').replace(')', ',0.22)')
: color;
for (const p of proj) {
const front = p.z > -0.04;
if ((pass === 0 && front) || (pass === 1 && !front)) { started = false; continue; }
if (!started) { ctx.moveTo(p.x, p.y); started = true; }
else ctx.lineTo(p.x, p.y);
}
if (started) ctx.stroke();
}
ctx.setLineDash([]);
}
// --- Draw data points, depth-sorted back-to-front ---
function drawPts(pts3D, color, r) {
const proj = pts3D.map(p => ({...project(p)})).sort((a, b) => a.z - b.z);
for (const p of proj) {
const alpha = 0.25 + 0.75 * Math.max(0, Math.min(1, (p.z + 1) / 2));
ctx.fillStyle = color.replace(/[\d.]+\)$/, alpha.toFixed(2) + ')');
ctx.beginPath();
ctx.arc(p.x, p.y, r * (0.6 + 0.4 * p.s), 0, 2 * Math.PI);
ctx.fill();
}
}
// --- Main render ---
function render() {
ctx.setTransform(dpr, 0, 0, dpr, 0, 0);
ctx.clearRect(0, 0, LW, LH);
// Background
ctx.fillStyle = '#f8f9fa';
ctx.fillRect(0, 0, LW, LH);
// Sphere halo / soft shadow. sr must be the silhouette radius of the
// perspective-projected sphere, R·d/√(d²−1), attained at cos θ = 1/d —
// otherwise grid lines and curves spill outside the shaded disc.
const sr = R * camDist / Math.sqrt(camDist * camDist - 1);
const grad1 = ctx.createRadialGradient(cx - 4, cy - 6, sr * 0.65, cx, cy, sr * 1.05);
grad1.addColorStop(0, '#e9ecef');
grad1.addColorStop(0.85, '#dee2e6');
grad1.addColorStop(1, '#ced4da');
ctx.fillStyle = grad1;
ctx.beginPath(); ctx.arc(cx, cy, sr * 1.03, 0, 2 * Math.PI); ctx.fill();
// Sphere body: radial gradient for 3D volume
const grad2 = ctx.createRadialGradient(cx - 0.12 * sr, cy - 0.22 * sr, sr * 0.06, cx, cy, sr);
grad2.addColorStop(0, 'rgba(255,255,255,0.94)');
grad2.addColorStop(0.40, 'rgba(222,235,247,0.70)');
grad2.addColorStop(1, 'rgba(165,192,218,0.48)');
ctx.fillStyle = grad2;
ctx.beginPath(); ctx.arc(cx, cy, sr, 0, 2 * Math.PI); ctx.fill();
// Wireframe grid
for (const p of parallels) drawCurve(p, 'rgb(173,181,189)', 0.55);
for (const m of meridians) drawCurve(m, 'rgb(173,181,189)', 0.50);
// Emphasize equator
drawCurve(makeParallel(0, 200), 'rgb(108,117,125)', 0.85);
// True geodesic (green dashed) and fitted geodesic (red solid)
drawCurve(data.geodPts.map(g => g.true), 'rgb(43,138,62)', 2.5, [7, 4]);
drawCurve(data.geodPts.map(g => g.fitted), 'rgb(224,49,49)', 2.7);
// Data points
drawPts(data.Y, 'rgba(116,192,252,0.72)', 3.6);
// North / South pole markers
for (const [py, label] of [[1, 'N'], [-1, 'S']]) {
const pp = project([0, py, 0]);
if (pp.z > -0.35) {
const a = pp.z > 0 ? 1 : 0.28;
ctx.fillStyle = `rgba(73,80,87,${a})`;
ctx.beginPath(); ctx.arc(pp.x, pp.y, 4.2, 0, 2 * Math.PI); ctx.fill();
ctx.fillStyle = `rgba(255,255,255,${a})`;
ctx.beginPath(); ctx.arc(pp.x, pp.y, 1.6, 0, 2 * Math.PI); ctx.fill();
ctx.fillStyle = `rgba(73,80,87,${a})`;
ctx.font = 'bold 11px system-ui, sans-serif';
ctx.fillText(label, pp.x + 7, pp.y + (py > 0 ? -5 : 14));
}
}
// --- Legend ---
const lx = 432, ly = 42;
ctx.font = '12px system-ui, sans-serif';
const subs = [];
if (data.vTrue) subs.push(`‖v_true‖ = ${Math.hypot(...data.vTrue).toFixed(3)} rad/unit`);
subs.push(`‖v_fit‖ = ${Math.hypot(...data.bestV).toFixed(3)} rad/unit`);
subs.push(`intrinsic R² = ${data.r2.toFixed(3)}`);
const items = [
{c: 'rgb(43,138,62)', d: [7, 4], t: data.isGeodesic ? 'True geodesic' : 'True curve', y: 0},
{c: 'rgb(224,49,49)', d: [], t: 'Fitted geodesic', y: 19},
{c: 'rgb(116,192,252)', t: 'Observed Yᵢ', y: 38, dot: 4.8},
...subs.map((t, k) => ({t, y: 65 + 18*k, sub: true})),
];
for (const it of items) {
const iy = ly + it.y;
if (it.dot) {
ctx.fillStyle = it.c;
ctx.beginPath(); ctx.arc(lx + 9, iy - 3, it.dot, 0, 2 * Math.PI); ctx.fill();
ctx.fillStyle = '#212529'; ctx.fillText(it.t, lx + 22, iy + 1);
} else if (it.c) {
ctx.strokeStyle = it.c; ctx.lineWidth = 2.2; ctx.setLineDash(it.d);
ctx.beginPath(); ctx.moveTo(lx, iy); ctx.lineTo(lx + 24, iy); ctx.stroke();
ctx.setLineDash([]);
ctx.fillStyle = '#212529'; ctx.fillText(it.t, lx + 32, iy + 4);
} else if (it.sub) {
ctx.fillStyle = '#495057'; ctx.font = '11px system-ui, sans-serif';
ctx.fillText(it.t, lx + 6, iy + 2);
ctx.font = '12px system-ui, sans-serif';
}
}
}
// --- Animation loop ---
function animate(ts) {
if (!document.body.contains(cv)) { animId = null; return; }
const dt = Math.min((ts - lastTime) / 1000, 0.1);
lastTime = ts;
if (!dragging && (Math.abs(velX) > 0.005 || Math.abs(velY) > 0.005)) {
rotX += velX * dt;
rotY += velY * dt;
rotX = Math.max(-Math.PI / 2.08, Math.min(Math.PI / 2.08, rotX));
velX *= Math.exp(-3.8 * dt);
velY *= Math.exp(-3.8 * dt);
render();
animId = requestAnimationFrame(animate);
} else if (!dragging) {
velX = 0; velY = 0;
animId = null;
} else {
render();
animId = requestAnimationFrame(animate);
}
}
function startAnim() {
if (!animId) { lastTime = performance.now(); animId = requestAnimationFrame(animate); }
}
// --- Mouse & touch input ---
function getPos(e) {
const r = cv.getBoundingClientRect();
const sx = cv.width / r.width, sy = cv.height / r.height;
const cx = e.touches ? e.touches[0].clientX : e.clientX;
const cy = e.touches ? e.touches[0].clientY : e.clientY;
return { x: (cx - r.left) * sx, y: (cy - r.top) * sy };
}
function onDown(e) {
dragging = true; velX = 0; velY = 0;
cv.style.cursor = 'grabbing';
const p = getPos(e);
dragStartX = p.x; dragStartY = p.y;
dragStartRotX = rotX; dragStartRotY = rotY;
lastDragX = p.x; lastDragY = p.y;
lastDragTime = performance.now();
// Track the drag even when the pointer leaves the canvas
window.addEventListener('mousemove', onMove);
window.addEventListener('mouseup', onUp);
startAnim();
e.preventDefault();
}
function onMove(e) {
if (!dragging) return;
const p = getPos(e);
const now = performance.now();
const ddt = Math.max(now - lastDragTime, 1) / 1000; // seconds
// Angular velocity in rad/s, with the same gain as the drag mapping below,
// so release momentum continues the drag motion seamlessly
velY = (p.x - lastDragX) / cv.width * Math.PI * 1.65 / ddt;
velX = (p.y - lastDragY) / cv.height * Math.PI * 1.65 / ddt;
lastDragX = p.x; lastDragY = p.y; lastDragTime = now;
rotY = dragStartRotY + (p.x - dragStartX) / cv.width * Math.PI * 1.65;
rotX = dragStartRotX + (p.y - dragStartY) / cv.height * Math.PI * 1.65;
rotX = Math.max(-Math.PI / 2.08, Math.min(Math.PI / 2.08, rotX));
e.preventDefault();
}
function onUp() {
if (!dragging) return;
dragging = false;
cv.style.cursor = 'grab';
window.removeEventListener('mousemove', onMove);
window.removeEventListener('mouseup', onUp);
// No fling if the pointer had already stopped moving before release
if (performance.now() - lastDragTime > 90) { velX = 0; velY = 0; }
if (Math.abs(velX) > 0.02 || Math.abs(velY) > 0.02) startAnim();
}
cv.addEventListener('mousedown', onDown);
cv.addEventListener('touchstart', onDown, {passive: false});
cv.addEventListener('touchmove', onMove, {passive: false});
cv.addEventListener('touchend', onUp);
// --- Initial render ---
render();
}
// Create the canvas element (HiDPI: backing store scaled by devicePixelRatio)
const W = 740, H = 430;
const dpr = Math.min(window.devicePixelRatio || 1, 2);
const canvas = document.createElement('canvas');
canvas.width = W * dpr;
canvas.height = H * dpr;
canvas.style.display = 'block';
canvas.style.width = '100%';
canvas.style.height = 'auto';
canvas.style.aspectRatio = `${W} / ${H}`;
canvas.style.border = '1px solid #dee2e6';
canvas.style.borderRadius = '4px';
canvas.style.cursor = 'grab';
canvas.style.background = '#f8f9fa';
canvas.style.touchAction = 'none';
initGlobe(canvas, grRes, W, H);
return html`
<div style="font-family:system-ui,sans-serif;max-width:760px;margin:0 auto;">
<h4>Geodesic Regression on S²</h4>
<div style="position:relative;">${canvas}</div>
<p style="font-size:0.85em;color:#868e96;margin-top:10px;">n=${grRes.n}, tangent-coordinate σ=${grRes.sigma} radians. True curve: <strong>${grRes.curveLabel}</strong>. <strong>Drag to rotate</strong> the 3D globe. The fit jointly estimates the intercept, tangent direction, and speed from the simulated observations; for a non-geodesic truth the fitted geodesic shows systematic lack of fit (lower intrinsic R²).</p>
</div>
`
}
```
## Application: Revisiting Market Risk with Geodesic Regression on SPD Manifolds
Lectures 4–8 developed a suite of regression and ANOVA tools for SPD-valued data using **metric-space** methods: Fréchet ANOVA, global Fréchet regression, kernel Fréchet regression, and local-linear Fréchet regression. All of these operate on the principle that, under the log-Euclidean metric $d_{\mathrm{LE}}(A, B) = \|\log A - \log B\|_F$, the weighted Fréchet mean has a closed form:
$$
\log \hat{\mu}(x) = \sum_{i=1}^n w_i(x)\, \log \Sigma_i, \qquad \hat{\mu}(x) = \exp\!\bigl(\log\hat{\mu}(x)\bigr).
$$
Within these log-Euclidean constructions, the main computational difference is **how the weights $w_i(x)$ are computed**. The response matrices then enter through squared distances and weighted averaging. This is flexible, but it does not introduce an explicit constant tangent slope or constrain the fitted trajectory to one geodesic.
This section revisits the market risk application from the **manifold perspective**. We now view $\mathcal{S}_{++}^p$ not merely as a metric space but as a **Riemannian symmetric space** (under the affine-invariant metric), and ask: what additional capabilities does the manifold structure unlock?
### From Weighted Averaging to Geodesic Modeling
In the metric-space approach, the regression function at stress level $x$ is a weighted Fréchet mean — a point on the manifold obtained by minimizing a weighted sum of squared distances. **There is no imposed parametric curve.** Fits at different target values share the same data and weight formula, but they are not constrained to lie on a single geodesic.
Geodesic regression takes the opposite approach. It posits that the conditional Fréchet mean follows a **single geodesic curve**:
$$
\mu(x) = \gamma_{p,v}(x) = \exp_p(x v), \qquad (p, v) \in T\mathcal{S}_{++}^p.
$$
This is a parametric model with a tangent-bundle parameter: the intercept $p \in \mathcal{S}_{++}^p$ (the typical covariance matrix at $x = 0$, which we center so that average stress corresponds to $x = 0$) and the slope $v \in T_p\mathcal{S}_{++}^p$ (a symmetric matrix describing the direction and rate of change of the covariance structure per unit stress).
### Geodesic Regression Under the Affine-Invariant Metric
Under the standard **affine-invariant metric** $d_{\mathrm{AI}}(A, B) = \|\log(A^{-1/2} B A^{-1/2})\|_F$, the exponential map is [@PennecFillardAyache2006]:
$$
\exp_P(V) = P^{1/2} \exp(P^{-1/2} V P^{-1/2}) P^{1/2},
$$
where $V \in T_P\mathcal{S}_{++}^p$ is a symmetric matrix (the tangent space at $P$ is the space of symmetric matrices). The geodesic from $P$ in direction $V$ is:
$$
\gamma_{P,V}(t) = P^{1/2} \exp(t\, P^{-1/2} V P^{-1/2}) P^{1/2}.
$$
Given observations $(\hat{\sigma}_i, \Sigma_i)$ of market stress and covariance matrices, the geodesic least-squares estimator solves:
$$
(\hat{p}, \hat{v}) = \arg\min_{(p,v) \in T\mathcal{S}_{++}^p} \frac{1}{2}\sum_{i=1}^n d_{\mathrm{AI}}^2\!\bigl(\Sigma_i, \gamma_{p,v}(\hat{\sigma}_i)\bigr).
$$
The key computational difference from the metric-space approach: this is a **joint optimization** over the tangent bundle. Minimization uses intrinsic gradient descent with gradients computed via the differentials of the exponential map (see the Gradients section above). There is no closed-form weighted-average solution.
::: {.callout-important title="The affine-invariant metric v.s. the log-Euclidean metric"}
The log-Euclidean metric made all three Fréchet regression methods computationally simple because $\mathcal{S}_{++}^p$ under $d_{\mathrm{LE}}$ is isometric to a Euclidean space of symmetric matrices. This metric is not invariant under arbitrary congruence transformations $\Sigma \mapsto A\Sigma A^\top$. The standard affine-invariant metric is congruence-invariant and is often natural for covariance matrices. It is not, however, the only congruence-invariant Riemannian metric: a broader invariant family exists [@ThanwerdasPennec2023]. Geodesic regression can be formulated for any chosen Riemannian metric, although the availability and cost of exponential maps, logarithms, and their differentials vary substantially.
:::
### Comparison: Metric-Space Regression vs. Geodesic Regression
| Aspect | Fréchet Regression (metric) | Geodesic Regression (manifold) |
|---|---|---|
| **Model** | $\mu(x) = \arg\min_\omega \mathbb{E}[s(X,x)d^2(Y,\omega)]$ | $\mu(x) = \exp_p(x v)$ |
| **Parameters** | Implicit (defined by weighted averaging) | Explicit: $(p, v) \in T\mathcal{M}$ |
| **Parameter interpretation** | Regression target defined implicitly; no single constant slope parameter | $p$ = intercept on manifold, $v$ = slope in $T_p\mathcal{M}$ |
| **Curve constraint** | None — each $x$ estimated independently | All $\mu(x)$ lie on a single geodesic |
| **Computational core** | Weighted Fréchet mean | Joint optimization over $T\mathcal{M}$ (iterative) |
| **Metric flexibility** | Broad metric-space definition; the SPD closed form used here is log-Euclidean | Any Riemannian metric in principle; computation is metric-specific |
| **MLE interpretation** | None in general | MLE under intrinsic Gaussian (symmetric spaces) |
| **Inference tools** | Depends on the model; bootstrap and asymptotics are common | Intrinsic $R^2$ and permutation test; likelihood tools under an explicit error model |
| **Generalization** | Any metric space | Riemannian manifolds with computable exp/log |
### What the Manifold Structure Adds
Beyond the comparison table, a geodesic parameterization makes several geometric operations directly available:
**1. Explicit tangent-space slope with geometric interpretation.** In geodesic regression, $\hat{v} \in T_{\hat{p}}\mathcal{S}_{++}^p$ is a symmetric matrix that encodes the *direction* in which the covariance structure changes with stress. Under the affine-invariant metric, the coordinate-free rate information is represented at the identity by $\hat{p}^{-1/2}\hat{v}\hat{p}^{-1/2}$; its eigendirections describe relative covariance-change modes. The raw Euclidean eigenvectors of $\hat v$ alone are basis-dependent. In contrast, the metric-space weights $s_{in}(x)$ are predictor-side quantities rather than an explicit constant response-space slope.
**2. Parallel transport for comparing slopes across groups.** Suppose we fit geodesic regression separately to data from three market regimes (calm, normal, stress). The estimated slopes $\hat{v}_{\text{calm}}, \hat{v}_{\text{normal}}, \hat{v}_{\text{stress}}$ live in different tangent spaces (at $\hat{p}_{\text{calm}}, \hat{p}_{\text{normal}}, \hat{p}_{\text{stress}}$, respectively). Parallel transport along the geodesic connecting the intercepts maps all slopes to a common tangent space, where they can be compared directly:
$$
\hat{v}_{\text{calm} \to \text{normal}} = \Gamma_{\hat{p}_{\text{calm}} \to \hat{p}_{\text{normal}}}(\hat{v}_{\text{calm}}).
$$
The norm $\|\hat{v}_{\text{calm} \to \text{normal}} - \hat{v}_{\text{normal}}\|$ is a descriptive effect size for the difference in structural change. Turning it into a statistical test additionally requires a sampling distribution, such as one obtained from an appropriate bootstrap or regression asymptotics.
**3. Intrinsic Gaussian MLE and likelihood-based inference.** On the symmetric space $\mathcal{S}_{++}^p$ (affine-invariant), minimizing geodesic least-squares is **equivalent** to maximum likelihood under the intrinsic Gaussian model:
$$
f(\Sigma \mid x; p, v, \sigma^2) \propto \exp\!\left(-\frac{d_{\mathrm{AI}}^2(\Sigma, \gamma_{p,v}(x))}{2\sigma^2}\right).
$$
The normalizing constant $C(\mu, \sigma)$ does **not** depend on $\mu$ on a symmetric space, so for fixed $\sigma$ the log-likelihood for $(p,v)$ differs from the negative sum of squared distances only by constants [@Fletcher2013]. This supplies a likelihood interpretation for least squares and may enable likelihood-ratio tests or confidence regions.
### Interactive Exploration: Geodesic vs. Global Fréchet Regression on Simulated SPD Data
The following demo compares **geodesic regression** (affine-invariant metric) with **global Fréchet regression** (log-Euclidean metric) on simulated $2 \times 2$ SPD matrices. The true relationship is a geodesic under the affine-invariant metric; the global Fréchet fit uses the log-Euclidean shortcut.
```{ojs}
//| label: fig-geod-vs-frechet
//| fig-cap: "Interactive: Geodesic regression (affine-invariant) vs. Global Fréchet regression (log-Euclidean) on 2×2 SPD matrices"
n_app18_control = Inputs.range([15, 100], {step: 5, value: 40, label: "Sample size n"})
noise_app18_control = Inputs.range([0.05, 0.6], {step: 0.05, value: 0.2, label: "Tangent noise σ"})
seed_app18_control = Inputs.range([1, 50], {step: 1, value: 7, label: "Random seed"})
n_app18 = Generators.input(n_app18_control)
noise_app18 = Generators.input(noise_app18_control)
seed_app18 = Generators.input(seed_app18_control)
geospd_controls_view = html`
<style>
.geospd-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; }
.geospd-slider-grid > * { flex:1 1 calc((100% - 40px)/3); min-width:0; margin:0; }
.geospd-slider-grid input[type="number"] { width:7.5rem !important; }
@container (max-width:700px) { .geospd-slider-grid > * { flex-basis:calc((100% - 20px)/2); } }
@container (max-width:480px) { .geospd-slider-grid > * { flex-basis:100%; } }
</style>
<div class="geospd-slider-grid">
<div>${n_app18_control}</div>
<div>${noise_app18_control}</div>
<div>${seed_app18_control}</div>
</div>`
// ---- SPD 2x2 utilities (affine-invariant) ----
function matLog(m) {
// Stable spectral formula, including diagonal and repeated-eigenvalue cases.
const mid=(m[0]+m[2])/2, h=(m[0]-m[2])/2;
const delta=Math.hypot(h,m[1]);
const lo=Math.max(1e-15,mid-delta), hi=Math.max(1e-15,mid+delta);
if (delta < 1e-12*Math.max(1,Math.abs(mid))) {
const z=Math.log(Math.max(mid,1e-15));
return [z,0,z];
}
const meanLog=(Math.log(lo)+Math.log(hi))/2;
const scale=(Math.log(hi)-Math.log(lo))/(2*delta);
return [meanLog+scale*h,scale*m[1],meanLog-scale*h];
}
function matExp(m) {
const mid=(m[0]+m[2])/2, h=(m[0]-m[2])/2;
const delta=Math.hypot(h,m[1]), em=Math.exp(mid);
const shc=delta < 1e-10 ? 1+delta*delta/6 : Math.sinh(delta)/delta;
const ch=Math.cosh(delta);
return [em*(ch+shc*h),em*shc*m[1],em*(ch-shc*h)];
}
function aiDist(a, b) {
// Eigenvalues of a^{-1}b equal those of a^{-1/2}ba^{-1/2}.
const detA = a[0]*a[2] - a[1]*a[1];
if (detA <= 0) return Infinity;
const aInv0 = a[2]/detA, aInv1 = -a[1]/detA, aInv2 = a[0]/detA;
// compute a^{-1/2} b a^{-1/2}
// For 2x2, use eigenvalue decomposition of aInv * b
const m00 = aInv0*b[0] + aInv1*b[1];
const m01 = aInv0*b[1] + aInv1*b[2];
const m10 = aInv1*b[0] + aInv2*b[1];
const m11 = aInv1*b[1] + aInv2*b[2];
// eigenvalues of [[m00,m01],[m10,m11]]
const trM = m00 + m11, detM = m00*m11 - m01*m10;
const discM = Math.sqrt(Math.max(0, trM*trM - 4*detM));
const l1 = Math.max(1e-15, (trM - discM)/2), l2 = Math.max(1e-15, (trM + discM)/2);
return Math.sqrt(Math.log(l1)**2 + Math.log(l2)**2);
}
function aiGeodesic(p, v, t) {
// gamma(t) = p^{1/2} exp(t * p^{-1/2} v p^{-1/2}) p^{1/2}
const pSqrt = matSqrt(p);
const pInvSqrt = matInvSqrt(p);
// Compute p^{-1/2} v p^{-1/2}
const vTransformed = matCongruence(pInvSqrt, v);
const expTV = matExp([t*vTransformed[0], t*vTransformed[1], t*vTransformed[2]]);
return matCongruence(pSqrt, expTV);
}
function matSqrt(m) {
return matExp([0.5*matLog(m)[0], 0.5*matLog(m)[1], 0.5*matLog(m)[2]]);
}
function matInvSqrt(m) {
const s = matLog(m);
return matExp([-0.5*s[0], -0.5*s[1], -0.5*s[2]]);
}
function matCongruence(a, x) {
// a * x * a^T (a, x symmetric 2x2)
const a00 = a[0], a01 = a[1], a10 = a[1], a11 = a[2];
const ax00 = a00*x[0] + a01*x[1], ax01 = a00*x[1] + a01*x[2];
const ax10 = a10*x[0] + a11*x[1], ax11 = a10*x[1] + a11*x[2];
return [
ax00*a00 + ax01*a10,
ax00*a01 + ax01*a11,
ax10*a01 + ax11*a11
];
}
function nelderMeadSPD18(f, start, steps, maxIter) {
const d=start.length;
let simplex=[start.slice()];
for (let j=0;j<d;j++) { const q=start.slice(); q[j]+=steps[j]; simplex.push(q); }
let vals=simplex.map(f);
for (let iter=0;iter<maxIter;iter++) {
const ord=Array.from({length:d+1},(_,i)=>i).sort((i,j)=>vals[i]-vals[j]);
simplex=ord.map(i=>simplex[i]); vals=ord.map(i=>vals[i]);
const avg=vals.reduce((a,b)=>a+b,0)/(d+1);
const spread=Math.sqrt(vals.reduce((s,z)=>s+(z-avg)**2,0)/(d+1));
if (spread < 1e-9*(1+Math.abs(avg))) break;
const cen=Array(d).fill(0);
for (let i=0;i<d;i++) for (let j=0;j<d;j++) cen[j]+=simplex[i][j]/d;
const trial=k=>cen.map((z,j)=>z+k*(z-simplex[d][j]));
const qr=trial(1), fr=f(qr);
if (fr<vals[0]) {
const qe=trial(2), fe=f(qe);
if (fe<fr) { simplex[d]=qe; vals[d]=fe; } else { simplex[d]=qr; vals[d]=fr; }
} else if (fr<vals[d-1]) {
simplex[d]=qr; vals[d]=fr;
} else {
const qc=fr<vals[d] ? trial(0.5) : cen.map((z,j)=>z+0.5*(simplex[d][j]-z));
const fc=f(qc);
if (fc<Math.min(fr,vals[d])) { simplex[d]=qc; vals[d]=fc; }
else {
for (let i=1;i<=d;i++) {
simplex[i]=simplex[0].map((z,j)=>z+0.5*(simplex[i][j]-z));
vals[i]=f(simplex[i]);
}
}
}
}
const best=vals.indexOf(Math.min(...vals));
return {q:simplex[best],value:vals[best]};
}
// ---- Simulation ----
function runGeoVsFrechet(n, noise, 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 || 42);
function rnorm() {
return Math.sqrt(-2*Math.log(Math.max(rng(),1e-15)))*Math.cos(2*Math.PI*rng());
}
const pTrue = [1.5,0.5,2.0];
const pTrueSqrt=matSqrt(pTrue);
// B is the slope in identity coordinates; V=P^{1/2}BP^{1/2}.
const bTrue=[0.15,0.32,0.50];
const vTrue=matCongruence(pTrueSqrt,bTrue);
const X = Array.from({length: n}, () => (rng() - 0.5) * 3);
const Y = X.map(xi => {
const mu = aiGeodesic(pTrue, vTrue, xi);
const s=[noise*rnorm(),noise*rnorm()/Math.sqrt(2),noise*rnorm()];
const eps=matCongruence(matSqrt(mu),s);
return aiGeodesic(mu,eps,1);
});
const logY = Y.map(y => matLog(y));
const meanX = X.reduce((a,b) => a+b, 0)/n;
const varX = X.reduce((s,xi) => s + (xi-meanX)**2, 0)/n;
const meanLog=[0,1,2].map(j=>logY.reduce((s,z)=>s+z[j],0)/n);
const leSlope=[0,1,2].map(j=>X.reduce((s,x,i)=>s+(x-meanX)*(logY[i][j]-meanLog[j]),0)/(n*varX));
const leIntercept=meanLog.map((z,j)=>z-meanX*leSlope[j]);
function globalFrechetFit(x) {
return matExp(leIntercept.map((z,j)=>z+x*leSlope[j]));
}
// Parameterize P=exp(S), V=P^{1/2}BP^{1/2}; all trial P remain SPD.
function decode(q) {
const s=q.slice(0,3), b=q.slice(3,6), p=matExp(s), ps=matExp(s.map(z=>z/2));
return {p,v:matCongruence(ps,b)};
}
const objective=q => {
const z=decode(q);
const value=X.reduce((s,x,i)=>s+aiDist(Y[i],aiGeodesic(z.p,z.v,x))**2,0);
return Number.isFinite(value) ? value : 1e100;
};
// Affine-invariant sample Fréchet mean for the constant-model baseline.
let frechetMean=matExp(meanLog);
for (let iter=0;iter<60;iter++) {
const invs=matInvSqrt(frechetMean), avg=[0,0,0];
for (const y of Y) {
const r=matLog(matCongruence(invs,y));
for (let j=0;j<3;j++) avg[j]+=r[j]/n;
}
const rnorm=Math.sqrt(avg[0]**2+2*avg[1]**2+avg[2]**2);
if (rnorm<1e-10) break;
frechetMean=matCongruence(matSqrt(frechetMean),matExp(avg));
}
const totalSSE=Y.reduce((s,y)=>s+aiDist(y,frechetMean)**2,0);
const constantStart=[...matLog(frechetMean),0,0,0];
const optLE=nelderMeadSPD18(objective,[...leIntercept,...leSlope],[0.08,0.05,0.08,0.08,0.05,0.08],260);
const opt0=nelderMeadSPD18(objective,constantStart,[0.05,0.03,0.05,0.05,0.03,0.05],260);
const opt=[optLE,opt0,{q:constantStart,value:totalSSE}].sort((a,b)=>a.value-b.value)[0];
const geoHat=decode(opt.q), pGeoHat=geoHat.p, vGeoHat=geoHat.v, bestGeoSSE=opt.value;
// ---- Predictions on a grid ----
const xGrid = [];
const nGrid = 40;
for (let j = 0; j <= nGrid; j++) xGrid.push(-1.5 + j*3/nGrid);
const trueCurve = xGrid.map(x => aiGeodesic(pTrue, vTrue, x));
const geoCurve = xGrid.map(x => aiGeodesic(pGeoHat, vGeoHat, x));
const frechetCurve = xGrid.map(x => globalFrechetFit(x));
const geoR2=totalSSE>0 ? 1-bestGeoSSE/totalSSE : 0;
const frechetSSE=X.reduce((s,x,i)=>s+aiDist(Y[i],globalFrechetFit(x))**2,0);
// Same AI-distance improvement score, but not an intrinsic R² for the LE fit.
const frechetR2=totalSSE>0 ? 1-frechetSSE/totalSSE : 0;
return { n, noise, X, Y, xGrid, trueCurve, geoCurve, frechetCurve, pTrue, vTrue, pGeoHat, vGeoHat, geoR2, frechetR2, bestGeoSSE };
}
geoVsFrechet = runGeoVsFrechet(n_app18, noise_app18, seed_app18);
html`
<div style="font-family: system-ui, sans-serif; max-width: 850px;">
<h4>Geodesic (Affine-Invariant) vs. Global Fréchet (Log-Euclidean) Regression</h4>
<p style="font-size:0.85em;color:#495057;">Showing determinant and off-diagonal entry of fitted SPD matrices over stress range</p>
<svg viewBox="0 0 820 380" style="display:block;width:100%;height:auto;border:1px solid #dee2e6;border-radius:4px;">
${(() => {
const mg = {top: 20, right: 150, bottom: 35, left: 55};
const pW = 820 - mg.left - mg.right, pH = 380 - mg.top - mg.bottom;
const xMin = -1.5, xMax = 1.5;
const xS = x => mg.left + (x - xMin)/(xMax - xMin) * pW;
const g = geoVsFrechet;
const det = m => m[0]*m[2] - m[1]*m[1];
const paddedRange = values => {
const lo=Math.min(...values), hi=Math.max(...values), pad=0.08*Math.max(hi-lo,1e-3);
return [lo-pad,hi+pad];
};
const detRange=paddedRange([...g.Y,...g.trueCurve,...g.geoCurve,...g.frechetCurve].map(det));
const offRange=paddedRange([...g.Y,...g.trueCurve,...g.geoCurve,...g.frechetCurve].map(m=>m[1]));
const ticks=range=>Array.from({length:5},(_,j)=>range[0]+j*(range[1]-range[0])/4);
// Panel A (top): determinant
const yMinA=detRange[0], yMaxA=detRange[1];
const ySA = y => mg.top + 140 - (y - yMinA)/(yMaxA - yMinA) * 120;
// Panel B (bottom): off-diagonal
const yMinB=offRange[0], yMaxB=offRange[1];
const ySB = y => mg.top + 325 - (y - yMinB)/(yMaxB - yMinB) * 120;
// Panel A lines: determinant
const trueA = g.xGrid.map((x,i) => `${i===0?'M':'L'} ${xS(x).toFixed(1)} ${ySA(det(g.trueCurve[i])).toFixed(1)}`).join(' ');
const geoA = g.xGrid.map((x,i) => `${i===0?'M':'L'} ${xS(x).toFixed(1)} ${ySA(det(g.geoCurve[i])).toFixed(1)}`).join(' ');
const frecA = g.xGrid.map((x,i) => `${i===0?'M':'L'} ${xS(x).toFixed(1)} ${ySA(det(g.frechetCurve[i])).toFixed(1)}`).join(' ');
// Panel B lines: off-diagonal
const trueB = g.xGrid.map((x,i) => `${i===0?'M':'L'} ${xS(x).toFixed(1)} ${ySB(g.trueCurve[i][1]).toFixed(1)}`).join(' ');
const geoB = g.xGrid.map((x,i) => `${i===0?'M':'L'} ${xS(x).toFixed(1)} ${ySB(g.geoCurve[i][1]).toFixed(1)}`).join(' ');
const frecB = g.xGrid.map((x,i) => `${i===0?'M':'L'} ${xS(x).toFixed(1)} ${ySB(g.frechetCurve[i][1]).toFixed(1)}`).join(' ');
// Scatter: observed determinants
const scatA = g.X.map((xi,i) => `<circle cx="${xS(xi).toFixed(1)}" cy="${ySA(det(g.Y[i])).toFixed(1)}" r="2.5" fill="#adb5bd" opacity="0.45"/>`).join('');
const scatB = g.X.map((xi,i) => `<circle cx="${xS(xi).toFixed(1)}" cy="${ySB(g.Y[i][1]).toFixed(1)}" r="2.5" fill="#adb5bd" opacity="0.45"/>`).join('');
return `
<!-- Panel A background -->
<rect x="${mg.left}" y="${mg.top}" width="${pW}" height="150" fill="#f8f9fa" rx="2"/>
<text x="${mg.left+10}" y="${mg.top+14}" font-size="10" fill="#495057" font-weight="bold">Panel A: Determinant |Σ|</text>
<!-- axes -->
<line x1="${mg.left}" y1="${mg.top+150}" x2="${mg.left+pW}" y2="${mg.top+150}" stroke="#adb5bd" stroke-width="0.5"/>
<line x1="${mg.left}" y1="${mg.top}" x2="${mg.left}" y2="${mg.top+150}" stroke="#adb5bd" stroke-width="0.5"/>
<text x="${mg.left+pW/2}" y="${mg.top+167}" text-anchor="middle" font-size="9">Stress (x)</text>
<!-- y ticks A -->
${ticks(detRange).map(y =>
`<text x="${mg.left-5}" y="${ySA(y)+4}" text-anchor="end" font-size="8">${y.toFixed(1)}</text>` +
`<line x1="${mg.left}" y1="${ySA(y)}" x2="${mg.left+pW}" y2="${ySA(y)}" stroke="#e9ecef" stroke-width="0.5"/>`
).join('')}
${scatA}
<path d="${trueA}" fill="none" stroke="#2b8a3e" stroke-width="1.8" stroke-dasharray="6,3"/>
<path d="${geoA}" fill="none" stroke="#1971c2" stroke-width="2.2"/>
<path d="${frecA}" fill="none" stroke="#e03131" stroke-width="2"/>
<!-- Panel B background -->
<rect x="${mg.left}" y="${mg.top+185}" width="${pW}" height="150" fill="#f8f9fa" rx="2"/>
<text x="${mg.left+10}" y="${mg.top+199}" font-size="10" fill="#495057" font-weight="bold">Panel B: Off-diagonal Σ₁₂</text>
<line x1="${mg.left}" y1="${mg.top+335}" x2="${mg.left+pW}" y2="${mg.top+335}" stroke="#adb5bd" stroke-width="0.5"/>
<line x1="${mg.left}" y1="${mg.top+185}" x2="${mg.left}" y2="${mg.top+335}" stroke="#adb5bd" stroke-width="0.5"/>
<text x="${mg.left+pW/2}" y="${mg.top+352}" text-anchor="middle" font-size="9">Stress (x)</text>
${ticks(offRange).map(y =>
`<text x="${mg.left-5}" y="${ySB(y)+4}" text-anchor="end" font-size="8">${y.toFixed(1)}</text>` +
`<line x1="${mg.left}" y1="${ySB(y)}" x2="${mg.left+pW}" y2="${ySB(y)}" stroke="#e9ecef" stroke-width="0.5"/>`
).join('')}
${scatB}
<path d="${trueB}" fill="none" stroke="#2b8a3e" stroke-width="1.8" stroke-dasharray="6,3"/>
<path d="${geoB}" fill="none" stroke="#1971c2" stroke-width="2.2"/>
<path d="${frecB}" fill="none" stroke="#e03131" stroke-width="2"/>
<!-- Legend -->
<g transform="translate(${mg.left+pW+5},${mg.top+15})">
<line x1="0" y1="4" x2="16" y2="4" stroke="#2b8a3e" stroke-width="1.8" stroke-dasharray="6,3"/><text x="20" y="8" font-size="9">True geodesic (AI)</text>
<line x1="0" y1="20" x2="16" y2="20" stroke="#1971c2" stroke-width="2.2"/><text x="20" y="24" font-size="9">Geodesic reg. (AI)</text>
<line x1="0" y1="36" x2="16" y2="36" stroke="#e03131" stroke-width="2"/><text x="20" y="40" font-size="9">Global Fréchet (LE)</text>
<circle cx="4" cy="52" r="2.5" fill="#adb5bd"/><text x="20" y="56" font-size="9">Observed Σᵢ</text>
</g>
<!-- R² display -->
<g transform="translate(${mg.left+pW+5},${mg.top+240})">
<text x="0" y="0" font-size="10" font-weight="bold">AI-distance improvement:</text>
<text x="0" y="16" font-size="10"> Geodesic: ${g.geoR2.toFixed(3)}</text>
<text x="0" y="32" font-size="10"> Fréchet: ${g.frechetR2.toFixed(3)}</text>
<text x="0" y="52" font-size="9" fill="#868e96">Relative to the AI Fréchet mean</text>
</g>
`;
})()}
</svg>
<p style="font-size:0.85em;color:#868e96;">
n = ${geoVsFrechet.n}, tangent-coordinate σ = ${geoVsFrechet.noise}.
The true data-generating process is a geodesic under the affine-invariant (AI) metric.
The global Fréchet fit uses the log-Euclidean (LE) shortcut — a straight line in the log domain — which is a <em>different</em> curve on the manifold.
The geodesic value is intrinsic R². The LE value is the analogous AI-distance improvement score and can be negative because that fit does not minimize AI loss.
</p>
</div>
`
```
**How to read the demo:** Because the data-generating mean is an affine-invariant geodesic, the correctly specified geodesic fit will generally have lower affine-invariant loss as $n$ grows. A finite noisy sample can occasionally favor the log-Euclidean fit. Both geometries can represent changing eigenvectors; the curves differ here because the intercept and direction do not commute, so flattening with the matrix logarithm does not preserve the affine-invariant geodesic.
## Key Takeaways
- **Geodesic regression** fits a geodesic curve $\exp_p(x v)$ to manifold-valued responses — the direct analogue of simple linear regression.
- **Gradients** are computed via differentials of the exponential map; estimation uses intrinsic gradient descent.
- **Intrinsic $R^2$** generalizes the Euclidean coefficient of determination; permutation tests assess significance.
- Under Fletcher's curvature and concentration conditions, **symmetric spaces** give a unique minimizer within the specified class of $r$-close geodesics.
- Under the associated conditions and intrinsic Gaussian error model, geodesic least-squares for $(p,v)$ agrees with maximum likelihood at fixed scale.
- **Metric-space vs. manifold regression:** the Fréchet methods used here define targets through predictor-side weights and distance minimization; geodesic regression imposes a parametric curve *on* the manifold, giving an explicit slope, transport-based comparisons, and a likelihood interpretation under Fletcher's error model.
## Exercises
1. **Euclidean recovery.** Show that on $\mathbb{R}^m$, geodesic regression reduces to ordinary least squares. <a href="javascript:void(0)" onclick="showSolution('l18-sol-1')" class="solution-link">📝 Show Solution</a>
2. **SPD inversion symmetry.** Verify that $\sigma_P(\Sigma) = P\Sigma^{-1}P$ is an involutive isometry for the affine-invariant metric. <a href="javascript:void(0)" onclick="showSolution('l18-sol-2')" class="solution-link">📝 Show Solution</a>
3. **Intrinsic $R^2$ interpretation.** What does $R^2 = 0$ and $R^2 = 1$ mean in geodesic regression? <a href="javascript:void(0)" onclick="showSolution('l18-sol-3')" class="solution-link">📝 Show Solution</a>
4. **Gradient interpretation.** Explain why $\nabla_p E_n$ involves the adjoint of the differential of $\exp$. <a href="javascript:void(0)" onclick="showSolution('l18-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="l18-sol-1"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 1</h4><button onclick="closeSolution('l18-sol-1')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Show geodesic regression = OLS on R^m.</div><div class="solution-answer"><strong>Solution:</strong> On $\mathbb{R}^m$, $\exp_p(v) = p + v$. The geodesic is $\gamma_{p,v}(x) = p + x v$. The objective is $E_n(p,v) = \frac{1}{2}\sum \|y_i - p - x_i v\|^2$. Setting derivatives to zero: $\sum(y_i - p - x_i v) = 0$ and $\sum x_i(y_i - p - x_i v) = 0$. These are exactly the normal equations for simple linear regression. The minimizer is $\hat{v} = \sum(x_i - \bar{x})(y_i - \bar{y}) / \sum(x_i - \bar{x})^2$, $\hat{p} = \bar{y} - \hat{v}\bar{x}$.</div></div></dialog>
<dialog id="l18-sol-2"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 2</h4><button onclick="closeSolution('l18-sol-2')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Verify $\sigma_P$ is involutive isometry for affine-invariant metric.</div><div class="solution-answer"><strong>Solution:</strong> **Involutive:** $\sigma_P(\sigma_P(\Sigma))=P(P\Sigma^{-1}P)^{-1}P=\Sigma$. **Isometry:** The affine-invariant distance is unchanged by congruence, $d_{\mathrm{AI}}(CAC^\top,CBC^\top)=d_{\mathrm{AI}}(A,B)$, and by inversion, $d_{\mathrm{AI}}(A^{-1},B^{-1})=d_{\mathrm{AI}}(A,B)$. (For inversion, the generalized eigenvalues are replaced by their reciprocals, whose logarithms change only sign.) Since $\sigma_P(\Sigma)=P\Sigma^{-1}P$ is inversion followed by congruence with $C=P$, their composition is an isometry. Together with involutivity and $\sigma_P(P)=P$, this is the geodesic symmetry at $P$.</div></div></dialog>
<dialog id="l18-sol-3"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 3</h4><button onclick="closeSolution('l18-sol-3')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Interpret $R^2 = 0$ and $R^2 = 1$ in geodesic regression.</div><div class="solution-answer"><strong>Solution:</strong> $R^2=0$ means the fitted geodesic attains no reduction in squared-distance loss relative to the best constant fit through a sample Fréchet mean. A zero-slope fit realizes this value, but $R^2=0$ need not force $\hat v=0$ if minimizers are nonunique. $R^2=1$ means zero residual loss: $y_i=\exp_{\hat p}(x_i\hat v)$ for every $i$. For a global minimizer, $R^2$ measures the fraction of sample Fréchet variance removed by the geodesic fit.</div></div></dialog>
<dialog id="l18-sol-4"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 4</h4><button onclick="closeSolution('l18-sol-4')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Why the adjoint in the gradient formula?</div><div class="solution-answer"><strong>Solution:</strong> The residual $\epsilon_i = \log_{\mu_i}(y_i)$ lives in $T_{\mu_i}\mathcal{M}$. The gradient with respect to $p$ must live in $T_p\mathcal{M}$. The differential $D_1(\exp)_{(p, x_i v)}$ maps $T_p\mathcal{M} \to T_{\mu_i}\mathcal{M}$ (it says how $\exp_p(x_i v)$ changes when $p$ moves). Its adjoint $D_1(\exp)^*$ maps back: $T_{\mu_i}\mathcal{M} \to T_p\mathcal{M}$. This is needed because the chain rule gives $d_{p}[d^2(y_i, \exp_p(x_i v))] = -2\langle \epsilon_i, D_1(\exp)_{(p, x_i v)}[\delta p]\rangle_{\mu_i} = -2\langle D_1(\exp)_{(p, x_i v)}^* \epsilon_i, \delta p\rangle_p$. So $\nabla_p E_n = -\sum D_1(\exp)^* \epsilon_i$ as claimed.</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
- @Fletcher2013 — Geodesic regression on Riemannian manifolds.
- @ShinOh2021GeodRegr — R package for geodesic regression.
- @Lee2018 — Riemannian geometry textbook.
- @PennecFillardAyache2006 — Manifold statistics on SPD matrices.
- @ThanwerdasPennec2023 — Families and invariance properties of Riemannian metrics on SPD matrices.
## Self-Assessment Quiz
Test your understanding of this lecture with the interactive MCQ quiz:
👉 **[Lecture 18 Quiz — 10 Multiple-Choice Questions](../quizzes/lecture-18-quiz.qmd)**