---
title: "Lecture 10: Logistic Regression with Metric-Space Covariates"
subtitle: "Binary classification when predictors live in a nonlinear metric space"
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:
- Explain how the roles of predictor and response are reversed in metric-space logistic regression compared to Fréchet regression.
- Derive the Alexandrov inner-product representation of the usual Euclidean logistic regression and explain why it suggests a metric-space generalization.
- State the metric-space logistic model using the Fréchet mean anchor and the coefficient object.
- Formulate the two-stage estimation procedure (anchor via Fréchet mean, coefficient via MLE).
- Describe the generalized LIPO algorithm for global optimization in metric spaces and explain how the Hölder constant $K$ controls the exploration–exploitation tradeoff.
- Interpret the fitted geodesic as a discriminating direction and relate it to log-odds.
- Visualize and interpret the decision boundary on $\mathbb{S}^2$ as a great circle orthogonal to the discriminating direction.
- Apply metric-space logistic regression to classify market regimes from portfolio covariance matrices (SPD-valued covariates).
- Summarize the convergence rates and classification risk bounds.
## Motivation: Reversing the Roles
So far in this course, the non-Euclidean structure has been placed on the **response** side — $Y$ lives in a metric space, $X$ is Euclidean. A complementary problem reverses the roles: the response is **binary**, while the covariate itself lives in a metric space. This setting is developed by @LinLin2025Binary, who propose a logistic-type regression model and classifier for metric-space-valued covariates.
::: {.callout-tip title="Why this matters"}
Many modern datasets produce covariates that are not vectors: brain connectivity matrices (SPD), phylogenetic trees (BHV space), or distribution-valued features (Wasserstein space). When the outcome is binary (disease/no disease, task A/task B), we need a classifier that respects the geometry of the covariate space.
:::
## From Euclidean to Metric-Space Logistic Regression
Recall the Euclidean model with a binary response $Y \in \{0, 1\}$ and a covariate $X \in \mathbb{R}^D$. Writing $P_X = \mathbb{P}(Y = 1 \mid X)$ and suppressing the intercept, the usual logistic regression model is
$$
\log\frac{P_X}{1 - P_X} = X^\top \beta^*.
$$
Geometrically, if one anchors the vectors at the origin, then
$$
X^\top \beta^* = \|X\|_2 \|\beta^*\|_2 \cos\angle_0(X, \beta^*).
$$
This representation depends only on a **distance from the anchor**, a **distance to the coefficient**, and an **angle** between the corresponding directions. That observation suggests a route to generalization in a uniquely geodesic metric space.
## The Metric-Space Logistic Model
Let $X$ now take values in a uniquely geodesic metric space $(\mathcal{M}, d)$, and let $Y \in \{0, 1\}$. Because a general metric space has no distinguished origin, one first chooses an **anchor point**. A natural and statistically meaningful choice is the Fréchet mean of the covariate distribution,
$$
\mu^* \in \arg\min_{\omega \in \mathcal{M}} \mathbb{E}\bigl[d^2(X, \omega)\bigr],
$$
which is unique in Hadamard spaces under a finite second-moment condition.
::: {.definition title="Euclidean comparison angle"}
For three points $p, q, r$ in a metric space $(\mathcal{M}, d)$, form a Euclidean comparison triangle $\triangle(\tilde{p}, \tilde{q}, \tilde{r})$ in $\mathbb{R}^2$ whose side lengths match: $\|\tilde{q} - \tilde{r}\| = d(q, r)$, $\|\tilde{p} - \tilde{q}\| = d(p, q)$, and $\|\tilde{p} - \tilde{r}\| = d(p, r)$. The **Euclidean comparison angle** $\overline{\angle}_p(q, r)$ is the angle at $\tilde{p}$ in this triangle. By the law of cosines in $\mathbb{R}^2$,
$$
\overline{\angle}_p(q, r) = \arccos\!\left(\frac{d(p,q)^2 + d(p,r)^2 - d(q,r)^2}{2\,d(p,q)\,d(p,r)}\right).
$$
The comparison angle is well defined whenever $d(p,q) > 0$ and $d(p,r) > 0$. It measures what the angle *would be* if the three points were embedded isometrically in a flat Euclidean plane.
:::
::: {#def-alexandrov-angle .definition title="Alexandrov angle and inner product"}
Let $(\mathcal{M}, d)$ be a uniquely geodesic space. For $p, q, r \in \mathcal{M}$, let $\gamma_q, \gamma_r : [0, \infty) \to \mathcal{M}$ be the unit-speed geodesics from $p$ to $q$ and from $p$ to $r$, respectively. The **Alexandrov angle** at $p$ between $q$ and $r$ is
$$
\angle_p(q, r) = \limsup_{s, t \downarrow 0}\; \overline{\angle}_p\bigl(\gamma_q(s), \gamma_r(t)\bigr)
$$
where $\overline{\angle}_p$ is the Euclidean comparison angle defined above. Intuitively, one follows the two geodesics a tiny distance $s$ and $t$ from $p$, computes the comparison angle of the resulting three points ($p$, $\gamma_q(s)$, $\gamma_r(t)$), and takes the limiting value as $s, t \to 0$. In a smooth Riemannian manifold, the Alexandrov angle coincides with the usual Riemannian angle between tangent vectors.
The **Alexandrov inner product** of $q$ and $r$ with anchor $p$ is
$$
\langle q, r \rangle_p = d(p, q)\, d(p, r) \cos\angle_p(q, r)
$$
In Euclidean space $(\mathbb{R}^D, \|\cdot\|_2)$ with $p = 0$, this reduces to the ordinary inner product $q^\top r$.
:::
<!-- Figure: Comparison angle and Alexandrov angle -->
```{ojs}
//| label: fig-alexandrov-angle
//| fig-cap: "Comparison angle (left) and Alexandrov angle (right). Left: a Euclidean comparison triangle with the same three side lengths as the metric-space triple. Right: two geodesics emanating from p along directions toward q and r in a curved metric space."
html`
<div style="display: flex; gap: 40px; flex-wrap: wrap; justify-content: center; align-items: flex-start; max-width: 780px; font-family: system-ui, sans-serif;">
<!-- Comparison Angle Panel -->
<div style="text-align: center;">
<svg width="280" height="260" style="border: 1px solid #dee2e6; border-radius: 6px; background: #fafafa;">
${(() => {
// Euclidean comparison triangle
var pts = { px: 70, py: 210, qx: 210, qy: 210, rx: 170, ry: 60 };
var margin = 30;
function line(x1,y1,x2,y2,c,w,dash) {
return '<line x1="' + x1 + '" y1="' + y1 +
'" x2="' + x2 + '" y2="' + y2 +
'" stroke="' + c + '" stroke-width="' + (w||1.8) + '"' +
(dash ? ' stroke-dasharray="' + dash + '"' : '') + '/>';
}
var parts = [];
// Side labels
var midPqX = (pts.px + pts.qx)/2;
var midPqY = (pts.py + pts.qy)/2 + 18;
parts.push('<text x="' + midPqX + '" y="' + midPqY + '" text-anchor="middle" font-size="11" fill="#495057">d(p,q)</text>');
var midPrX = (pts.px + pts.rx)/2 - 28;
var midPrY = (pts.py + pts.ry)/2;
parts.push('<text x="' + midPrX + '" y="' + midPrY + '" text-anchor="middle" font-size="11" fill="#495057">d(p,r)</text>');
var midQrX = (pts.qx + pts.rx)/2 + 20;
var midQrY = (pts.qy + pts.ry)/2;
parts.push('<text x="' + midQrX + '" y="' + midQrY + '" text-anchor="middle" font-size="11" fill="#495057">d(q,r)</text>');
// Triangle edges
parts.push(line(pts.px, pts.py, pts.qx, pts.qy, '#212529', 1.8));
parts.push(line(pts.px, pts.py, pts.rx, pts.ry, '#212529', 1.8));
parts.push(line(pts.qx, pts.qy, pts.rx, pts.ry, '#212529', 1.8));
// Angle arc at p
var arcR = 25;
var vx1 = pts.qx - pts.px, vy1 = pts.qy - pts.py;
var vx2 = pts.rx - pts.px, vy2 = pts.ry - pts.py;
var a1 = Math.atan2(vy1, vx1);
var a2 = Math.atan2(vy2, vx2);
parts.push('<path d="M ' + (pts.px+arcR*Math.cos(a1)).toFixed(1) + ' ' + (pts.py+arcR*Math.sin(a1)).toFixed(1) +
' A ' + arcR + ' ' + arcR + ' 0 0 1 ' + (pts.px+arcR*Math.cos(a2)).toFixed(1) + ' ' + (pts.py+arcR*Math.sin(a2)).toFixed(1) +
'" fill="none" stroke="#e03131" stroke-width="2.2"/>');
var angleMid = (a1 + a2) / 2;
var labelR = arcR + 18;
parts.push('<text x="' + (pts.px + labelR*Math.cos(angleMid)).toFixed(0) +
'" y="' + (pts.py + labelR*Math.sin(angleMid)).toFixed(0) +
'" text-anchor="middle" font-size="13" dy="4" fill="#e03131" font-weight="bold">' +
'<tspan text-decoration="overline">∠</tspan></text>');
// Vertex labels
parts.push('<circle cx="' + pts.px + '" cy="' + pts.py + '" r="4" fill="#1971c2"/>');
parts.push('<text x="' + (pts.px-12) + '" y="' + (pts.py+18) + '" font-size="13" fill="#1971c2" font-weight="bold">p̃</text>');
parts.push('<circle cx="' + pts.qx + '" cy="' + pts.qy + '" r="4" fill="#2b8a3e"/>');
parts.push('<text x="' + (pts.qx+8) + '" y="' + (pts.qy+18) + '" font-size="13" fill="#2b8a3e" font-weight="bold">q̃</text>');
parts.push('<circle cx="' + pts.rx + '" cy="' + pts.ry + '" r="4" fill="#f08c00"/>');
parts.push('<text x="' + (pts.rx+8) + '" y="' + (pts.ry-8) + '" font-size="13" fill="#f08c00" font-weight="bold">r̃</text>');
// Subtitle
parts.push('<text x="140" y="20" text-anchor="middle" font-size="13" fill="#495057" font-weight="600">Euclidean Comparison Triangle</text>');
parts.push('<text x="140" y="36" text-anchor="middle" font-size="11" fill="#868e96">' +
'<tspan text-decoration="overline">∠</tspan>' +
'<tspan baseline-shift="sub" font-size="8">p</tspan>' +
'<tspan baseline-shift="baseline" font-size="11">(q,r) = arccos((d²+d²−d²)/(2·d·d))</tspan></text>');
return parts.join('\n');
})()}
</svg>
</div>
<!-- Alexandrov Angle Panel -->
<div style="text-align: center;">
<svg width="280" height="260" style="border: 1px solid #dee2e6; border-radius: 6px; background: #fafafa;">
${(() => {
var parts = [];
var cx = 100, cy = 195, R = 160;
// Curved background (metric space)
parts.push('<ellipse cx="' + cx + '" cy="' + cy + '" rx="' + R + '" ry="' + (R*0.6) +
'" fill="none" stroke="#dee2e6" stroke-width="1" opacity="0.7"/>');
parts.push('<ellipse cx="' + cx + '" cy="' + cy + '" rx="' + (R*0.6) + '" ry="' + (R*0.35) +
'" fill="none" stroke="#dee2e6" stroke-width="0.7" opacity="0.5"/>');
// Point p (anchor)
var px = 125, py = 190;
// Geodesic to q — quadratic bezier curve
var qx = 210, qy = 65;
var ctrlQx = 200, ctrlQy = 125;
parts.push('<path d="M ' + px + ' ' + py + ' Q ' + ctrlQx + ' ' + ctrlQy + ' ' + qx + ' ' + qy +
'" fill="none" stroke="#2b8a3e" stroke-width="2.2"/>');
// Geodesic to r — quadratic bezier curve
var rx = 55, ry = 80;
var ctrlRx = 80, ctrlRy = 135;
parts.push('<path d="M ' + px + ' ' + py + ' Q ' + ctrlRx + ' ' + ctrlRy + ' ' + rx + ' ' + ry +
'" fill="none" stroke="#f08c00" stroke-width="2.2"/>');
// Angle arc at p
var tanQx = ctrlQx - px, tanQy = ctrlQy - py;
var tanRx = ctrlRx - px, tanRy = ctrlRy - py;
var aQ = Math.atan2(tanQy, tanQx);
var aR = Math.atan2(tanRy, tanRx);
var arcR = 22;
parts.push('<path d="M ' + (px+arcR*Math.cos(aQ)).toFixed(1) + ' ' + (py+arcR*Math.sin(aQ)).toFixed(1) +
' A ' + arcR + ' ' + arcR + ' 0 0 0 ' + (px+arcR*Math.cos(aR)).toFixed(1) + ' ' + (py+arcR*Math.sin(aR)).toFixed(1) +
'" fill="none" stroke="#e03131" stroke-width="2.2"/>');
var angleMid = (aQ + aR) / 2;
var labelR = arcR + 18;
parts.push('<text x="' + (px + labelR*Math.cos(angleMid)).toFixed(0) +
'" y="' + (py + labelR*Math.sin(angleMid)).toFixed(0) +
'" text-anchor="middle" font-size="14" dy="4" fill="#e03131" font-weight="bold">∠</text>');
// Points
parts.push('<circle cx="' + px + '" cy="' + py + '" r="4.5" fill="#1971c2" stroke="#fff" stroke-width="1.2"/>');
parts.push('<text x="' + (px-16) + '" y="' + (py+20) + '" font-size="13" fill="#1971c2" font-weight="bold">p</text>');
parts.push('<circle cx="' + qx + '" cy="' + qy + '" r="4" fill="#2b8a3e" stroke="#fff" stroke-width="1"/>');
parts.push('<text x="' + (qx+8) + '" y="' + (qy-2) + '" font-size="13" fill="#2b8a3e" font-weight="bold">q</text>');
parts.push('<circle cx="' + rx + '" cy="' + ry + '" r="4" fill="#f08c00" stroke="#fff" stroke-width="1"/>');
parts.push('<text x="' + (rx-16) + '" y="' + (ry-4) + '" font-size="13" fill="#f08c00" font-weight="bold">r</text>');
// Dashed tangent directions
var sLen = 35;
var qsx = px + sLen * Math.cos(aQ);
var qsy = py + sLen * Math.sin(aQ);
var rsx = px + sLen * Math.cos(aR);
var rsy = py + sLen * Math.sin(aR);
parts.push('<line x1="' + px + '" y1="' + py +
'" x2="' + qsx.toFixed(0) + '" y2="' + qsy.toFixed(0) +
'" stroke="#adb5bd" stroke-width="0.8" stroke-dasharray="4,3"/>');
parts.push('<line x1="' + px + '" y1="' + py +
'" x2="' + rsx.toFixed(0) + '" y2="' + rsy.toFixed(0) +
'" stroke="#adb5bd" stroke-width="0.8" stroke-dasharray="4,3"/>');
// Subtitle
parts.push('<text x="140" y="20" text-anchor="middle" font-size="13" fill="#495057" font-weight="600">Alexandrov Angle in Metric Space</text>');
parts.push('<text x="140" y="36" text-anchor="middle" font-size="11" fill="#868e96">' +
'<tspan>∠</tspan><tspan baseline-shift="sub" font-size="8">p</tspan>' +
'<tspan baseline-shift="baseline" font-size="11">(q,r) = limsup </tspan>' +
'<tspan text-decoration="overline">∠</tspan>' +
'<tspan baseline-shift="sub" font-size="8">p</tspan>' +
'<tspan baseline-shift="baseline" font-size="11">(γ</tspan>' +
'<tspan baseline-shift="sub" font-size="8">q</tspan>' +
'<tspan baseline-shift="baseline" font-size="11">(s), γ</tspan>' +
'<tspan baseline-shift="sub" font-size="8">r</tspan>' +
'<tspan baseline-shift="baseline" font-size="11">(t))</tspan></text>');
return parts.join('\n');
})()}
</svg>
</div>
</div>
`
```
With these definitions in hand, the **metric-space logistic model** is
$$
\log\frac{P_X}{1 - P_X} = h(\beta^*; X, \mu^*) = d(\mu^*, X)\, d(\mu^*, \beta^*) \cos\angle_{\mu^*}(X, \beta^*).
$$
The right-hand side is the **Alexandrov inner product** $\langle X, \beta^* \rangle_{\mu^*}$ defined above. In Euclidean space it reduces to
$$
h(\beta^*; X, \mu^*) = (X - \mu^*)^\top(\beta^* - \mu^*),
$$
so the model is a genuine extension of logistic regression.
::: {.callout-note title="Geometric interpretation"}
The geodesic from $\mu^*$ to $\beta^*$ plays the role of a **discriminating direction**: points $X$ whose geodesic from $\mu^*$ has a large positive projection onto this direction increase the log odds of $Y = 1$, while points pointing in the opposite direction decrease it. The scalar $d(\mu^*, \beta^*)$ controls the **magnitude** of this effect.
:::
## Estimation
Given i.i.d. observations $(X_1, Y_1), \ldots, (X_n, Y_n)$, estimation proceeds in two stages:
**Stage 1 — Anchor estimation.** Estimate the anchor by the sample Fréchet mean:
$$
\hat{\mu} \in \arg\min_{\omega \in \mathcal{M}} \frac{1}{n}\sum_{i=1}^n d^2(X_i, \omega).
$$
**Stage 2 — Coefficient estimation.** Define the population log-likelihood
$$
L(\beta) = \mathbb{E}\Bigl[Y\, h(\beta; X, \mu^*) - \log\bigl(1 + e^{h(\beta; X, \mu^*)}\bigr)\Bigr]
$$
and its empirical counterpart
$$
L_n(\beta) = \frac{1}{n}\sum_{i=1}^n \Bigl[Y_i\, h(\beta; X_i, \hat{\mu}) - \log\bigl(1 + e^{h(\beta; X_i, \hat{\mu})}\bigr)\Bigr].
$$
The estimator of the coefficient object is the metric-space maximum likelihood estimator:
$$
\hat{\beta} \in \arg\max_{\beta \in \mathcal{M}} L_n(\beta).
$$
Under correct specification, the true parameter maximizes the population log-likelihood.
::: {#prp-logistic-pop .proposition title="Population likelihood characterization"}
If $(X, Y)$ follows the metric-space logistic model with coefficient $\beta^*$ and anchor $\mu^*$, then
$$
\beta^* \in \arg\max_{\beta \in \mathcal{M}} L(\beta).
$$
If the maximizer of $L$ is unique, then $\beta^*$ is exactly that unique maximizer.
:::
## Computation: Generalized LIPO Algorithm
The new computational challenge is that $L_n$ is usually **neither convex nor differentiable** in any linear sense on a general metric space. To address this, @LinLin2025Binary extend the LIPO (Lipschitz optimization) idea to metric spaces.
Assume $L_n$ satisfies a Hölder continuity bound:
$$
|L_n(\beta_1) - L_n(\beta_2)| \le K\, d(\beta_1, \beta_2)^{\alpha_U}, \qquad \alpha_U \in (0, 1],
$$
for some known $K > 0$.
::: {#alg-metric-lipo .algorithm title="Generalized LIPO algorithm on a metric space"}
Fix an iteration budget $T \in \mathbb{N}$, a Hölder constant $K > 0$, and an exponent $\alpha_U \in (0, 1]$.
1. Sample $\beta_{(0)}$ uniformly from $\mathcal{M}$.
2. Set $t \leftarrow 0$.
3. While $t < T$, repeat:
- Sample a proposal $\beta$ uniformly from $\mathcal{M}$.
- If $\min_{j=0,\ldots,t}\bigl(L_n(\beta_{(j)}) + K\, d(\beta, \beta_{(j)})^{\alpha_U}\bigr) \ge \max_{j=0,\ldots,t} L_n(\beta_{(j)})$,
then accept: $\beta_{(t+1)} = \beta$, $t \leftarrow t + 1$.
4. Return $\beta_{(\tau)}$, where $\tau \in \arg\max_{0 \le t \le T} L_n(\beta_{(t)})$.
:::
The acceptance rule keeps only those candidates whose optimistic Hölder upper envelope can still match or beat the best objective value seen so far. In Hadamard spaces one may take $C_U = 5$ and $\alpha_U = 1$, with
$$
K = \frac{2 C_U}{n}\sum_{i=1}^n d(X_i, \hat{\mu}).
$$
::: {.callout-tip}
When additional Riemannian structure is available, one may also exploit manifold optimization methods, but the appeal of the LIPO-type strategy is that it **only uses the metric**.
:::
### Interactive Demo: LIPO Algorithm in Action
The following demo visualizes the generalized LIPO algorithm on the circle $S^1$, where the parameter $\beta$ is an angle. The log-likelihood $L_n(\beta)$ is shown as a black curve. As the algorithm runs, it evaluates proposals and builds a **Hölder upper envelope** (shaded region). A proposal is accepted only if its upper envelope can still beat the current best value — this prunes the search space adaptively.
```{ojs}
//| label: fig-lipo-demo
//| fig-cap: "Interactive: Generalized LIPO algorithm searching for the MLE on S¹"
lipo_n_control = Inputs.range([30, 300], {step: 10, value: 150, label: "Sample size n"})
lipo_K_factor_control = Inputs.range([0.2, 3.0], {step: 0.1, value: 1.0, label: "Hölder constant factor (× true K)"})
lipo_budget_control = Inputs.range([5, 80], {step: 1, value: 25, label: "LIPO budget T"})
lipo_seed_control = Inputs.range([1, 50], {step: 1, value: 7, label: "Random seed"})
lipo_n = Generators.input(lipo_n_control)
lipo_K_factor = Generators.input(lipo_K_factor_control)
lipo_budget = Generators.input(lipo_budget_control)
lipo_seed = Generators.input(lipo_seed_control)
lipo_controls_view = html`
<style>
.lipo-controls-shell {
width:100%;
margin:0 0 14px;
font-family:system-ui,sans-serif;
container-type:inline-size;
}
.lipo-controls-grid {
display:grid;
grid-template-columns:repeat(2,minmax(0,1fr));
gap:8px 24px;
width:100%;
}
.lipo-controls-grid > * {
min-width:0;
margin:0;
}
.lipo-controls-grid input[type="number"] {
width:7.5rem !important;
}
@container (max-width:560px) {
.lipo-controls-grid {
grid-template-columns:1fr;
}
}
</style>
<div class="lipo-controls-shell">
<div class="lipo-controls-grid">
<div>${lipo_n_control}</div>
<div>${lipo_K_factor_control}</div>
<div>${lipo_budget_control}</div>
<div>${lipo_seed_control}</div>
</div>
</div>
`
function runLIPODemo(n, Kfactor, T, seed) {
// ---- PRNG ----
function mulberry32(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;
};
}
var rng = mulberry32(seed);
// ---- Generate data on S¹ ----
// True parameters
var muStar = 0;
var betaTrue = 1.2; // ~69 degrees
// Circular utilities
function circDist(a, b) {
var d = Math.abs(a - b) % (2 * Math.PI);
return Math.min(d, 2 * Math.PI - d);
}
function alexIP(beta, x, anchor) {
var dx = circDist(anchor, x);
var db = circDist(anchor, beta);
// Alexandrov angle on S¹
var dirX = x - anchor;
if (dirX > Math.PI) dirX -= 2 * Math.PI;
if (dirX < -Math.PI) dirX += 2 * Math.PI;
var dirBeta = beta - anchor;
if (dirBeta > Math.PI) dirBeta -= 2 * Math.PI;
if (dirBeta < -Math.PI) dirBeta += 2 * Math.PI;
var cosAngle = Math.cos(dirX) * Math.cos(dirBeta) + Math.sin(dirX) * Math.sin(dirBeta);
return dx * db * cosAngle;
}
// Generate covariates concentrated around anchor
var X = [];
var kappa = 3.0;
for (var i = 0; i < n; i++) {
// von Mises sampling (simple rejection)
var theta;
var done = false;
while (!done) {
theta = rng() * 2 * Math.PI;
var vmDens = Math.exp(kappa * Math.cos(theta - muStar));
if (rng() < vmDens / Math.exp(kappa)) done = true;
}
X.push(theta);
}
// Generate binary responses
var Y = [];
for (var i = 0; i < n; i++) {
var logit = alexIP(betaTrue, X[i], muStar);
var prob = 1 / (1 + Math.exp(-logit));
Y.push(rng() < prob ? 1 : 0);
}
// Sample Fréchet mean (anchor)
var muHat = 0; // known for this demo
for (var iter = 0; iter < 20; iter++) {
var sumDir = 0;
for (var i = 0; i < n; i++) {
var d = X[i] - muHat;
if (d > Math.PI) d -= 2 * Math.PI;
if (d < -Math.PI) d += 2 * Math.PI;
sumDir += d;
}
muHat = (muHat + sumDir / n + 2 * Math.PI) % (2 * Math.PI);
}
// ---- Log-likelihood on a fine grid ----
var nFine = 360;
var llGrid = [];
var betaGrid = [];
for (var g = 0; g < nFine; g++) {
var beta = g * 2 * Math.PI / nFine;
var ll = 0;
for (var i = 0; i < n; i++) {
var h = alexIP(beta, X[i], muHat);
ll += Y[i] * h - Math.log(1 + Math.exp(h));
}
ll /= n;
betaGrid.push(beta);
llGrid.push(ll);
}
var llMax = Math.max.apply(null, llGrid);
var llMin = Math.min.apply(null, llGrid);
var bestIdx = llGrid.indexOf(llMax);
var betaMLE = betaGrid[bestIdx];
// ---- True Hölder constant K ----
// K = (2*C_U / n) * sum_i d(X_i, muHat) with C_U=1 on S¹
var Ktrue = 0;
for (var i = 0; i < n; i++) {
Ktrue += circDist(X[i], muHat);
}
Ktrue = 2 * Ktrue / n;
var Kused = Kfactor * Ktrue;
// ---- Run LIPO ----
var evaluated = []; // {beta, ll}
var accepted = []; // indices into evaluated
var rejected = []; // {beta, upperBound}
// Start with random point
var b0 = rng() * 2 * Math.PI;
var ll0 = 0;
for (var i = 0; i < n; i++) {
var h = alexIP(b0, X[i], muHat);
ll0 += Y[i] * h - Math.log(1 + Math.exp(h));
}
ll0 /= n;
evaluated.push({beta: b0, ll: ll0});
accepted.push(0);
var t = 0;
var totalProposals = 0;
while (t < T && totalProposals < 2000) {
totalProposals++;
var prop = rng() * 2 * Math.PI;
// Compute upper envelope: min_j (ll_j + K * d(prop, beta_j)^alpha)
var bestVal = -Infinity;
for (var j = 0; j < evaluated.length; j++) {
if (evaluated[j].ll > bestVal) bestVal = evaluated[j].ll;
}
var upperBound = Infinity;
for (var j = 0; j < evaluated.length; j++) {
var ub = evaluated[j].ll + Kused * circDist(prop, evaluated[j].beta);
if (ub < upperBound) upperBound = ub;
}
if (upperBound >= bestVal) {
// Accept
var llProp = 0;
for (var i = 0; i < n; i++) {
var h = alexIP(prop, X[i], muHat);
llProp += Y[i] * h - Math.log(1 + Math.exp(h));
}
llProp /= n;
evaluated.push({beta: prop, ll: llProp});
accepted.push(evaluated.length - 1);
t++;
} else {
rejected.push({beta: prop, upperBound: upperBound});
}
}
// LIPO best
var lipoBest = -Infinity, lipoBestBeta = 0;
for (var a = 0; a < accepted.length; a++) {
var ev = evaluated[accepted[a]];
if (ev.ll > lipoBest) { lipoBest = ev.ll; lipoBestBeta = ev.beta; }
}
// ---- Build upper envelope curve for display ----
var envelopeCurve = [];
for (var g = 0; g < nFine; g++) {
var beta = g * 2 * Math.PI / nFine;
var ub = Infinity;
for (var j = 0; j < evaluated.length; j++) {
var u = evaluated[j].ll + Kused * circDist(beta, evaluated[j].beta);
if (u < ub) ub = u;
}
envelopeCurve.push(ub);
}
return {
X: X, Y: Y, muHat: muHat, betaTrue: betaTrue,
betaGrid: betaGrid, llGrid: llGrid, llMin: llMin, llMax: llMax,
betaMLE: betaMLE,
evaluated: evaluated, accepted: accepted, rejected: rejected,
envelopeCurve: envelopeCurve,
lipoBest: lipoBest, lipoBestBeta: lipoBestBeta,
Ktrue: Ktrue, Kused: Kused, T: T, totalProposals: totalProposals
};
}
lipoResult = runLIPODemo(lipo_n, lipo_K_factor, lipo_budget, lipo_seed);
html`
<div style="font-family: system-ui, sans-serif; max-width: 880px;">
<div style="display: flex; gap: 20px; flex-wrap: wrap;">
<!-- Main LIPO visualization -->
<div>
<h4>LIPO Search on S¹</h4>
<svg width="480" height="310" style="border: 1px solid #dee2e6; border-radius: 4px;">
${(() => {
var margin = {top: 20, right: 20, bottom: 35, left: 50};
var pw = 480 - margin.left - margin.right;
var ph = 310 - margin.top - margin.bottom;
function xS(b) { return margin.left + (b / (2 * Math.PI)) * pw; }
var res = lipoResult;
var vMin = res.llMin, vMax = res.llMax;
var range = vMax - vMin || 0.01;
function yS(v) { return margin.top + ph - ((v - vMin) / range) * ph; }
var parts = [];
// Background
parts.push('<rect x="' + margin.left + '" y="' + margin.top + '" width="' + pw +
'" height="' + ph + '" fill="#fafafa" rx="2"/>');
// Upper envelope (shaded)
var envPath = '';
for (var g = 0; g < res.envelopeCurve.length; g++) {
var bx = xS(res.betaGrid[g]);
var by = yS(res.envelopeCurve[g]);
envPath += (g === 0 ? 'M' : 'L') + ' ' + bx.toFixed(1) + ' ' + by.toFixed(1);
}
// Fill from envelope down to bottom
parts.push('<path d="' + envPath + ' L ' + xS(2*Math.PI).toFixed(1) + ' ' + yS(vMin).toFixed(1) +
' L ' + xS(0).toFixed(1) + ' ' + yS(vMin).toFixed(1) + ' Z" fill="#ffd43b" opacity="0.18"/>');
parts.push('<path d="' + envPath + '" fill="none" stroke="#f08c00" stroke-width="1.2" stroke-dasharray="5,3" opacity="0.8"/>');
// True log-likelihood curve
var llPath = '';
for (var g = 0; g < res.betaGrid.length; g++) {
llPath += (g === 0 ? 'M' : 'L') + ' ' + xS(res.betaGrid[g]).toFixed(1) + ' ' + yS(res.llGrid[g]).toFixed(1);
}
parts.push('<path d="' + llPath + '" fill="none" stroke="#212529" stroke-width="2.2"/>');
// Rejected proposals (gray X marks at their upper bound estimate)
for (var r = 0; r < res.rejected.length; r++) {
var rx = xS(res.rejected[r].beta);
var ry = yS(res.rejected[r].upperBound);
var s = 3;
parts.push('<line x1="' + (rx-s) + '" y1="' + (ry-s) + '" x2="' + (rx+s) + '" y2="' + (ry+s) + '" stroke="#adb5bd" stroke-width="1"/>');
parts.push('<line x1="' + (rx-s) + '" y1="' + (ry+s) + '" x2="' + (rx+s) + '" y2="' + (ry-s) + '" stroke="#adb5bd" stroke-width="1"/>');
}
// Evaluated/accepted points
for (var a = 0; a < res.accepted.length; a++) {
var ev = res.evaluated[res.accepted[a]];
var ex = xS(ev.beta);
var ey = yS(ev.ll);
var isBest = ev.ll >= res.lipoBest - 1e-10;
if (a === 0) {
// First point: initial random
parts.push('<circle cx="' + ex + '" cy="' + ey + '" r="5" fill="#1971c2" stroke="#0b3d91" stroke-width="1.5"/>');
} else if (isBest) {
parts.push('<circle cx="' + ex + '" cy="' + ey + '" r="5.5" fill="#e03131" stroke="#a61e24" stroke-width="1.5"/>');
} else {
parts.push('<circle cx="' + ex + '" cy="' + ey + '" r="4" fill="#1971c2" opacity="0.7"/>');
}
}
// True β* marker
var tx = xS(res.betaTrue);
parts.push('<line x1="' + tx + '" y1="' + margin.top + '" x2="' + tx + '" y2="' + (margin.top+ph) +
'" stroke="#2b8a3e" stroke-width="1.2" stroke-dasharray="4,4"/>');
parts.push('<text x="' + tx + '" y="' + (margin.top-5) + '" text-anchor="middle" font-size="9" fill="#2b8a3e">β*</text>');
// Current best line
var bestY = yS(res.lipoBest);
parts.push('<line x1="' + margin.left + '" y1="' + bestY + '" x2="' + (margin.left+pw) +
'" y2="' + bestY + '" stroke="#e03131" stroke-width="0.8" stroke-dasharray="3,3"/>');
// Axes
parts.push('<line x1="' + margin.left + '" y1="' + (margin.top+ph) + '" x2="' + (margin.left+pw) +
'" y2="' + (margin.top+ph) + '" stroke="#adb5bd"/>');
parts.push('<line x1="' + margin.left + '" y1="' + margin.top + '" x2="' + margin.left +
'" y2="' + (margin.top+ph) + '" stroke="#adb5bd"/>');
// X-axis labels
var labels = [{v:0,l:'0'},{v:Math.PI/2,l:'π/2'},{v:Math.PI,l:'π'},{v:3*Math.PI/2,l:'3π/2'},{v:2*Math.PI,l:'2π'}];
for (var li = 0; li < labels.length; li++) {
parts.push('<text x="' + xS(labels[li].v) + '" y="' + (margin.top+ph+14) +
'" text-anchor="middle" font-size="9">' + labels[li].l + '</text>');
}
parts.push('<text x="' + (margin.left+pw/2) + '" y="' + (margin.top+ph+28) +
'" text-anchor="middle" font-size="10">β (radians)</text>');
parts.push('<text x="' + (margin.left-42) + '" y="' + (margin.top+ph/2) +
'" text-anchor="middle" font-size="10" transform="rotate(-90,' + (margin.left-42) + ',' + (margin.top+ph/2) + ')">Log-likelihood</text>');
// Legend
var lx = margin.left + pw - 195, ly = margin.top + 5;
parts.push('<circle cx="' + (lx+4) + '" cy="' + (ly) + '" r="3.5" fill="#1971c2" opacity="0.7"/>');
parts.push('<text x="' + (lx+10) + '" y="' + (ly+4) + '" font-size="8" fill="#616161">Evaluated</text>');
parts.push('<circle cx="' + (lx+68) + '" cy="' + (ly) + '" r="4" fill="#e03131"/>');
parts.push('<text x="' + (lx+74) + '" y="' + (ly+4) + '" font-size="8" fill="#616161">Best</text>');
parts.push('<line x1="' + (lx+106) + '" y1="' + (ly) + '" x2="' + (lx+118) + '" y2="' + (ly) +
'" stroke="#f08c00" stroke-width="1" stroke-dasharray="4,2"/>');
parts.push('<text x="' + (lx+121) + '" y="' + (ly+4) + '" font-size="8" fill="#616161">Envelope</text>');
parts.push('<rect x="' + margin.left + '" y="' + margin.top + '" width="' + pw +
'" height="' + ph + '" fill="none" stroke="#dee2e6"/>');
return parts.join('\n');
})()}
</svg>
</div>
<!-- Right panel: stats -->
<div style="flex: 1; min-width: 270px;">
<div style="padding: 12px; background: #f8f9fa; border-radius: 6px;">
<table style="width: 100%; border-collapse: collapse; font-size: 0.9em;">
<tr><td style="padding: 3px 8px;">True β*</td>
<td style="padding: 3px 8px; text-align: right;">${(lipoResult.betaTrue * 180 / Math.PI).toFixed(1)}°</td></tr>
<tr><td style="padding: 3px 8px;">MLE (grid search)</td>
<td style="padding: 3px 8px; text-align: right;">${(lipoResult.betaMLE * 180 / Math.PI).toFixed(1)}°</td></tr>
<tr><td style="padding: 3px 8px; font-weight: bold; color: #e03131;">LIPO best β̂</td>
<td style="padding: 3px 8px; text-align: right; font-weight: bold; color: #e03131;">${(lipoResult.lipoBestBeta * 180 / Math.PI).toFixed(1)}°</td></tr>
<tr><td colspan="2"><hr style="margin: 4px 0;"></td></tr>
<tr><td style="padding: 3px 8px;">Accepted proposals</td>
<td style="padding: 3px 8px; text-align: right;">${lipoResult.accepted.length}</td></tr>
<tr><td style="padding: 3px 8px;">Rejected proposals</td>
<td style="padding: 3px 8px; text-align: right;">${lipoResult.rejected.length}</td></tr>
<tr><td style="padding: 3px 8px;">Total proposals</td>
<td style="padding: 3px 8px; text-align: right;">${lipoResult.totalProposals}</td></tr>
<tr><td style="padding: 3px 8px;">Acceptance rate</td>
<td style="padding: 3px 8px; text-align: right;">${(lipoResult.accepted.length / lipoResult.totalProposals * 100).toFixed(1)}%</td></tr>
<tr><td colspan="2"><hr style="margin: 4px 0;"></td></tr>
<tr><td style="padding: 3px 8px;">True K</td>
<td style="padding: 3px 8px; text-align: right;">${lipoResult.Ktrue.toFixed(3)}</td></tr>
<tr><td style="padding: 3px 8px;">Used K (factor × true)</td>
<td style="padding: 3px 8px; text-align: right;">${lipoResult.Kused.toFixed(3)}</td></tr>
<tr><td style="padding: 3px 8px;">max L_n (LIPO)</td>
<td style="padding: 3px 8px; text-align: right; font-family: monospace;">${lipoResult.lipoBest.toFixed(4)}</td></tr>
<tr><td style="padding: 3px 8px;">max L_n (grid)</td>
<td style="padding: 3px 8px; text-align: right; font-family: monospace;">${lipoResult.llMax.toFixed(4)}</td></tr>
</table>
</div>
</div>
</div>
<p style="margin-top: 12px; font-size: 0.9em; color: #495057;">
<b>How to read the plot:</b> The black curve is the true log-likelihood $L_n(\beta)$.
Orange dashed curve = Hölder upper envelope built from evaluated points. Blue dots = accepted
LIPO evaluations (the first is the initial random point). Red dot = current best. Gray × marks = rejected
proposals (their upper envelope bound fell below the current best, so they can't improve on it).
The green dashed line marks the true $\beta^*$.
</p>
</div>
`
```
::: {.callout-tip title="Try these experiments"}
- **Set K factor = 0.5:** The envelope is too tight — the algorithm rejects proposals that could actually beat the current best, potentially missing the global maximum (premature convergence).
- **Set K factor = 2.5:** The envelope is too loose — almost every proposal passes the test, and the algorithm degenerates to uniform random search (acceptance rate near 100%).
- **Increase T:** More evaluations refine the envelope and drive the best found value closer to the true MLE.
- **Compare LIPO best vs. grid MLE:** With enough budget and a well-chosen K, LIPO approaches the grid-search optimum without evaluating every grid point.
:::
## Application: Task-Based fMRI Classification
Section 5 of @LinLin2025Binary illustrates the method using task-based fMRI data from the Human Connectome Project.
### Data
- After preprocessing BOLD signals from 8 regions of interest, each subject is represented by a **functional-connectivity matrix** $X_i \in \mathcal{S}_{++}^8$.
- $Y_i = 0$: motor task ($n_0 = 201$ subjects)
- $Y_i = 1$: language-processing task ($n_1 = 193$ subjects)
- Total $n = 394$.
The model is fitted using the **Log-Cholesky metric** on $\mathcal{S}_{++}^8$ [@LinSPDCholesky2019]; the log-Euclidean and affine-invariant metrics lead to similar conclusions.
### The Fitted Geodesic
::: {.text-center}
{width=80%}
{width=60%}
:::
The fitted geodesic is **normal to the empirical decision boundary**: the positive direction points toward the language-processing task, while the negative direction points toward the motor task. Replacing the connectivity matrix $\gamma_{\hat{\mu}}^{\hat{\beta}}(t)$ by $\gamma_{\hat{\mu}}^{\hat{\beta}}(t + \delta)$ with $\delta = 0.00375$ increases the fitted log odds by $1.486$.
### Neuroanatomical Interpretability
The fitted effect is anatomically interpretable:
- **Region 1** (superior temporal gyrus) becomes markedly more active in the positive direction — consistent with its role in auditory and language processing.
- The correlation between **Region 1** and **Region 4** (inferior frontal gyrus, associated with speech/language) becomes visibly stronger in the positive direction.
- Correlations between Region 1 and Regions 5–8 (primary motor and visual cortices) become **weaker** in the positive direction — suggesting that language processing is accompanied by relative suppression of connectivity between superior temporal gyrus and motor/visual areas.
### Classification Performance
{width=60%}
The in-sample ROC curve above shows that the proposed model is adequate to classify the two tasks based on functional connectivity. For a comprehensive evaluation, the proposed classifier is compared with a diverse set of methods using **10-fold cross-validated** point measures (accuracy, sensitivity, specificity) and summary measures (AUC, Brier score, logarithmic score).
**Competing methods.** The comparison includes:
- **SPDNet** [@huang2017riemannian]: a neural network designed specifically for SPD-matrix inputs.
- **$k$-NN**: $k$-nearest neighbors classifier on the SPD manifold [@chaudhuri2014rates].
- **kSVM**: kernel support vector machine with a Gaussian kernel adapted to SPD matrices [@jayasumana2013kernel].
- **KDC**: kernel density classifier [@hastie2009elements].
- **RF**: random forest using the vectorized upper triangle of the matrix logarithm.
- **Logistic**: Euclidean logistic regression on the vectorized upper triangle of the matrix logarithm.
- **Logistic-BOLD**: Euclidean logistic regression using only the diagonal entries (BOLD signals) of each $X_i$, ignoring inter-regional connectivity.
- **GLMAF-LR / GLMAF-RF**: task-related GLM activation features [@barch2013function] fed into logistic regression and random forest.
- **ALFF-LR / ALFF-RF**: amplitude of low-frequency fluctuation features [@yu2007altered] fed into logistic regression and random forest.
**Results.** The full comparison is summarized in @tbl-fmri-classification.
| Method | Accuracy | Sensitivity | Specificity | AUC | Brier | LogS |
|---------------|:--------:|:-----------:|:-----------:|:------:|:-------:|:-------:|
| SPDNet | <u>0.8667</u> | 0.8613 | **0.8690** | **0.9205** | −0.1278 | −0.4600 |
| $k$-NN | 0.8333 | 0.8401 | 0.8297 | 0.8349 | −0.0062 | −0.4382 |
| kSVM | 0.8462 | <u>0.8667</u> | 0.8266 | 0.8466 | <u>−0.0039</u> | **−0.3458** |
| KDC | 0.7744 | 0.8313 | 0.7309 | 0.7811 | −0.0072 | −0.8431 |
| RF | 0.8179 | 0.8536 | 0.7930 | 0.8913 | −0.0067 | −0.4487 |
| Logistic | 0.8205 | 0.8617 | 0.7878 | 0.8647 | −0.0065 | −0.8565 |
| Logistic-BOLD | 0.7641 | 0.8490 | 0.6984 | 0.8120 | −0.0109 | −0.8002 |
| GLMAF-LR | 0.7241 | 0.7467 | 0.7147 | 0.7902 | −0.0063 | −0.5843 |
| GLMAF-RF | 0.7215 | 0.7400 | 0.7088 | 0.8019 | −0.0056 | −0.5528 |
| ALFF-LR | 0.7266 | 0.7109 | 0.7459 | 0.7881 | −0.0047 | −0.6455 |
| ALFF-RF | 0.7557 | 0.7284 | 0.7784 | 0.8478 | **−0.0024** | −0.4894 |
| **Proposed** | **0.8718** | **0.9158** | <u>0.8376</u> | <u>0.9060</u> | −0.0060 | <u>−0.4200</u> |
: Comparison of classification performance for task-based fMRI data. **Bold** = best, <u>underlined</u> = second best. Source: @LinLin2025Binary. {#tbl-fmri-classification}
**Key observations:**
1. **The proposed classifier achieves the best accuracy (0.8718) and sensitivity (0.9158)**, and ranks second in specificity (0.8376), AUC (0.9060), and logarithmic score (−0.4200). It outperforms all non-neural-network competitors across accuracy, sensitivity, and specificity.
2. **SPDNet achieves the best AUC (0.9205) and specificity (0.8690)**, but the proposed method is competitive — trailing by only 0.0145 in AUC and 0.0314 in specificity — while offering substantially greater geometric interpretability.
3. **Nonparametric classifiers suffer from the curse of dimensionality.** The manifold $\mathcal{S}_{++}^8$ has intrinsic dimension 36, which degrades the performance of $k$-NN, KDC, and even the kernel SVM relative to the parametric proposed model.
4. **Connectivity structure matters decisively.** The Logistic-BOLD baseline (accuracy 0.7641, AUC 0.8120) ignores off-diagonal entries and performs substantially worse than the proposed model (accuracy 0.8718, AUC 0.9060). The gap of over 10 percentage points in accuracy demonstrates that **inter-regional functional connectivity** — encoded in the off-diagonal entries of each $X_i$ — carries the discriminating information between motor and language tasks.
5. **Alternative feature extraction underperforms.** Methods based on GLM activation features (GLMAF-LR, GLMAF-RF) and ALFF features (ALFF-LR, ALFF-RF) all achieve substantially lower accuracy (0.72–0.76) than the proposed functional-connectivity-based classifier, confirming that the SPD-matrix representation is more informative for this task.
::: {.callout-important title="Key finding"}
Even relative to SPDNet — a purpose-built neural network for SPD inputs — the metric-space logistic classifier is **competitive in predictive performance** while remaining **far more interpretable**: the fitted geodesic $\gamma_{\hat{\mu}}^{\hat{\beta}}$ directly reveals *which* connectivity patterns drive the classification decision and *how*. This interpretability is a central advantage of the model-based approach.
:::
## Application: Classifying Market Regimes from Portfolio Covariance
The portfolio stress-testing framework developed in Lectures 6–9 models how the **covariance matrix** of asset returns varies with a market stress indicator. But a complementary question arises naturally: given only a covariance matrix (estimated from a window of returns), can we **classify** whether the market was in a crisis regime or a normal regime? This is a binary classification problem where the covariate $X$ lives in $\mathcal{S}_{++}^p$ and the response $Y$ is the regime label.
### Setup
Consider a portfolio of $p = 3$ assets (e.g., equity index, bond index, commodity index). During crisis periods, correlations spike and volatilities increase; during calm periods, correlations are moderate and volatilities are lower. We observe a $3 \times 3$ SPD covariance matrix $\Sigma_t$ for each time window $t$, and we want to classify the regime.
Under the **log-Euclidean metric**, the SPD manifold is isometric to Euclidean space (the space of symmetric matrices), so the Alexandrov inner product simplifies. For $X = \log\Sigma$ (the matrix logarithm) and anchor $\mu^*$ (the Fréchet mean of log-covariances):
$$
h(\beta^*; X, \mu^*) = \langle X - \mu^*, \beta^* - \mu^* \rangle_F = \operatorname{tr}\bigl((X - \mu^*)^\top(\beta^* - \mu^*)\bigr),
$$
where $\langle\cdot,\cdot\rangle_F$ is the Frobenius inner product. This is ordinary logistic regression in the log-domain — but the geometry of the SPD manifold shapes the interpretation: the fitted coefficient $\hat{\beta}$ (exponentiated) is an SPD matrix whose eigendirections reveal which portfolio weights drive the classification.
### Interactive Demo: Regime Classification on Simulated SPD Data
The following demo simulates $3 \times 3$ SPD covariance matrices from two regimes (normal and crisis), fits the metric-space logistic classifier, and visualizes the discriminating direction in terms of its effect on portfolio variance.
```{ojs}
//| label: fig-portfolio-classify
//| fig-cap: "Interactive: Metric-space logistic classification of market regimes from 3×3 SPD covariance matrices"
pc_n_control = Inputs.range([40, 300], {step: 10, value: 120, label: "Sample size per regime"})
pc_sep_control = Inputs.range([0.1, 2.0], {step: 0.1, value: 0.7, label: "Regime separation"})
pc_noise_control = Inputs.range([0.05, 0.4], {step: 0.02, value: 0.15, label: "Noise level σ"})
pc_seed_control = Inputs.range([1, 50], {step: 1, value: 3, label: "Random seed"})
pc_n = Generators.input(pc_n_control)
pc_sep = Generators.input(pc_sep_control)
pc_noise = Generators.input(pc_noise_control)
pc_seed = Generators.input(pc_seed_control)
pc_controls_view = html`
<style>
.pc-slider-grid {
display: grid;
grid-template-columns: repeat(2, minmax(0, 1fr));
gap: 6px 20px;
max-width: 620px;
margin: 0 auto 12px;
font-family: system-ui, sans-serif;
font-size: 0.85em;
}
.pc-slider-grid > * {
min-width: 0;
margin: 0;
}
.pc-slider-grid input[type="number"] {
width: 7.5rem !important;
}
@media (max-width: 500px) {
.pc-slider-grid {
grid-template-columns: 1fr;
}
}
</style>
<div class="pc-slider-grid">
<div>${pc_n_control}</div>
<div>${pc_sep_control}</div>
<div>${pc_noise_control}</div>
<div>${pc_seed_control}</div>
</div>
`
function runPortfolioClassify(nPerRegime, sep, noise, 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);
function randn() {
var u1 = rng(), u2 = rng();
return Math.sqrt(-2 * Math.log(Math.max(u1, 1e-15))) * Math.cos(2 * Math.PI * u2);
}
// 3x3 symmetric matrix stored as [a,b,c, d,e, f] for [[a,b,c],[b,d,e],[c,e,f]]
function symFrobSq(m) {
return m[0]*m[0]+m[1]*m[1]+m[2]*m[2]+m[3]*m[3]+m[4]*m[4]+m[5]*m[5] + 2*(m[1]*m[1]+m[2]*m[2]+m[4]*m[4]);
}
function symFrob(m) { return Math.sqrt(symFrobSq(m)); }
function symDot(a, b) {
return a[0]*b[0]+a[3]*b[3]+a[5]*b[5] + 2*(a[1]*b[1]+a[2]*b[2]+a[4]*b[4]);
}
function symScale(m, s) { return [m[0]*s, m[1]*s, m[2]*s, m[3]*s, m[4]*s, m[5]*s]; }
function symAdd(a, b) { return [a[0]+b[0], a[1]+b[1], a[2]+b[2], a[3]+b[3], a[4]+b[4], a[5]+b[5]]; }
function symSub(a, b) { return [a[0]-b[0], a[1]-b[1], a[2]-b[2], a[3]-b[3], a[4]-b[4], a[5]-b[5]]; }
// Regime centers (in log-domain, i.e., symmetric matrices)
// Normal regime: moderate variances, low correlations
var regime0 = [0.3, 0.05, 0.05, 0.35, 0.0, 0.4];
// Crisis regime: high variances, high correlations
var regime1 = [0.3 + sep, 0.05 + 0.8*sep, 0.05 + 0.7*sep, 0.35 + sep, 0.0 + 0.5*sep, 0.4 + sep];
// Generate data
var X = [], Y = [];
for (var i = 0; i < nPerRegime; i++) {
// Regime 0
var x0 = [regime0[0]+noise*randn(), regime0[1]+noise*randn(), regime0[2]+noise*randn(),
regime0[3]+noise*randn(), regime0[4]+noise*randn(), regime0[5]+noise*randn()];
X.push(x0); Y.push(0);
// Regime 1
var x1 = [regime1[0]+noise*randn(), regime1[1]+noise*randn(), regime1[2]+noise*randn(),
regime1[3]+noise*randn(), regime1[4]+noise*randn(), regime1[5]+noise*randn()];
X.push(x1); Y.push(1);
}
var n = X.length;
// Fréchet mean (anchor) = sample mean in log-domain
var muHat = [0,0,0,0,0,0];
for (var i = 0; i < n; i++) {
muHat = symAdd(muHat, X[i]);
}
muHat = symScale(muHat, 1/n);
// Fit logistic regression: max L_n(beta) = (1/n) sum [Y_i * <X_i-mu, beta-mu> - log(1+exp(...))]
// Gradient ascent in the 6-dimensional log-domain
var betaHat = symAdd(muHat, [0.01,0,0,0.01,0,0.01]); // initialize near mu
var lr = 0.05, nIter = 500;
for (var iter = 0; iter < nIter; iter++) {
var grad = [0,0,0,0,0,0];
var ll = 0;
for (var i = 0; i < n; i++) {
var h = symDot(symSub(X[i], muHat), symSub(betaHat, muHat));
var p = 1 / (1 + Math.exp(-h));
ll += Y[i] * h - Math.log(1 + Math.exp(h));
var resid = Y[i] - p;
grad = symAdd(grad, symScale(symSub(X[i], muHat), resid));
}
grad = symScale(grad, 1/n);
betaHat = symAdd(betaHat, symScale(grad, lr));
// Normalize step size decay
if (iter % 100 === 99) lr *= 0.5;
}
// Compute log-likelihood at betaHat
var llFinal = 0;
for (var i = 0; i < n; i++) {
var h = symDot(symSub(X[i], muHat), symSub(betaHat, muHat));
llFinal += Y[i] * h - Math.log(1 + Math.exp(h));
}
llFinal /= n;
// Classification accuracy
var correct = 0;
var probs = [];
for (var i = 0; i < n; i++) {
var h = symDot(symSub(X[i], muHat), symSub(betaHat, muHat));
var pred = h >= 0 ? 1 : 0;
if (pred === Y[i]) correct++;
probs.push(1 / (1 + Math.exp(-h)));
}
var accuracy = correct / n;
// The discriminating direction (in log-domain)
var discDir = symSub(betaHat, muHat);
var discNorm = symFrob(discDir);
// Project data onto the discriminating direction for 1D visualization
var projections = [];
for (var i = 0; i < n; i++) {
var centered = symSub(X[i], muHat);
if (discNorm > 1e-10) {
projections.push(symDot(centered, discDir) / discNorm);
} else {
projections.push(0);
}
}
// ROC data: sort by predicted prob
var idx = Array.from({length: n}, function(_, k) { return k; });
idx.sort(function(a, b) { return probs[a] - probs[b]; });
var tpr = [], fpr = [], auc = 0;
var n0 = 0, n1 = 0;
for (var i = 0; i < n; i++) { if (Y[i] === 0) n0++; else n1++; }
var tp = 0, fp = 0, prevFpr = 0;
for (var k = n - 1; k >= 0; k--) {
var i = idx[k];
if (Y[i] === 1) {
tp++;
} else {
fp++;
auc += (tp / n1) * ((fp / n0) - prevFpr);
prevFpr = fp / n0;
}
tpr.push(tp / n1);
fpr.push(fp / n0);
}
auc += (tp / n1) * (1 - prevFpr);
return {
X: X, Y: Y, n: n, muHat: muHat, betaHat: betaHat,
discDir: discDir, discNorm: discNorm,
projections: projections, probs: probs,
accuracy: accuracy, auc: auc, llFinal: llFinal,
regime0: regime0, regime1: regime1
};
}
pcResult = runPortfolioClassify(pc_n, pc_sep, pc_noise, pc_seed);
html`
<div style="font-family: system-ui, sans-serif; max-width: 920px;">
<div style="display: flex; gap: 20px; flex-wrap: wrap;">
<!-- Projection histogram -->
<div>
<h4>Data Projected onto Discriminating Direction</h4>
<svg width="420" height="280" style="border: 1px solid #dee2e6; border-radius: 4px;">
${(() => {
var res = pcResult;
var margin = {top: 20, right: 15, bottom: 35, left: 50};
var pw = 420 - margin.left - margin.right;
var ph = 280 - margin.top - margin.bottom;
// Separate projections by class
var proj0 = [], proj1 = [];
for (var i = 0; i < res.n; i++) {
if (res.Y[i] === 0) proj0.push(res.projections[i]);
else proj1.push(res.projections[i]);
}
var allP = res.projections;
var pMin = Math.min.apply(null, allP), pMax = Math.max.apply(null, allP);
var pRange = pMax - pMin || 0.01;
var pad = pRange * 0.1;
pMin -= pad; pMax += pad; pRange = pMax - pMin;
function xS(v) { return margin.left + ((v - pMin) / pRange) * pw; }
// Histogram bins
var nBins = 40;
var bins0 = new Array(nBins).fill(0);
var bins1 = new Array(nBins).fill(0);
for (var i = 0; i < proj0.length; i++) {
var bi = Math.floor((proj0[i] - pMin) / pRange * nBins);
if (bi >= 0 && bi < nBins) bins0[bi]++;
}
for (var i = 0; i < proj1.length; i++) {
var bi = Math.floor((proj1[i] - pMin) / pRange * nBins);
if (bi >= 0 && bi < nBins) bins1[bi]++;
}
var maxCount = Math.max(Math.max.apply(null, bins0), Math.max.apply(null, bins1));
if (maxCount === 0) maxCount = 1;
function yS(c) { return margin.top + ph - (c / maxCount) * (ph * 0.85); }
var parts = [];
parts.push('<rect x="' + margin.left + '" y="' + margin.top + '" width="' + pw +
'" height="' + ph + '" fill="#fafafa" rx="2"/>');
// Decision boundary (projection = 0)
var dbX = xS(0);
parts.push('<line x1="' + dbX + '" y1="' + margin.top + '" x2="' + dbX +
'" y2="' + (margin.top+ph) + '" stroke="#212529" stroke-width="1.5" stroke-dasharray="6,3"/>');
parts.push('<text x="' + dbX + '" y="' + (margin.top-5) + '" text-anchor="middle" font-size="9" fill="#212529">Decision boundary</text>');
// Bars
var barW = pw / nBins - 1;
for (var b = 0; b < nBins; b++) {
if (bins0[b] > 0) {
var bx = xS(pMin + (b + 0.5) * (pRange / nBins)) - barW/2;
var bh = (bins0[b] / maxCount) * ph * 0.85;
parts.push('<rect x="' + bx + '" y="' + (margin.top+ph-bh) + '" width="' + barW +
'" height="' + bh + '" fill="#1971c2" opacity="0.7"/>');
}
if (bins1[b] > 0) {
var bx2 = xS(pMin + (b + 0.5) * (pRange / nBins)) - barW/2;
var bh2 = (bins1[b] / maxCount) * ph * 0.85;
parts.push('<rect x="' + bx2 + '" y="' + (margin.top+ph-bh2) + '" width="' + barW +
'" height="' + bh2 + '" fill="#e03131" opacity="0.7"/>');
}
}
// Axes
parts.push('<line x1="' + margin.left + '" y1="' + (margin.top+ph) + '" x2="' + (margin.left+pw) +
'" y2="' + (margin.top+ph) + '" stroke="#adb5bd"/>');
parts.push('<line x1="' + margin.left + '" y1="' + margin.top + '" x2="' + margin.left +
'" y2="' + (margin.top+ph) + '" stroke="#adb5bd"/>');
parts.push('<text x="' + (margin.left+pw/2) + '" y="' + (margin.top+ph+28) +
'" text-anchor="middle" font-size="10">Projection onto ⟨β̂ − μ̂⟩</text>');
parts.push('<text x="' + (margin.left-44) + '" y="' + (margin.top+ph/2) +
'" text-anchor="middle" font-size="10" transform="rotate(-90,' + (margin.left-44) + ',' + (margin.top+ph/2) + ')">Count</text>');
// Legend
parts.push('<rect x="' + (margin.left+10) + '" y="' + (margin.top+5) + '" width="10" height="10" fill="#1971c2" opacity="0.7"/>');
parts.push('<text x="' + (margin.left+24) + '" y="' + (margin.top+13) + '" font-size="9" fill="#616161">Regime 0 (Normal)</text>');
parts.push('<rect x="' + (margin.left+150) + '" y="' + (margin.top+5) + '" width="10" height="10" fill="#e03131" opacity="0.7"/>');
parts.push('<text x="' + (margin.left+164) + '" y="' + (margin.top+13) + '" font-size="9" fill="#616161">Regime 1 (Crisis)</text>');
parts.push('<rect x="' + margin.left + '" y="' + margin.top + '" width="' + pw +
'" height="' + ph + '" fill="none" stroke="#dee2e6"/>');
return parts.join('\n');
})()}
</svg>
</div>
<!-- ROC curve -->
<div>
<h4>ROC Curve (In-sample)</h4>
<svg width="260" height="280" style="border: 1px solid #dee2e6; border-radius: 4px;">
${(() => {
var res = pcResult;
var margin = {top: 20, right: 15, bottom: 35, left: 45};
var pw = 260 - margin.left - margin.right;
var ph = 280 - margin.top - margin.bottom;
// Recompute ROC
var idx = Array.from({length: res.n}, function(_, k) { return k; });
idx.sort(function(a, b) { return res.probs[a] - res.probs[b]; });
var n0 = 0, n1 = 0;
for (var i = 0; i < res.n; i++) { if (res.Y[i] === 0) n0++; else n1++; }
var rocX = [0], rocY = [0];
var tp = 0, fp = 0;
for (var k = res.n - 1; k >= 0; k--) {
var i = idx[k];
if (res.Y[i] === 1) tp++; else fp++;
rocX.push(fp / n0);
rocY.push(tp / n1);
}
function xS(v) { return margin.left + v * pw; }
function yS(v) { return margin.top + ph - v * ph; }
var parts = [];
parts.push('<rect x="' + margin.left + '" y="' + margin.top + '" width="' + pw +
'" height="' + ph + '" fill="#fafafa" rx="2"/>');
// Diagonal
parts.push('<line x1="' + margin.left + '" y1="' + (margin.top+ph) + '" x2="' + (margin.left+pw) +
'" y2="' + margin.top + '" stroke="#adb5bd" stroke-width="1" stroke-dasharray="4,4"/>');
// ROC curve
var d = '';
for (var j = 0; j < rocX.length; j++) {
d += (j===0?'M':'L') + ' ' + xS(rocX[j]).toFixed(1) + ' ' + yS(rocY[j]).toFixed(1);
}
parts.push('<path d="' + d + '" fill="none" stroke="#1971c2" stroke-width="2.2"/>');
// Axes
parts.push('<line x1="' + margin.left + '" y1="' + (margin.top+ph) + '" x2="' + (margin.left+pw) +
'" y2="' + (margin.top+ph) + '" stroke="#adb5bd"/>');
parts.push('<line x1="' + margin.left + '" y1="' + margin.top + '" x2="' + margin.left +
'" y2="' + (margin.top+ph) + '" stroke="#adb5bd"/>');
parts.push('<text x="' + (margin.left+pw/2) + '" y="' + (margin.top+ph+28) +
'" text-anchor="middle" font-size="10">FPR</text>');
parts.push('<text x="' + (margin.left-38) + '" y="' + (margin.top+ph/2) +
'" text-anchor="middle" font-size="10" transform="rotate(-90,' + (margin.left-38) + ',' + (margin.top+ph/2) + ')">TPR</text>');
// AUC label
parts.push('<text x="' + (margin.left+pw*0.7) + '" y="' + (margin.top+ph*0.8) +
'" font-size="11" fill="#1971c2" font-weight="bold">AUC = ' + res.auc.toFixed(3) + '</text>');
parts.push('<rect x="' + margin.left + '" y="' + margin.top + '" width="' + pw +
'" height="' + ph + '" fill="none" stroke="#dee2e6"/>');
return parts.join('\n');
})()}
</svg>
</div>
</div>
<!-- Summary stats -->
<div style="display: flex; gap: 15px; flex-wrap: wrap; margin-top: 12px;">
<div style="padding: 10px 16px; background: #f8f9fa; border-radius: 6px; font-size: 0.9em;">
<table style="border-collapse: collapse;">
<tr><td style="padding: 2px 8px;">Classification accuracy</td>
<td style="padding: 2px 8px; text-align: right; font-weight: bold;">${(pcResult.accuracy * 100).toFixed(1)}%</td></tr>
<tr><td style="padding: 2px 8px;">AUC</td>
<td style="padding: 2px 8px; text-align: right; font-weight: bold;">${pcResult.auc.toFixed(3)}</td></tr>
<tr><td style="padding: 2px 8px;">Log-likelihood</td>
<td style="padding: 2px 8px; text-align: right; font-family: monospace;">${pcResult.llFinal.toFixed(4)}</td></tr>
<tr><td style="padding: 2px 8px;">‖β̂ − μ̂‖_F</td>
<td style="padding: 2px 8px; text-align: right; font-family: monospace;">${pcResult.discNorm.toFixed(3)}</td></tr>
<tr><td style="padding: 2px 8px;">Sample size</td>
<td style="padding: 2px 8px; text-align: right;">${pcResult.n}</td></tr>
</table>
</div>
<div style="padding: 10px 16px; background: #f8f9fa; border-radius: 6px; font-size: 0.85em;">
<b>Discriminating direction</b> (log-domain, upper triangle):<br>
<span style="font-family: monospace;">
Δ = [${pcResult.discDir[0].toFixed(3)}, ${pcResult.discDir[1].toFixed(3)}, ${pcResult.discDir[2].toFixed(3)};<br>
${pcResult.discDir[1].toFixed(3)}, ${pcResult.discDir[3].toFixed(3)}, ${pcResult.discDir[4].toFixed(3)};<br>
${pcResult.discDir[2].toFixed(3)}, ${pcResult.discDir[4].toFixed(3)}, ${pcResult.discDir[5].toFixed(3)}]
</span>
</div>
</div>
</div>
`
```
::: {.callout-note title="How to read the plots"}
The histogram shows data projected onto the estimated discriminating direction $\langle \hat{\beta} - \hat{\mu} \rangle$. Larger positive values indicate higher predicted probability of crisis (Regime 1). The decision boundary (dashed vertical line) separates the two regimes. The off-diagonal entries of the discriminating direction matrix reveal which **asset pairs** drive the classification — analogously to how the fMRI application identified which brain-region connections distinguish tasks.
:::
::: {.callout-tip title="Try these experiments"}
- **Increase regime separation:** Larger separation makes the two regimes more distinct — observe how the histogram separates and accuracy rises.
- **Increase noise:** Higher noise makes the distributions overlap more; the decision boundary becomes less reliable and AUC drops.
- **Inspect the discriminating direction:** The off-diagonal entries (positions 1, 2, 4 in the compact storage) capture which asset pairs change most across regimes. In a real application, this tells the risk manager *which correlations* are most diagnostic of a crisis.
- **Compare to the fMRI application:** Both applications use the same model — an SPD-matrix covariate and a binary response — but the interpretation shifts from "which brain connections" to "which asset correlations" distinguish the classes.
:::
### Connection to the Portfolio Stress-Testing Framework
The metric-space logistic classifier complements the Fréchet regression framework developed in earlier lectures:
| Method | Goal | Predictor | Response |
|--------|------|-----------|----------|
| Global Fréchet (L6) | Estimate $\mathbb{E}[\Sigma \mid \text{stress}]$ | Stress level (scalar) | Covariance matrix (SPD) |
| Kernel/Local-linear (L7–8) | Smooth estimate of $\mathbb{E}[\Sigma \mid \text{stress}]$ | Stress level (scalar) | Covariance matrix (SPD) |
| TV-regularized (L9) | Piecewise-constant $\mathbb{E}[\Sigma \mid \text{time}]$ | Time (ordered scalar) | Covariance matrix (SPD) |
| **Metric logistic (L10)** | **Classify regime from Σ** | **Covariance matrix (SPD)** | **Binary regime label** |
The last row reverses the roles: instead of predicting the covariance matrix from a scalar, we predict a binary label from the covariance matrix. In practice, a risk manager might use **both** approaches:
1. Fit a smooth Fréchet regression model to estimate how covariance evolves with a continuous stress indicator.
2. Fit the metric-space logistic classifier to identify which features of the covariance matrix are most diagnostic of crisis vs. calm regimes.
3. Use the discriminating direction from (2) to construct a **monitoring index** that tracks, in real time, how far the current covariance matrix is from the anchor in the crisis direction.
## Asymptotic Theory
We now present the formal asymptotic theory of the metric-space logistic estimator, following the framework of @LinLin2025Binary. The analysis proceeds from five assumptions — two geometric, two distributional, and one complexity-theoretic — together with bounds on estimation and classification risk.
::: {.callout-note title="Notation"}
For two non-negative sequences $\{a_n\}$ and $\{b_n\}$, we write $a_n \lesssim b_n$ (resp. $a_n \gtrsim b_n$) if there exists a constant $c > 0$, independent of $n$, such that $a_n \le c b_n$ (resp. $a_n \ge c b_n$) for all sufficiently large $n$. We write $a_n \asymp b_n$ if both $a_n \lesssim b_n$ and $a_n \gtrsim b_n$ hold.
:::
### Assumptions
The first two assumptions are geometric; the next two are distributional.
::: {#as-space .callout-note title="Assumption 1: Geometric regularity"}
$(\mathcal{M}, d)$ is a complete uniquely geodesic metric space.
:::
This standard condition ensures that any pair of points in $\mathcal{M}$ can be connected by a unique geodesic.
::: {#as-lip .callout-note title="Assumption 2: Hölder continuity of Alexandrov projections"}
There exist constants $C_U > 0$ and $\alpha_U \in (0, 1]$ such that, for all $p, q, r, u \in \mathcal{M}$,
$$
\bigl|d(p, q)\cos\angle_p(u, q) - d(p, r)\cos\angle_p(u, r)\bigr| \le C_U\, d(q, r)^{\alpha_U}.
$$
:::
Assumption 2 is a curvature regularity condition analogous to the weak quadruple inequality in [@schotz2019convergence]. It is satisfied with $\alpha_U = 1$ in all Hadamard spaces and in certain subspaces of positively curved Alexandrov spaces. Concrete examples include the SPD manifold with the Log-Cholesky or affine-invariant metric, the positive orthant of the unit sphere, the Billera--Holmes--Vogtmann tree space, and the 2-Wasserstein space of distributions on a compact interval [see Example 1 of @LinLin2025Binary].
::: {#as-existence .callout-note title="Assumption 3: Existence of basic estimators"}
The distribution of $X \in \mathcal{M}$ has bounded support. The population Fréchet mean $\mu^*$ and the sample Fréchet mean $\hat{\mu}$ exist and are unique, the latter almost surely. Moreover, the maximum likelihood estimator $\hat{\beta}$ exists almost surely.
:::
::: {#as-margin .callout-note title="Assumption 4: Margin / identification condition"}
There exist constants $\lambda_X > 0$, $C_L > 0$, and $\alpha_L > 1 + 2/m$ for some $m \ge 4$ with $\mathbb{E}[d^m(\mu^*, X)] < \infty$, such that for all $\beta \in \mathcal{M}$,
$$
\mathbb{E}\bigl[|h(\beta; X, \mu^*) - h(\beta^*; X, \mu^*)|^2\bigr] \ge C_L \lambda_X\, d(\beta, \beta^*)^{\alpha_L}.
$$
:::
Assumption 4 controls the behavior of $L_n - L$ near $\beta^*$ and is standard in M-estimation theory [@ahidar2020convergence]. In Euclidean space $\mathcal{M} = \mathbb{R}^D$, it reduces to the familiar requirement that $\mathbb{E}[(X - \mu^*)(X - \mu^*)^{\top}]$ be positive definite, with $\lambda_X$ being its smallest eigenvalue. The moment condition $m \ge 4$ is needed for the expected *squared* estimation error; for the first-moment error, $m \ge 2$ suffices. When $\mathcal{M}$ is bounded, the constraint $\alpha_L > 1 + 2/m$ can be relaxed to $\alpha_L > 1$, since $\mathbb{E}[d^m(\mu^*, X)] < \infty$ holds for any arbitrarily large fixed $m$.
::: {.callout-note title="Discussion of existence and uniqueness"}
Some remarks on Assumptions 3–4 (see Web Appendix E of @LinLin2025Binary for details):
- **Fréchet means:** In Hadamard spaces, $\mu^*$ and $\hat{\mu}$ exist and are unique whenever $\mathbb{E}[d^2(X, u)] < \infty$ for some $u \in \mathcal{M}$, regardless of the distribution of $X$ [@sturm2003probability]. For Riemannian manifolds, sufficient conditions involve data residing in a geodesically convex region [@Bhattacharya2003; @Afsari2011; @lin2021total].
- **Existence of $\hat{\beta}$:** Proposition 1 of @LinLin2025Binary shows that $\hat{\beta}$ exists when $\mathcal{M}$ is complete and totally bounded. For unbounded spaces — even in $\mathbb{R}^D$ — additional distributional conditions are required, such as the overlap condition [@albert1984existence]. In the SPD case with the Log-Euclidean metric, an analogous overlap condition can be formulated.
- **Uniqueness of $\beta^*$:** Proposition 2 of @LinLin2025Binary establishes uniqueness of $\beta^*$ under Assumption 4.
- **Measurability:** The theoretical results hold for *any* measurable choice of $\hat{\beta}$ maximizing $L_n$; uniqueness of $\hat{\beta}$ is not required.
:::
### Local Metric Entropy
The convergence rate of $\hat{\beta}$ depends on the complexity of the metric space near $\beta^*$, as measured by covering numbers. Let $B(\beta^*, \delta) \subset \mathcal{M}$ be the ball of radius $\delta > 0$ centered at $\beta^*$, and let $N(\varepsilon, B(\beta^*, \delta), d)$ be its covering number using balls of radius $\varepsilon > 0$ [@wainwright2019high, Section 5.1]. The *metric entropy* is $\log N(\varepsilon, B(\beta^*, \delta), d)$.
::: {#as-entropy .callout-note title="Assumption 5: Local metric entropy"}
There exists a constant $D > 0$, such that one of the following holds:
(a) **(Log-polynomial metric entropy)**
$$
N(\varepsilon, B(\beta^*, \delta), d) \lesssim \Bigl(\frac{\delta}{\varepsilon}\Bigr)^{D}.
$$
(b) **(Polynomial metric entropy)**
$$
\log N(\varepsilon, B(\beta^*, \delta), d) \lesssim \Bigl(\frac{\delta}{\varepsilon}\Bigr)^{D}.
$$
:::
Part (a) is typical of **finite-dimensional** metric spaces: it holds for all finite-dimensional compact Riemannian manifolds, and more generally for spaces admitting an Ahlfors--David regular measure. In this regime the entropy grows as $D \log(\delta / \varepsilon)$. Part (b) is much weaker and accommodates spaces of much higher intrinsic complexity. A notable example is the 2-Wasserstein space of distributions on a compact interval, which satisfies (b) with $D = 1$ but not (a) [@yang1999minimax; @petersen2019frechet].
### Estimation Rates
Recall that $\alpha_U$ (from Assumption 2) is related to curvature, while $\alpha_L$ (from Assumption 4) controls the local behavior of the loss near $\beta^*$. The **benchmark rate** $v_{L,n}$ is defined piecewise according to the entropy regime:
$$
\begin{aligned}
&\text{Under Assumption 5(a) — Log-polynomial entropy:}
& v_{L,n} &= \Bigl(\frac{D}{n}\Bigr)^{\frac{1}{2(\alpha_L - \alpha_U)}}. \\[12pt]
&\text{Under Assumption 5(b) — Polynomial entropy:}
& v_{L,n} &=
\begin{cases}
\Bigl(\dfrac{1}{n}\Bigr)^{\frac{1}{2(\alpha_L - \alpha_U)}}, & \text{if } D < 2\alpha_U, \\[10pt]
\Bigl(\dfrac{\log^2 n}{n}\Bigr)^{\frac{1}{2(\alpha_L - \alpha_U)}}, & \text{if } D = 2\alpha_U, \\[10pt]
\Bigl(\dfrac{1}{n}\Bigr)^{\frac{\alpha_U}{D(\alpha_L - \alpha_U)}}, & \text{if } D > 2\alpha_U.
\end{cases}
\end{aligned}
$$
The three-case structure under polynomial entropy reflects a competition between the intrinsic dimension $D$ and the curvature exponent $\alpha_U$:
- **$D < 2\alpha_U$ (low-dimensional regime):** The rate is dimension-free — the estimator achieves the same rate as when the space is low-dimensional.
- **$D = 2\alpha_U$ (critical regime):** A logarithmic penalty appears, reflecting the boundary behavior of the entropy integral.
- **$D > 2\alpha_U$ (high-dimensional regime):** The dimension $D$ enters the exponent, degrading the rate; the estimator pays a genuine curse-of-dimensionality cost relative to the log-polynomial case.
::: {.callout-note title="The common case: $\alpha_L = 2\alpha_U$"}
An important and common setting is $\alpha_L = 2\alpha_U$. For example, on a Hadamard manifold where the expected Hessian of $h$ at $\beta^*$ is positive definite, one has $\alpha_L = 2$ and $\alpha_U = 1$. In this case, under log-polynomial entropy, $v_{L,n} = \sqrt{D/n}$, matching the familiar parametric rate in dimension $D$.
:::
::: {#thm-logistic-rates .theorem title="Estimation rate — Upper bound"}
Suppose Assumptions 1–5 hold and the metric-space logistic model is well specified. Then there exists a constant $c > 0$, depending on the constants in the assumptions, such that for any measurable sequence $\hat{\beta}$ determined by $\hat{\beta} \in \arg\max_{\beta \in \mathcal{M}} L_n(\beta)$,
$$
\mathbb{E}\bigl[d^2(\beta^*, \hat{\beta})\bigr] \le c\, v_{L,n}^2.
$$
:::
::: {#thm-logistic-lower .theorem title="Estimation rate — Minimax lower bound"}
Suppose $\sigma_X^2 = \mathbb{E}[d^2(\mu^*, X)] < \infty$ and Assumptions 1–2 hold. If $(X_i, Y_i)_{i=1}^n$ are i.i.d. from the metric-space logistic model, then
$$
\inf_{\hat{\beta}} \sup_{\beta \in \mathcal{M}} \mathbb{E}_\beta\bigl[d^2(\beta, \hat{\beta})\bigr] \gtrsim
\begin{cases}
\Bigl(\dfrac{D}{\sigma_X^2 n}\Bigr)^{\frac{1}{\alpha_U}}, & \text{when } N(\varepsilon, B(\beta^*, \delta), d) \asymp \Bigl(\dfrac{\delta}{\varepsilon}\Bigr)^{D}, \\[10pt]
\Bigl(\dfrac{1}{\sigma_X^2 n}\Bigr)^{\frac{1}{\alpha_U}}, & \text{when } \log N(\varepsilon, B(\beta^*, \delta), d) \asymp \Bigl(\dfrac{\delta}{\varepsilon}\Bigr)^{D}.
\end{cases}
$$
Here the infimum is taken over all estimators $\hat{\beta}$, and $\mathbb{E}_\beta$ denotes the expectation under the model with $\beta^* = \beta$.
:::
**Optimality summary.** When the common quadratic-margin relation $\alpha_L = 2\alpha_U$ holds (as is the case on Hadamard manifolds with a positive-definite expected Hessian):
- Under log-polynomial entropy (Assumption 5a), the upper and lower bounds match, establishing that $\mathbb{E}[d^2(\beta^*, \hat{\beta})] \asymp D/n$. The estimator is **rate-optimal**.
- Under polynomial entropy (Assumption 5b) with $D \le 2\alpha_U$, the estimator is also rate-optimal, up to a logarithmic factor when $D = 2\alpha_U$.
- When $D > 2\alpha_U$, a gap between the upper and lower bound remains.
### Classification Risk
The regression model induces a plug-in classifier. Defining the estimated success probability
$$
\hat{P}(x) = \bigl[1 + \exp\{-h(\hat{\beta}; x, \hat{\mu})\}\bigr]^{-1},
$$
the classifier is
$$
\hat{s}(x) = \begin{cases} 1, & \hat{P}(x) \ge 1/2, \\ 0, & \hat{P}(x) < 1/2. \end{cases}
$$
The Bayes classifier $s^*$ — minimizing the error probability under the model — is defined analogously with $(\mu^*, \beta^*)$ in place of $(\hat{\mu}, \hat{\beta})$. The **excess risk** of $\hat{s}$ is
$$
\mathcal{E}(\hat{s}, s^*) = \mathbb{E}[\ell(\hat{s})] - \mathbb{E}[\ell(s^*)],
\qquad \ell(s) = \mathbf{1}\{s(X) \ne Y\}.
$$
The excess risk measures how much accuracy is lost relative to the optimal (Bayes) rule.
::: {#thm-classification-metric .theorem title="Classification risk — Alexandrov space bound"}
Suppose $\mathcal{M}$ is a complete Alexandrov space with curvature upper bounded by $\kappa$ and diameter at most $D_\kappa/2$, where $D_\kappa = \pi/\sqrt{\kappa}$ if $\kappa > 0$ and $D_\kappa = \infty$ otherwise. Assume further that the Alexandrov inner product $h(\beta^*; x, \mu)$ is regular near $\mu^*$ (i.e., Lipschitz in both arguments locally). Then, under Assumption 3, there exists a constant $c > 0$ such that
$$
\mathcal{E}(\hat{s}, s^*) \le c\,
\Bigl(\mathbb{E}[d^2(\mu^*, \hat{\mu})] + \mathbb{E}[d^2(\beta^*, \hat{\beta})]\Bigr)^{1/2}.
$$
:::
Under standard regularity conditions, the Fréchet mean converges at rate $\mathbb{E}[d^2(\mu^*, \hat{\mu})] \lesssim n^{-1}$ [@schotz2019convergence; @ahidar2020convergence], which is dominated by the estimation error of $\hat{\beta}$. Hence the excess risk inherits the rates established above for $\hat{\beta}$.
::: {#thm-classification-manifold .theorem title="Classification risk — Riemannian manifold bound and optimality"}
For a $D$-dimensional Hadamard manifold of bounded sectional curvature, suppose the distribution of $X$ has bounded support, $\lambda_X > 0$, and $\hat{\beta}$ exists almost surely. Then there exists a constant $c > 0$ such that
$$
\mathcal{E}(\hat{s}, s^*) \le c\,\sqrt{\frac{D}{n}}.
$$
The same conclusion holds for a complete, simply connected $D$-dimensional Riemannian manifold with non-negative sectional curvature bounded above by $\kappa > 0$ and diameter at most $D_\kappa/2$, under the same Hessian non-degeneracy condition as in the estimation corollary.
Moreover, this rate is **minimax optimal**: for any $D$-dimensional Riemannian manifold $\mathcal{M}$,
$$
\inf_{\widehat{s}} \sup_{s^* \in \mathcal{C}} \mathcal{E}(\hat{s}, s^*) \gtrsim \sqrt{\frac{D}{n}},
$$
where $\mathcal{C} = \{s(X) = \mathbf{1}\{h(\beta; X, \mu^*) \ge 0\} : \beta \in \mathcal{M}\}$ is the class of half-space-type classifiers induced by the model, and the infimum is taken over all classifiers.
:::
::: {.callout-tip title="Interpretation of the $\sqrt{D/n}$ rate"}
The classification risk bound $\mathcal{E}(\hat{s}, s^*) \lesssim \sqrt{D/n}$ reveals three important features:
1. **Parametric rate.** The $n^{-1/2}$ dependence is optimal even in $\mathbb{R}^D$, confirming that the metric-space logistic classifier does not pay a nonparametric penalty.
2. **Dimension dependence.** The $\sqrt{D}$ factor reflects the intrinsic dimension of the covariate space — not an ambient or embedding dimension — and matches the classical Euclidean result.
3. **Minimax optimality.** The matching lower bound means that no classifier can achieve a faster rate uniformly over the class $\mathcal{C}$, regardless of whether it exploits metric, Riemannian, or Euclidean structure. The proposed estimator attains the fundamental limit.
:::
## Interactive Exploration: Logistic Regression on the 2D Sphere
The following demo implements metric-space logistic regression on the 2D sphere $\mathbb{S}^2$, where the covariate $X$ lies on the unit sphere in $\mathbb{R}^3$ and the response is binary. Compared to the circle $S^1$, the sphere offers a richer geometry: the decision boundary is a **great circle** orthogonal to the discriminating direction at the anchor, and the estimated parameters $\hat{\mu}$ and $\hat{\beta}$ are visualized directly on the sphere.
The **class separation** slider directly sets the geodesic distance between the two class centers. Larger values move the red and blue point clouds apart across the true great-circle boundary; smaller values move both clouds toward the boundary and create more overlap.
**Visual guide:**
- **Wireframe sphere** = $\mathbb{S}^2$, with the north pole at top
- **Red circles** = $Y = 1$, **Blue triangles** = $Y = 0$
- **Green sphere** = estimated anchor $\hat{\mu}$ (Fréchet mean)
- **Blue arrow** = estimated discriminating direction $\hat{\beta}$
- **Dashed red great circle** = estimated decision boundary (where $h(\hat{\beta}; x, \hat{\mu}) = 0$)
- **Dashed green great circle** = true decision boundary
```{ojs}
//| label: fig-metric-logistic-s2
//| fig-cap: "Interactive: Metric-space logistic regression on the 2D sphere S²"
s2_n_control = Inputs.range([30, 300], {step: 10, value: 300, label: "Sample size n"})
s2_separation_control = Inputs.range([10, 120], {step: 5, value: 60, label: "Class separation (geodesic °)"})
s2_conc_control = Inputs.range([2, 30], {step: 1, value: 3, label: "Covariate concentration κ"})
s2_rotY_control = Inputs.range([0, 360], {step: 5, value: 210, label: "View rotation (azimuth °)"})
s2_rotX_control = Inputs.range([-60, 60], {step: 5, value: 25, label: "View tilt (elevation °)"})
s2_n = Generators.input(s2_n_control)
s2_separation = Generators.input(s2_separation_control)
s2_conc = Generators.input(s2_conc_control)
s2_rotY = Generators.input(s2_rotY_control)
s2_rotX = Generators.input(s2_rotX_control)
s2_controls_view = html`
<style>
.s2-slider-grid {
display: grid;
grid-template-columns: repeat(2, minmax(0, 1fr));
gap: 6px 20px;
max-width: 620px;
margin: 0 auto 12px;
font-family: system-ui, sans-serif;
font-size: 0.85em;
}
.s2-slider-grid > * {
min-width: 0;
margin: 0;
}
.s2-slider-grid input[type="number"] {
width: 7.5rem !important;
}
@media (max-width: 500px) {
.s2-slider-grid {
grid-template-columns: 1fr;
}
}
</style>
<div class="s2-slider-grid">
<div>${s2_n_control}</div>
<div>${s2_separation_control}</div>
<div>${s2_conc_control}</div>
<div>${s2_rotY_control}</div>
<div>${s2_rotX_control}</div>
</div>
`
function attachS2DragRotation(root, selector, azimuthLabel, elevationLabel, minElevation, maxElevation, data, kind) {
var svg = root.querySelector(selector);
if (!svg) return root;
var scene = svg.querySelector('g[transform^="translate"]');
if (!scene) return root;
function findRange(labelText) {
var labels = Array.from(document.querySelectorAll("label"));
for (var i = 0; i < labels.length; i++) {
if (labels[i].textContent.indexOf(labelText) >= 0) {
return labels[i].querySelector('input[type="range"]') ||
(labels[i].parentElement && labels[i].parentElement.querySelector('input[type="range"]'));
}
}
return null;
}
var az = findRange(azimuthLabel);
var el = findRange(elevationLabel);
var state = {
rotY: az ? Number(az.value) : 210,
rotX: el ? Number(el.value) : 25
};
function clamp(x, lo, hi) {
return Math.max(lo, Math.min(hi, x));
}
function makeProjector(rotY, rotX, R) {
var ry = rotY * Math.PI / 180;
var rx = rotX * Math.PI / 180;
return function(v) {
var x1 = v[0]*Math.cos(ry) + v[2]*Math.sin(ry);
var y1 = v[1];
var z1 = -v[0]*Math.sin(ry) + v[2]*Math.cos(ry);
var x2 = x1;
var y2 = y1*Math.cos(rx) - z1*Math.sin(rx);
var z2 = y1*Math.sin(rx) + z1*Math.cos(rx);
return {x: x2*R, y: -y2*R, z: z2};
};
}
function renderLogisticS2(res, rotY, rotX) {
var parts = [];
var R = 180;
var proj = makeProjector(rotY, rotX, R);
function vdot_local(a, b) { return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; }
for (var p = 0; p <= 18; p++) {
var z = 1 - p * 2/18;
var r = Math.sqrt(Math.max(0, 1 - z*z));
var d = '';
for (var j = 0; j <= 120; j++) {
var phi = j * 2 * Math.PI / 120;
var pt = proj([r*Math.cos(phi), r*Math.sin(phi), z]);
d += (j===0 ? 'M' : 'L') + ' ' + pt.x.toFixed(1) + ' ' + pt.y.toFixed(1);
}
var alpha = p === 0 || p === 12 ? 0.25 : 0.12;
parts.push('<path d="' + d + '" fill="none" stroke="#adb5bd" stroke-width="0.6" opacity="' + alpha + '"/>');
}
for (var m = 0; m < 24; m++) {
var mphi = m * Math.PI / 12;
var md = '';
for (var mj = 0; mj <= 90; mj++) {
var theta = mj * Math.PI / 90;
var mpt = proj([Math.sin(theta)*Math.cos(mphi), Math.sin(theta)*Math.sin(mphi), Math.cos(theta)]);
md += (mj===0 ? 'M' : 'L') + ' ' + mpt.x.toFixed(1) + ' ' + mpt.y.toFixed(1);
}
parts.push('<path d="' + md + '" fill="none" stroke="#adb5bd" stroke-width="0.6" opacity="0.12"/>');
}
parts.push('<circle cx="0" cy="0" r="' + R + '" fill="none" stroke="#dee2e6" stroke-width="1.2"/>');
var dPath = '';
for (var dj = 0; dj < res.decPoints.length; dj++) {
var dpt = proj(res.decPoints[dj]);
dPath += (dj===0 ? 'M' : 'L') + ' ' + dpt.x.toFixed(1) + ' ' + dpt.y.toFixed(1);
}
parts.push('<path d="' + dPath + '" fill="none" stroke="#e03131" stroke-width="2" stroke-dasharray="10,5" opacity="0.85"/>');
var tdPath = '';
for (var tj = 0; tj < res.trueDecPoints.length; tj++) {
var tdpt = proj(res.trueDecPoints[tj]);
tdPath += (tj===0 ? 'M' : 'L') + ' ' + tdpt.x.toFixed(1) + ' ' + tdpt.y.toFixed(1);
}
parts.push('<path d="' + tdPath + '" fill="none" stroke="#2b8a3e" stroke-width="1.5" stroke-dasharray="5,4" opacity="0.6"/>');
var dataPts = [];
for (var i = 0; i < res.n; i++) {
var pt = proj(res.X[i]);
dataPts.push({x: pt.x, y: pt.y, z: pt.z, cls: res.Y[i]});
}
dataPts.sort(function(a, b) { return a.z - b.z; });
for (var di = 0; di < dataPts.length; di++) {
var dp = dataPts[di];
var col = dp.cls === 1 ? '#e03131' : '#1971c2';
var opacity = 0.35 + 0.45 * ((dp.z + 1) / 2);
if (dp.cls === 1) {
parts.push('<circle cx="' + dp.x.toFixed(1) + '" cy="' + dp.y.toFixed(1) + '" r="3.5" fill="' + col + '" opacity="' + opacity.toFixed(2) + '"/>');
} else {
var s = 5, px = dp.x, py = dp.y;
parts.push('<polygon points="' + px.toFixed(1) + ',' + (py-s).toFixed(1) + ' ' +
(px-s*0.866).toFixed(1) + ',' + (py+s*0.5).toFixed(1) + ' ' +
(px+s*0.866).toFixed(1) + ',' + (py+s*0.5).toFixed(1) +
'" fill="' + col + '" opacity="' + opacity.toFixed(2) + '"/>');
}
}
var muP = proj(res.muHat);
var betaP = proj(res.betaHat);
var tb = [res.betaHat[0]-res.muHat[0]*vdot_local(res.muHat, res.betaHat),
res.betaHat[1]-res.muHat[1]*vdot_local(res.muHat, res.betaHat),
res.betaHat[2]-res.muHat[2]*vdot_local(res.muHat, res.betaHat)];
var ntb = Math.sqrt(tb[0]*tb[0]+tb[1]*tb[1]+tb[2]*tb[2]);
if (ntb > 1e-10) {
tb = [tb[0]/ntb, tb[1]/ntb, tb[2]/ntb];
var arrowLen = 0.55;
var arcPath = '';
for (var a = 0; a <= 20; a++) {
var t = a * arrowLen / 20;
var interp = [res.muHat[0] + t*tb[0], res.muHat[1] + t*tb[1], res.muHat[2] + t*tb[2]];
var ni = Math.sqrt(interp[0]*interp[0]+interp[1]*interp[1]+interp[2]*interp[2]);
interp = [interp[0]/ni, interp[1]/ni, interp[2]/ni];
var ip = proj(interp);
arcPath += (a===0 ? 'M' : 'L') + ' ' + ip.x.toFixed(1) + ' ' + ip.y.toFixed(1);
}
parts.push('<path d="' + arcPath + '" fill="none" stroke="#1971c2" stroke-width="3" stroke-linecap="round"/>');
var tip = [res.muHat[0] + arrowLen*tb[0], res.muHat[1] + arrowLen*tb[1], res.muHat[2] + arrowLen*tb[2]];
var nTip = Math.sqrt(tip[0]*tip[0]+tip[1]*tip[1]+tip[2]*tip[2]);
tip = [tip[0]/nTip, tip[1]/nTip, tip[2]/nTip];
var ah = proj(tip);
parts.push('<circle cx="' + ah.x.toFixed(1) + '" cy="' + ah.y.toFixed(1) + '" r="5.5" fill="#1971c2"/>');
}
parts.push('<circle cx="' + muP.x.toFixed(1) + '" cy="' + muP.y.toFixed(1) + '" r="6" fill="#2b8a3e" stroke="#fff" stroke-width="1.5"/>');
parts.push('<circle cx="' + betaP.x.toFixed(1) + '" cy="' + betaP.y.toFixed(1) + '" r="4.5" fill="#1971c2" stroke="#fff" stroke-width="1"/>');
var betaTrueP = proj(res.betaStarDir);
parts.push('<circle cx="' + betaTrueP.x.toFixed(1) + '" cy="' + betaTrueP.y.toFixed(1) + '" r="4" fill="#2b8a3e" stroke="#fff" stroke-width="1" opacity="0.7"/>');
return parts.join('\n');
}
function render() {
if (kind === "logistic") scene.innerHTML = renderLogisticS2(data, state.rotY, state.rotX);
}
var dragging = false;
var lastX = 0;
var 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;
var dx = event.clientX - lastX;
var dy = event.clientY - lastY;
lastX = event.clientX;
lastY = event.clientY;
state.rotY = (state.rotY + dx * 0.55 + 360) % 360;
// Match vertical rotation to the pointer's direction.
state.rotX = clamp(state.rotX + dy * 0.55, minElevation, maxElevation);
render();
});
function stopDrag(event) {
if (!dragging) return;
dragging = false;
svg.style.cursor = "grab";
if (event.pointerId !== undefined) svg.releasePointerCapture(event.pointerId);
}
svg.addEventListener("pointerup", stopDrag);
svg.addEventListener("pointercancel", stopDrag);
svg.addEventListener("pointerleave", stopDrag);
return root;
}
function runLogisticS2(n, rotYdeg, rotXdeg, kappa, separation, seed) {
// ---- PRNG ----
function mulberry32(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;
};
}
var rng = mulberry32(seed);
// ---- 3D vector utilities ----
function vec3(x, y, z) { return [x, y, z]; }
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 vsub(a, b) { return [a[0]-b[0], a[1]-b[1], a[2]-b[2]]; }
function vadd(a, b) { return [a[0]+b[0], a[1]+b[1], a[2]+b[2]]; }
function vnormalize(a) {
var n = vnorm(a);
if (n < 1e-14) return [0, 0, 1];
return [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]];
}
// ---- Spherical coordinates to 3D ----
function sphTo3D(theta, phi) {
return [Math.sin(theta)*Math.cos(phi), Math.sin(theta)*Math.sin(phi), Math.cos(theta)];
}
function toSph(v) {
var r = vnorm(v);
var theta = Math.acos(Math.max(-1, Math.min(1, v[2]/r)));
var phi = Math.atan2(v[1], v[0]);
return {theta: theta, phi: phi};
}
// Geodesic distance on unit sphere
function geoDist(a, b) {
var dot = vdot(a, b);
return Math.acos(Math.max(-1, Math.min(1, dot)));
}
// ---- Generate two spherical class clouds ----
// Their centers are symmetric about muStar, so the population anchor remains
// the north pole and the true decision boundary remains the great circle x=0.
var muStar = [0, 0, 1]; // north pole as anchor
function rvMFAround(center, kap) {
// Exact inverse-CDF sampler for the S² von Mises-Fisher distribution
// followed by a rotation from the north pole to the requested center.
var u = rng();
var z;
if (kap < 1e-8) {
z = 2*u - 1;
} else {
z = 1 + Math.log(u + (1 - u) * Math.exp(-2 * kap)) / kap;
}
var r = Math.sqrt(Math.max(0, 1 - z*z));
var phi = rng() * 2 * Math.PI;
var reference = Math.abs(center[2]) < 0.9 ? [0, 0, 1] : [0, 1, 0];
var e1 = vnormalize(vcross(reference, center));
var e2 = vnormalize(vcross(center, e1));
return vnormalize(vadd(
vadd(vscale(e1, r * Math.cos(phi)), vscale(e2, r * Math.sin(phi))),
vscale(center, z)
));
}
// Put the two centers on opposite sides of the boundary, separated by the
// slider's geodesic angle.
var halfSeparation = separation * Math.PI / 360;
var center0 = sphTo3D(halfSeparation, Math.PI);
var center1 = sphTo3D(halfSeparation, 0);
var X = [], Y = [];
for (var i = 0; i < n; i++) {
var cls = i % 2;
X.push(rvMFAround(cls === 1 ? center1 : center0, kappa));
Y.push(cls);
}
// ---- True coefficient direction ----
// β* points from the anchor toward the Y=1 class center.
var betaStarThetaDeg = 90;
var betaStarDir = sphTo3D(betaStarThetaDeg * Math.PI / 180, 0);
betaStarDir = vnormalize(betaStarDir);
// ---- Alexandrov inner product on S² ----
function alexIPS2(beta, x, anchor) {
// d(anchor, x)
var dx = geoDist(anchor, x);
var db = geoDist(anchor, beta);
if (dx < 1e-12 || db < 1e-12) return 0;
// Tangent vector at anchor toward x
var vx = vsub(x, vscale(anchor, vdot(anchor, x)));
var nvx = vnorm(vx);
if (nvx < 1e-12) return 0;
vx = vscale(vx, 1/nvx);
// Tangent vector at anchor toward beta
var vb = vsub(beta, vscale(anchor, vdot(anchor, beta)));
var nvb = vnorm(vb);
if (nvb < 1e-12) return 0;
vb = vscale(vb, 1/nvb);
var cosAngle = vdot(vx, vb);
return dx * db * cosAngle;
}
// ---- Estimate Fréchet mean (anchor) ----
// Gradient descent on S²
var muHat = [0, 0, 1];
var lrMu = 0.35;
for (var iter = 0; iter < 120; iter++) {
var descent = [0, 0, 0];
for (var i = 0; i < n; i++) {
var d = geoDist(muHat, X[i]);
if (d < 1e-12) continue;
// The negative gradient points from muHat toward each observation.
var tdir = vsub(X[i], vscale(muHat, vdot(muHat, X[i])));
var ntd = vnorm(tdir);
if (ntd < 1e-12) continue;
tdir = vscale(tdir, 1/ntd);
descent = vadd(descent, vscale(tdir, d));
}
descent = vscale(descent, 1/n);
var nd = vnorm(descent);
if (nd < 1e-8) break;
var step = Math.min(lrMu * nd, 0.25);
var dir = vscale(descent, 1/nd);
muHat = vadd(vscale(muHat, Math.cos(step)), vscale(dir, Math.sin(step)));
muHat = vnormalize(muHat);
}
// ---- Fit logistic model: grid search over S² ----
// Fibonacci sphere for uniform-ish grid
var nGrid = 300;
var betaCandidates = [];
var phi_golden = Math.PI * (3 - Math.sqrt(5));
for (var g = 0; g < nGrid; g++) {
var y = 1 - (g / (nGrid - 1)) * 2; // y from 1 to -1
var radius = Math.sqrt(1 - y * y);
var theta = phi_golden * g;
betaCandidates.push([Math.cos(theta)*radius, Math.sin(theta)*radius, y]);
}
var bestLL = -Infinity, betaHat = betaCandidates[0];
var logLikGrid = [];
for (var g = 0; g < nGrid; g++) {
var beta = betaCandidates[g];
var ll = 0;
for (var i = 0; i < n; i++) {
var h = alexIPS2(beta, X[i], muHat);
ll += Y[i] * h - Math.log(1 + Math.exp(h));
}
ll /= n;
logLikGrid.push({beta: beta, ll: ll});
if (ll > bestLL) { bestLL = ll; betaHat = beta; }
}
// ---- Classification accuracy ----
var correct = 0;
var probs = [];
for (var i = 0; i < n; i++) {
var hVal = alexIPS2(betaHat, X[i], muHat);
var pred = hVal >= 0 ? 1 : 0;
if (pred === Y[i]) correct++;
probs.push(1 / (1 + Math.exp(-hVal)));
}
var accuracy = correct / n;
// ---- Decision boundary great circle ----
// Points on S² where alexIPS2(betaHat, x, muHat) = 0
// This is the great circle orthogonal to the tangent direction of betaHat at muHat
var tb = vsub(betaHat, vscale(muHat, vdot(muHat, betaHat)));
var ntb = vnorm(tb);
var decPlaneNormal = [0, 0, 1];
if (ntb > 1e-10) {
tb = vscale(tb, 1/ntb);
// h(betaHat; x, muHat) = 0 iff x is orthogonal to this tangent direction.
// Thus the decision boundary is the great circle in the plane normal to tb.
decPlaneNormal = tb;
}
// True decision boundary
var tbTrue = vsub(betaStarDir, vscale(muStar, vdot(muStar, betaStarDir)));
var ntbTrue = vnorm(tbTrue);
var trueDecNormal = [0, 1, 0]; // default
if (ntbTrue > 1e-10) {
tbTrue = vscale(tbTrue, 1/ntbTrue);
trueDecNormal = tbTrue;
}
// Sample dec boundary points
var decPoints = [];
var nDec = 120;
for (var d = 0; d < nDec; d++) {
var angle = d * 2 * Math.PI / nDec;
// Point on great circle: rotate muHat by 90° around decPlaneNormal
var axis = decPlaneNormal;
var cosa = Math.cos(angle), sina = Math.sin(angle);
// Rodrigues rotation around axis
// Start from a point orthogonal to both muHat and decNormal
var startPt = vcross(decPlaneNormal, muHat);
startPt = vnormalize(startPt);
var pt = vadd(vadd(vscale(startPt, cosa), vscale(vcross(axis, startPt), sina)),
vscale(axis, vdot(axis, startPt) * (1 - cosa)));
decPoints.push(vnormalize(pt));
}
var trueDecPoints = [];
for (var d = 0; d < nDec; d++) {
var angle = d * 2 * Math.PI / nDec;
var startPt = vcross(trueDecNormal, muStar);
startPt = vnormalize(startPt);
var axis = trueDecNormal;
var cosa = Math.cos(angle), sina = Math.sin(angle);
var pt = vadd(vadd(vscale(startPt, cosa), vscale(vcross(axis, startPt), sina)),
vscale(axis, vdot(axis, startPt) * (1 - cosa)));
trueDecPoints.push(vnormalize(pt));
}
// ---- Rotation matrices ----
var ry = rotYdeg * Math.PI / 180;
var rx = rotXdeg * Math.PI / 180;
function rotate(v) {
// Rotate around Y axis
var x1 = v[0]*Math.cos(ry) + v[2]*Math.sin(ry);
var y1 = v[1];
var z1 = -v[0]*Math.sin(ry) + v[2]*Math.cos(ry);
// Rotate around X axis
var x2 = x1;
var y2 = y1*Math.cos(rx) - z1*Math.sin(rx);
var z2 = y1*Math.sin(rx) + z1*Math.cos(rx);
return [x2, y2, z2];
}
function project(v) {
// Orthographic projection: (x, y) with z for depth
return {x: v[0], y: v[1], z: v[2]};
}
return {
X: X, Y: Y, probs: probs,
muStar: muStar, betaStarDir: betaStarDir,
muHat: muHat, betaHat: betaHat,
decPoints: decPoints, trueDecPoints: trueDecPoints,
decPlaneNormal: decPlaneNormal, trueDecNormal: trueDecNormal,
betaStarThetaDeg: betaStarThetaDeg,
accuracy: accuracy, bestLL: bestLL,
n: n, kappa: kappa, separation: separation,
center0: center0, center1: center1,
rotate: rotate, project: project
};
}
s2result = runLogisticS2(s2_n, s2_rotY, s2_rotX, s2_conc, s2_separation, 42);
(() => {
const view = html`
<div style="font-family: system-ui, sans-serif; max-width: 920px;">
<div style="display: flex; gap: 20px; flex-wrap: wrap;">
<!-- Sphere visualization -->
<div>
<h4>Classification on 𝕊²</h4>
<svg class="s2-logistic-globe" width="440" height="440" style="border: 1px solid #dee2e6; border-radius: 4px; cursor: grab; touch-action: none; user-select: none; background: linear-gradient(180deg, #fbfdff 0%, #f5f8fb 100%);">
<g transform="translate(220, 220)">
${(() => {
var res = s2result;
var parts = [];
var R = 180;
// Helper: dot product (defined inside this scope)
function vdot_local(a, b) { return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; }
// Helper: project + depth sort
function proj(v) {
var rv = res.rotate(v);
return {x: rv[0]*R, y: -rv[1]*R, z: rv[2]};
}
// Sphere wireframe: parallels and meridians
// Parallels (constant z)
for (var p = 0; p <= 18; p++) {
var z = 1 - p * 2/18;
var r = Math.sqrt(Math.max(0, 1 - z*z));
var d = '';
var nPts = 120;
for (var j = 0; j <= nPts; j++) {
var phi = j * 2 * Math.PI / nPts;
var pt = proj([r*Math.cos(phi), r*Math.sin(phi), z]);
d += (j===0 ? 'M' : 'L') + ' ' + pt.x.toFixed(1) + ' ' + pt.y.toFixed(1);
}
var alpha = p === 0 || p === 12 ? 0.25 : 0.12;
parts.push('<path d="' + d + '" fill="none" stroke="#adb5bd" stroke-width="0.6" opacity="' + alpha + '"/>');
}
// Meridians (constant phi)
for (var m = 0; m < 24; m++) {
var phi = m * Math.PI / 12;
var d = '';
var nPts = 90;
for (var j = 0; j <= nPts; j++) {
var theta = j * Math.PI / nPts;
var pt = proj([Math.sin(theta)*Math.cos(phi), Math.sin(theta)*Math.sin(phi), Math.cos(theta)]);
d += (j===0 ? 'M' : 'L') + ' ' + pt.x.toFixed(1) + ' ' + pt.y.toFixed(1);
}
parts.push('<path d="' + d + '" fill="none" stroke="#adb5bd" stroke-width="0.6" opacity="0.12"/>');
}
// Sphere outline
parts.push('<circle cx="0" cy="0" r="' + R + '" fill="none" stroke="#dee2e6" stroke-width="1.2"/>');
// Decision boundary great circle (estimated)
var dPath = '';
var decPts = res.decPoints;
for (var j = 0; j < decPts.length; j++) {
var pt = proj(decPts[j]);
dPath += (j===0 ? 'M' : 'L') + ' ' + pt.x.toFixed(1) + ' ' + pt.y.toFixed(1);
}
parts.push('<path d="' + dPath + '" fill="none" stroke="#e03131" stroke-width="2" stroke-dasharray="10,5" opacity="0.85"/>');
// True decision boundary
var tdPath = '';
for (var j = 0; j < res.trueDecPoints.length; j++) {
var pt = proj(res.trueDecPoints[j]);
tdPath += (j===0 ? 'M' : 'L') + ' ' + pt.x.toFixed(1) + ' ' + pt.y.toFixed(1);
}
parts.push('<path d="' + tdPath + '" fill="none" stroke="#2b8a3e" stroke-width="1.5" stroke-dasharray="5,4" opacity="0.6"/>');
// Data points (depth sorted)
var dataPts = [];
for (var i = 0; i < res.n; i++) {
var pt = proj(res.X[i]);
dataPts.push({x: pt.x, y: pt.y, z: pt.z, cls: res.Y[i], idx: i});
}
dataPts.sort(function(a, b) { return a.z - b.z; }); // far first
for (var i = 0; i < dataPts.length; i++) {
var dp = dataPts[i];
var col = dp.cls === 1 ? '#e03131' : '#1971c2';
var opacity = 0.35 + 0.45 * ((dp.z + 1) / 2); // brighter in front
if (dp.cls === 1) {
parts.push('<circle cx="' + dp.x.toFixed(1) + '" cy="' + dp.y.toFixed(1) +
'" r="3.5" fill="' + col + '" opacity="' + opacity.toFixed(2) + '"/>');
} else {
var s = 5;
var px = dp.x, py = dp.y;
parts.push('<polygon points="' + px.toFixed(1) + ',' + (py-s).toFixed(1) + ' ' +
(px-s*0.866).toFixed(1) + ',' + (py+s*0.5).toFixed(1) + ' ' +
(px+s*0.866).toFixed(1) + ',' + (py+s*0.5).toFixed(1) +
'" fill="' + col + '" opacity="' + opacity.toFixed(2) + '"/>');
}
}
// β̂ direction (arrow from origin through muHat pointing toward betaHat)
var muP = proj(res.muHat);
var betaP = proj(res.betaHat);
// Draw line from muHat toward betaHat direction
// The discriminating direction at muHat: tangent toward betaHat
var tb = [res.betaHat[0]-res.muHat[0]*vdot_local(res.muHat, res.betaHat),
res.betaHat[1]-res.muHat[1]*vdot_local(res.muHat, res.betaHat),
res.betaHat[2]-res.muHat[2]*vdot_local(res.muHat, res.betaHat)];
var ntb = Math.sqrt(tb[0]*tb[0]+tb[1]*tb[1]+tb[2]*tb[2]);
if (ntb > 1e-10) {
tb = [tb[0]/ntb, tb[1]/ntb, tb[2]/ntb];
// Extend along tangent: muHat + t * tb (projected back to sphere)
var arrowLen = 0.55;
var tip = [res.muHat[0] + arrowLen*tb[0], res.muHat[1] + arrowLen*tb[1], res.muHat[2] + arrowLen*tb[2]];
var nTip = Math.sqrt(tip[0]*tip[0]+tip[1]*tip[1]+tip[2]*tip[2]);
tip = [tip[0]/nTip, tip[1]/nTip, tip[2]/nTip];
var tipP = proj(tip);
// Draw an arc from muHat to tip
var arcPts = [];
var nArc = 20;
for (var a = 0; a <= nArc; a++) {
var t = a * arrowLen / nArc;
var interp = [res.muHat[0] + t*tb[0], res.muHat[1] + t*tb[1], res.muHat[2] + t*tb[2]];
var ni = Math.sqrt(interp[0]*interp[0]+interp[1]*interp[1]+interp[2]*interp[2]);
interp = [interp[0]/ni, interp[1]/ni, interp[2]/ni];
var ip = proj(interp);
arcPts.push(ip);
}
var arcPath = '';
for (var a = 0; a < arcPts.length; a++) {
arcPath += (a===0 ? 'M' : 'L') + ' ' + arcPts[a].x.toFixed(1) + ' ' + arcPts[a].y.toFixed(1);
}
parts.push('<path d="' + arcPath + '" fill="none" stroke="#1971c2" stroke-width="3" stroke-linecap="round"/>');
// Arrowhead
var ah = proj(tip);
parts.push('<circle cx="' + ah.x.toFixed(1) + '" cy="' + ah.y.toFixed(1) + '" r="5.5" fill="#1971c2"/>');
}
// μ̂ marker (anchor)
parts.push('<circle cx="' + muP.x.toFixed(1) + '" cy="' + muP.y.toFixed(1) +
'" r="6" fill="#2b8a3e" stroke="#fff" stroke-width="1.5"/>');
// β̂ marker
parts.push('<circle cx="' + betaP.x.toFixed(1) + '" cy="' + betaP.y.toFixed(1) +
'" r="4.5" fill="#1971c2" stroke="#fff" stroke-width="1"/>');
// True β* marker
var betaTrueP = proj(res.betaStarDir);
parts.push('<circle cx="' + betaTrueP.x.toFixed(1) + '" cy="' + betaTrueP.y.toFixed(1) +
'" r="4" fill="#2b8a3e" stroke="#fff" stroke-width="1" opacity="0.7"/>');
return parts.join('\n');
})()}
</g>
<!-- Legend -->
<g transform="translate(15, 400)">
<circle cx="0" cy="0" r="4" fill="#e03131" opacity="0.7"/><text x="8" y="4" font-size="10" fill="#495057">Y=1</text>
<polygon points="70,-3.5 67,3 73,3" fill="#1971c2" opacity="0.7"/><text x="78" y="4" font-size="10" fill="#495057">Y=0</text>
<circle cx="145" cy="0" r="4.5" fill="#2b8a3e" stroke="#fff" stroke-width="1"/><text x="153" y="4" font-size="10" fill="#495057">μ̂</text>
<circle cx="175" cy="0" r="4" fill="#1971c2" stroke="#fff" stroke-width="1"/><text x="183" y="4" font-size="10" fill="#495057">β̂</text>
<line x1="208" y1="0" x2="222" y2="0" stroke="#e03131" stroke-width="1.5" stroke-dasharray="6,3"/>
<text x="226" y="4" font-size="10" fill="#495057">Decision boundary (est.)</text>
</g>
<text x="220" y="426" text-anchor="middle" font-size="11" fill="#6c757d">Drag the sphere to rotate and tilt</text>
</svg>
</div>
<!-- Right panel: stats and info -->
<div style="flex: 1; min-width: 280px;">
<h4>Model Summary</h4>
<div style="padding: 12px; background: #f8f9fa; border-radius: 6px;">
<table style="width: 100%; border-collapse: collapse; font-size: 0.9em;">
<tr><td style="padding: 3px 8px;">True β*</td>
<td style="padding: 3px 8px; text-align: right; font-family: monospace;">θ=${s2result.betaStarThetaDeg.toFixed(0)}°, φ=0°</td></tr>
<tr><td style="padding: 3px 8px;">Estimated β̂</td>
<td style="padding: 3px 8px; text-align: right; font-weight: bold; color: #1971c2; font-family: monospace;">
θ=${(() => { var v = s2result.betaHat; return (Math.acos(Math.max(-1,Math.min(1,v[2])))*180/Math.PI).toFixed(0); })()}°,
φ=${(() => { var v = s2result.betaHat; return (Math.atan2(v[1],v[0])*180/Math.PI).toFixed(0); })()}°</td></tr>
<tr><td style="padding: 3px 8px;">‖μ̂ − μ*‖ (geodesic)</td>
<td style="padding: 3px 8px; text-align: right; font-family: monospace;">${(() => { var v = s2result.muHat; return (Math.acos(Math.max(-1,Math.min(1,v[2])))*180/Math.PI).toFixed(2); })()}°</td></tr>
<tr><td style="padding: 3px 8px;">‖β̂ − β*‖ (geodesic)</td>
<td style="padding: 3px 8px; text-align: right; font-family: monospace;">${(() => { function gd(a,b){return Math.acos(Math.max(-1,Math.min(1,a[0]*b[0]+a[1]*b[1]+a[2]*b[2])));} return (gd(s2result.betaHat,s2result.betaStarDir)*180/Math.PI).toFixed(2); })()}°</td></tr>
<tr><td style="padding: 3px 8px;">Classification accuracy</td>
<td style="padding: 3px 8px; text-align: right; font-weight: bold;">${(s2result.accuracy * 100).toFixed(1)}%</td></tr>
<tr><td style="padding: 3px 8px;">Max log-likelihood</td>
<td style="padding: 3px 8px; text-align: right; font-family: monospace;">${s2result.bestLL.toFixed(4)}</td></tr>
<tr><td style="padding: 3px 8px;">Sample size</td>
<td style="padding: 3px 8px; text-align: right;">${s2result.n}</td></tr>
<tr><td style="padding: 3px 8px;">Concentration κ</td>
<td style="padding: 3px 8px; text-align: right;">${s2result.kappa}</td></tr>
<tr><td style="padding: 3px 8px;">Class separation</td>
<td style="padding: 3px 8px; text-align: right;">${s2result.separation.toFixed(0)}°</td></tr>
</table>
</div>
</div>
</div>
`
return attachS2DragRotation(view, ".s2-logistic-globe", "View rotation", "View tilt", -60, 60, s2result, "logistic");
})()
// Helper: vdot for use in the svg generation scope
function vdot_s2(a, b) { return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; }
```
::: {.callout-note title="Geometric interpretation"}
- The **green dot** is the estimated anchor $\hat{\mu}$ (Fréchet mean of covariates).
- The **blue arrow** shows the geodesic from $\hat{\mu}$ toward $\hat{\beta}$ — this is the **discriminating direction**.
- The **dashed red great circle** is the estimated decision boundary — points whose geodesic from $\hat{\mu}$ is orthogonal to the $\hat{\beta}$ direction.
- Points on the blue-arrow side of the boundary get higher predicted probability of $Y = 1$.
- Drag the sphere to rotate and tilt the view directly; the sliders remain available for preset camera angles.
:::
::: {.callout-note title="How to read the plot"}
The anchor $\hat{\mu}$ is the Fréchet mean of the covariates (green dot). The geodesic from $\hat{\mu}$ toward $\hat{\beta}$ (blue arc with arrow) is the discriminating direction. The dashed red great circle is the estimated decision boundary — points with positive Alexandrov inner product (on the $\hat{\beta}$ side) are classified as $Y = 1$; points on the opposite side as $Y = 0$. Depth-cued rendering (lighter = closer to viewer) helps perceive the 3D structure.
:::
::: {.callout-tip title="Try these experiments"}
- **Rotate the view:** Drag the sphere left/right to rotate and up/down to tilt. The 3D structure — especially the great-circle decision boundary — is much clearer when viewed from multiple perspectives.
- **Adjust class separation:** Increase it to make the red and blue classes more distinct across the true boundary; decrease it to create greater overlap and a harder classification problem.
- **Increase $\kappa$:** Each class forms a tighter cloud around its own center, reducing overlap and making classification easier.
- **Increase $n$:** The estimated $\hat{\beta}$ converges to $\beta^*$ and the decision boundary aligns with the true boundary (green dashed).
- **Compare to S¹:** On the sphere, the decision boundary is a great circle (a 1D curve) rather than two points — this richer geometry allows more flexible classification while remaining fully interpretable.
- **Observe $\hat{\mu}$ estimation:** The estimated anchor (green dot) should be near the north pole (the true $\mu^*$). The small geodesic error $\|\hat{\mu} - \mu^*\|$ confirms that the Fréchet mean is well estimated even with moderate $n$.
:::
## Key Takeaways
- **Metric-space logistic regression** reverses the roles compared to Fréchet regression: the response is binary ($Y \in \{0,1\}$) and the **covariate** $X$ lives in a metric space.
- The model replaces the Euclidean inner product $X^\top\beta^*$ with the **Alexandrov inner product** $d(\mu^*, X) d(\mu^*, \beta^*) \cos\angle_{\mu^*}(X, \beta^*)$ anchored at the Fréchet mean $\mu^*$ of the covariates.
- **Estimation** is two-stage: estimate the anchor via the sample Fréchet mean, then estimate the coefficient $\beta$ via maximum likelihood.
- **LIPO global optimization** handles the non-convex, non-differentiable likelihood landscape using only metric distances and a Hölder continuity bound. The Hölder constant $K$ is critical: too small leads to premature convergence, too large degenerates to random search.
- The fitted geodesic from $\hat{\mu}$ to $\hat{\beta}$ is **geometrically interpretable**: on $\mathbb{S}^2$ the decision boundary is a great circle orthogonal to the discriminating direction; on SPD manifolds it identifies which feature pairs (brain connections, asset correlations) drive the classification.
- The **plug-in classifier** achieves competitive accuracy while remaining geometrically transparent — the decision boundary is orthogonal to $\hat{\beta}$ at $\hat{\mu}$.
- **Convergence rates** and **classification risk bounds** depend on local metric entropy and the margin exponent $\alpha_L$.
## Exercises
1. **Alexandrov inner product in Euclidean space.** Verify that for $\mathcal{M} = \mathbb{R}^D$ with the Euclidean metric and anchor $\mu^* = 0$, the metric-space logistic model reduces exactly to ordinary logistic regression. <a href="javascript:void(0)" onclick="showSolution('l10-sol-1')" class="solution-link">📝 Show Solution</a>
2. **Decision boundary geometry.** Show that on $\mathbb{S}^2$, the decision boundary $\{x : h(\hat{\beta}; x, \hat{\mu}) = 0\}$ is a great circle lying in the plane orthogonal to the tangent direction of $\hat{\beta}$ at $\hat{\mu}$.<a href="javascript:void(0)" onclick="showSolution('l10-sol-4')" class="solution-link">📝 Show Solution</a>
3. **Interpretability vs. black-box.** The metric-space logistic model achieves slightly lower AUC than SPDNet (0.9060 vs 0.9205) but much higher interpretability. Discuss the tradeoff: in what applications would you prefer the interpretable model? <a href="javascript:void(0)" onclick="showSolution('l10-sol-5')" class="solution-link">📝 Show Solution</a>
4. **Implementation.** Implement the metric-space logistic model on a manifold (e.g., hyperbolic space) and evaluate its performance on a synthetic dataset.
<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-original strong { color: #495057; }
.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="l10-sol-1">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 1: Euclidean Recovery</h4>
<button onclick="closeSolution('l10-sol-1')" class="solution-dialog-close">✕ Close</button>
</div>
<div class="solution-original"><strong>Exercise:</strong> Verify the metric-space logistic model reduces to ordinary logistic regression when $\mathcal{M} = \mathbb{R}^D$.</div>
<div class="solution-answer">
<strong>Solution:</strong>
For $\mathcal{M} = \mathbb{R}^D$ with $d(x, y) = \|x - y\|_2$ and anchor $\mu^* = 0$:
$$
d(\mu^*, X) = \|X\|_2, \qquad d(\mu^*, \beta^*) = \|\beta^*\|_2.
$$
The Alexandrov angle at $0$ between $X$ and $\beta^*$ is the Euclidean angle:
$$
\cos\angle_0(X, \beta^*) = \frac{X^\top\beta^*}{\|X\|_2\|\beta^*\|_2}.
$$
Therefore,
$$
h(\beta^*; X, \mu^*) = \|X\|_2 \|\beta^*\|_2 \cdot \frac{X^\top\beta^*}{\|X\|_2\|\beta^*\|_2} = X^\top\beta^*.
$$
The model becomes $\log\frac{P_X}{1-P_X} = X^\top\beta^*$, which is exactly ordinary logistic regression (without intercept). The anchor $\mu^* = 0 = \mathbb{E}[X]$ for centered data, and the coefficient $\beta^*$ is the usual logistic coefficient vector.
</div>
</div>
</dialog>
<dialog id="l10-sol-4">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 4: Decision Boundary on 𝕊²</h4>
<button onclick="closeSolution('l10-sol-4')" class="solution-dialog-close">✕ Close</button>
</div>
<div class="solution-original"><strong>Exercise:</strong> Show that the decision boundary on $\mathbb{S}^2$ is a great circle orthogonal to $\hat{\beta}$.</div>
<div class="solution-answer">
<strong>Solution:</strong>
On $\mathbb{S}^2$, the decision boundary is $\{x \in \mathbb{S}^2 : h(\hat{\beta}; x, \hat{\mu}) = 0\}$.
The Alexandrov inner product is:
$$
h(\hat{\beta}; x, \hat{\mu}) = d(\hat{\mu}, x) d(\hat{\mu}, \hat{\beta}) \cos\angle_{\hat{\mu}}(x, \hat{\beta}).
$$
Assuming $d(\hat{\mu}, \hat{\beta}) > 0$, this is zero when either:
1. $d(\hat{\mu}, x) = 0$ (i.e., $x = \hat{\mu}$, a single point — the anchor), or
2. $\cos\angle_{\hat{\mu}}(x, \hat{\beta}) = 0$ (i.e., $\angle_{\hat{\mu}}(x, \hat{\beta}) = \pi/2$).
**Why it's a great circle:** The condition $\angle_{\hat{\mu}}(x, \hat{\beta}) = \pi/2$ means the tangent vector at $\hat{\mu}$ toward $x$ is orthogonal to the tangent vector at $\hat{\mu}$ toward $\hat{\beta}$. Let $v_\beta$ be the unit tangent vector at $\hat{\mu}$ pointing toward $\hat{\beta}$ (i.e., the direction of the geodesic $\gamma_{\hat{\mu} \to \hat{\beta}}$). The set of all unit tangent vectors at $\hat{\mu}$ orthogonal to $v_\beta$ forms a 1-dimensional circle in the tangent space $T_{\hat{\mu}}\mathbb{S}^2 \cong \mathbb{R}^2$. Following each of these directions along a geodesic for distance $\pi/2$ from $\hat{\mu}$ sweeps out exactly a **great circle** on $\mathbb{S}^2$.
Equivalently, the great circle lies in the plane whose normal vector is $v_\beta$ (embedded in $\mathbb{R}^3$). The plane passes through the origin and is orthogonal to the direction of $\hat{\beta}$ at $\hat{\mu}$. The intersection of this plane with the unit sphere is the decision boundary great circle.
</div>
</div>
</dialog>
<dialog id="l10-sol-5">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 5: Interpretability Tradeoff</h4>
<button onclick="closeSolution('l10-sol-5')" class="solution-dialog-close">✕ Close</button>
</div>
<div class="solution-original"><strong>Exercise:</strong> Discuss when you would prefer the interpretable metric-space logistic model over a black-box classifier like SPDNet.</div>
<div class="solution-answer">
<strong>Solution:</strong>
**Prefer the interpretable model when:**
1. **Scientific discovery is the goal.** The fMRI application is a perfect example — the goal is not just to classify tasks but to understand *which brain regions and connections* drive the difference. The fitted geodesic directly identifies these.
2. **Regulatory or ethical requirements.** In medical diagnosis, credit scoring, or criminal justice, black-box models face regulatory scrutiny (e.g., GDPR "right to explanation"). An interpretable model can justify its decisions.
3. **Small sample sizes.** The metric-space logistic model has only one parameter object $\beta$ (same "dimension" as a single covariate). SPDNet and other deep learning methods have thousands of parameters and require large $n$.
4. **Generalization across domains.** A geometrically motivated model may generalize better when the data-generating process changes slightly, because it encodes structural knowledge (the metric geometry) rather than learning everything from data.
5. **Debugging and iteration.** When the model makes mistakes, you can inspect *why*: which geodesic direction was involved? This enables iterative model improvement.
**Prefer a black-box model when:**
1. **Pure prediction** is the sole objective (e.g., winning a Kaggle competition).
2. **The decision boundary is highly nonlinear** and cannot be well-approximated by a single geodesic half-space.
3. **Very large $n$** is available, making deep learning feasible and likely to outperform.
**A better approach** is often to fit both: use the interpretable model for scientific insight and the black-box model as a performance benchmark.
</div>
</div>
</dialog>
<script>
function showSolution(id) {
const dialog = document.getElementById(id);
if (dialog) {
dialog.showModal();
if (window.MathJax && MathJax.typesetPromise) {
MathJax.typesetPromise([dialog]).catch(function(err) { console.log('MathJax typeset error:', err); });
}
}
}
function closeSolution(id) {
const dialog = document.getElementById(id);
if (dialog) { dialog.close(); }
}
document.addEventListener('click', function(e) {
if (e.target.tagName === 'DIALOG') { e.target.close(); }
});
</script>
## Self-Assessment Quiz
Test your understanding of this lecture with the interactive MCQ quiz:
👉 **[Lecture 10 Quiz — 10 Multiple-Choice Questions](../quizzes/lecture-10-quiz.qmd)**
## Further Reading
- @LinLin2025Binary — The foundational paper on metric-space logistic regression and classification.
- @LinSPDCholesky2019 — Log-Cholesky metric on SPD matrices, used in the fMRI application.
- @ArsignyEtAl2007 — Log-Euclidean metrics for SPD matrices.
- @Moakher2005GeometricMean — Affine-invariant geometric mean on SPD matrices.
- @Sturm2003 — Hadamard space theory: convexity and uniqueness results.
