---
title: "Lecture 19: Riemannian Manifolds — Local Polynomial Regression"
subtitle: "Intrinsic local polynomial fitting for manifold-valued responses"
format:
html:
code-fold: true
code-tools: true
code-copy: true
pdf:
documentclass: scrartcl
pdf-engine: xelatex
toc: true
number-sections: true
geometry:
- margin=1in
colorlinks: true
bibliography: source/ref.bib
---
## Learning Goals
By the end of this lecture, learners should be able to:
- Contrast local Fréchet regression (metric-only) with intrinsic local polynomial regression (manifold structure).
- Write the intrinsic mean-zero error model $Y = \exp_{\mu(X)}(\varepsilon)$.
- Formulate the order-$\ell$ local polynomial curve expansion in a tangent space.
- Define the population and sample local polynomial criteria.
- Explain how $\ell = 0$ recovers local-constant (kernel Fréchet) and $\ell = 1$ gives local-linear regression.
- Discuss computational strategies and asymptotic theory.
- Apply intrinsic local polynomial regression to SPD covariance data and interpret tangent-space slopes as rates of change of risk structure.
- Compare local-constant (metric-space kernel Fréchet) and local-linear (manifold intrinsic) estimators on simulated portfolio stress-testing data, including boundary behavior and MSE.
## Motivation: Metric vs. Manifold Methods
Local Fréchet regression (Lectures 7–8) works in any metric space — it only needs distances. When the response space is a **smooth Riemannian manifold**, we can do more: use the exponential and logarithm maps to linearize locally, fit a polynomial in a **single tangent space**, and map back.
This is the approach of @YuanZhuLinMarron2012, who developed the methodology for $\mathcal{S}_{++}^m$ under the affine-invariant and log-Euclidean metrics.
## Euclidean Local Polynomial Regression
Before introducing manifold-valued responses, it is useful to recall the Euclidean construction that local polynomial regression generalizes. Suppose that
$$
Y = m(X) + \varepsilon, \qquad \mathbb{E}(\varepsilon \mid X)=0,
$$
where $X\in\mathbb{R}$ and $Y\in\mathbb{R}^q$; the scalar-response case is obtained by taking $q=1$. Our target is the conditional mean curve $m(x)=\mathbb{E}(Y\mid X=x)$ and, when needed, its derivatives $m^{(j)}(x)$.
### Why Fit a Polynomial Locally?
A single global polynomial imposes one shape over the entire predictor range and can behave poorly when the true regression curve has spatially varying curvature. A local fit instead uses only observations with $X_i$ close to a target point $x$. It therefore combines three ideas:
1. **Taylor approximation:** every sufficiently smooth curve is approximately polynomial in a small neighborhood;
2. **Kernel localization:** observations nearer to $x$ receive more weight;
3. **Repeated pointwise fitting:** moving $x$ across the predictor range produces a smooth estimated curve without specifying a global parametric model.
Local polynomial regression contains several familiar smoothers. Order $\ell=0$ gives the local-constant or Nadaraya--Watson estimator, $\ell=1$ gives local-linear regression, and $\ell=2$ gives local-quadratic regression. The local-linear fit is especially important: it reproduces linear trends exactly and removes the leading boundary bias of a local-constant smoother under standard smoothness and design conditions.
### From Taylor Expansion to a Local Criterion
Fix a target $x$ and write $u=z-x$. If $m$ is $\ell+1$ times differentiable near $x$, then
$$
m(x+u)
= \sum_{j=0}^{\ell}\frac{u^j}{j!}m^{(j)}(x)
+ R_{\ell+1}(u),
\qquad R_{\ell+1}(u)=O(|u|^{\ell+1}).
$$
This suggests the local polynomial model
$$
p_{\beta,x}(z)
= \sum_{j=0}^{\ell}\frac{(z-x)^j}{j!}\beta_j,
\qquad \beta_j\in\mathbb{R}^q,
$$
where $\beta_0$ approximates $m(x)$ and $\beta_j$ approximates $m^{(j)}(x)$. Centering the powers at $x$ is crucial: the fitted value at the target is simply the intercept,
$$
p_{\beta,x}(x)=\beta_0.
$$
Let $K$ be a nonnegative kernel, usually symmetric with $\int K(u)\,du=1$, and define $K_h(u)=K(u/h)/h$ for bandwidth $h>0$. The **population local-polynomial target** is any minimizer
$$
\beta_{\ell,h}(x)
= \arg\min_{\beta_0,\ldots,\beta_\ell}
\mathbb{E}\!\left[
K_h(X-x)
\left\|Y-\sum_{j=0}^{\ell}\frac{(X-x)^j}{j!}\beta_j\right\|^2
\right].
$$
The kernel restricts attention to a neighborhood of width roughly $h$. As $h\to0$, the Taylor approximation becomes more accurate, but fewer observations contribute appreciably and the variance increases.
::: {.callout-note title="Population target versus regression function"}
For fixed $h$, $\beta_{0,\ell,h}(x)$ is a bandwidth-dependent smoothed target and need not equal $m(x)$ exactly. Under the usual smoothness, design-density, and kernel assumptions, it approaches $m(x)$ as $h\to0$. Consistency also requires the effective local sample size to grow, typically $nh\to\infty$ for a scalar predictor.
:::
### Sample Estimator as Weighted Least Squares
Given independent observations $(X_i,Y_i)$, $i=1,\ldots,n$, define
$$
u_i=X_i-x,
\qquad
r_\ell(u_i)=\left(1,u_i,\frac{u_i^2}{2!},\ldots,\frac{u_i^\ell}{\ell!}\right)^\top,
\qquad
w_i(x)=K_h(u_i).
$$
Stack the coefficient vectors as rows of $B=(\beta_0,\ldots,\beta_\ell)^\top\in\mathbb{R}^{(\ell+1)\times q}$. The estimator minimizes
$$
\widehat B_{\ell,h}(x)
=\arg\min_{B\in\mathbb{R}^{(\ell+1)\times q}}
\sum_{i=1}^n w_i(x)
\left\|Y_i-B^\top r_\ell(u_i)\right\|^2,
$$
and
$$
\widehat m(x)=\widehat B_{\ell,h}(x)^\top e_0.
$$
Here $e_j$ selects row $j+1$.
### Estimation Choices and Statistical Tradeoffs
**Polynomial order.** Local linear ($\ell=1$) is a strong default for estimating $m(x)$: it corrects boundary bias without the instability of a high-order fit. Local quadratic or cubic fits are useful when curvature or higher derivatives are scientifically meaningful, but they require more data within each neighborhood and can have larger variance.
**Bandwidth.** The bandwidth usually matters more than the polynomial order. For local linear regression at an interior point, the leading mean-estimation bias is typically $O(h^2)$ and the variance is $O\{(nh)^{-1}\}$ for a scalar predictor. Balancing squared bias and variance gives the familiar scale $h\asymp n^{-1/5}$, although constants depend on the kernel, design density, noise variance, and curvature of $m$. For the $j$th derivative, variance is of order $O\{(nh^{2j+1})^{-1}\}$, so derivative estimation is more variable and requires a bandwidth criterion tailored to the derivative target.
**Data-driven bandwidth selection.** A practical prediction-oriented choice is leave-one-out cross-validation,
$$
\widehat h_{\mathrm{CV}}
=\arg\min_{h\in\mathcal H}
\frac{1}{n}\sum_{i=1}^n
\left\|Y_i-\widehat m_{-i,h}(X_i)\right\|^2,
$$
where $\widehat m_{-i,h}$ is fitted without observation $i$. Plug-in selectors are another option when reliable estimates of curvature and noise are available. If the inferential target is a derivative rather than the mean curve, the bandwidth should be selected for that derivative-specific loss.
**Local identifiability.** The matrix $R_x^\top W_xR_x$ must have full rank. Failure occurs when too few distinct predictor values receive positive weight, when $h$ is extremely small, or when the local polynomial order is too high for the available design. In computation, inspecting its condition number is preferable to blindly applying a matrix inverse.
::: {.callout-tip title="Bridge to manifold-valued responses"}
The Euclidean estimator uses two operations that are not intrinsically available on a curved manifold: subtracting $Y_i-p_{\beta,x}(X_i)$ and adding polynomial coefficient vectors based at different points. The manifold construction below replaces the Euclidean polynomial by one in a **single tangent space**, maps it to the manifold with $\exp$, and replaces squared Euclidean residuals by squared geodesic distances. The local intercept, slope, bandwidth, and polynomial-order interpretations remain the same.
:::
## Local Polynomial Regression on Manifolds
**Intrinsic Error Model.**
For scalar predictor $X \in \mathbb{R}$ and manifold response $Y \in \mathcal{M}$, assume the conditional Fréchet mean $\mu(x) = \arg\min_\omega \mathbb{E}\{d^2(Y, \omega) \mid X = x\}$ is unique. If the conditional support lies in a normal neighborhood of $\mu(x)$ and differentiation under the conditional expectation is valid, its first-order condition is
$$
\mathbb{E}\{\log_{\mu(x)}(Y) \mid X = x\} = 0.
$$
This is the manifold analogue of $\mathbb{E}\{Y - \mu(X) \mid X\} = 0$. A convenient error representation is
$$
Y = \exp_{\mu(X)}(\varepsilon), \qquad \varepsilon \in T_{\mu(X)}\mathcal{M}, \qquad \mathbb{E}(\varepsilon \mid X) = 0.
$$
**Tangent-Space Taylor Expansion.**
Fix a target point $x$. Near $x$, the mean curve $\mu$ can be expanded in $T_{\mu(x)}\mathcal{M}$:
$$
\log_{\mu(x)}(\mu(x + u)) = u\,\dot{\mu}(x) + \frac{u^2}{2}\,\ddot{\mu}(x) + o(u^2), \qquad u \to 0,
$$
where $\dot{\mu}(x)$ is the velocity and $\ddot{\mu}(x)=\nabla_{\dot\mu}\dot\mu(x)$ is the covariant acceleration, both represented in $T_{\mu(x)}\mathcal{M}$. The expansion requires that $\mu(x+u)$ remain inside the normal neighborhood used by the logarithm map.
**Order-$\ell$ Local Polynomial Curves.**
For parameters $\theta = (\omega, v_1, \ldots, v_\ell)$ with $\omega \in \mathcal{M}$, $v_j \in T_\omega\mathcal{M}$, define:
$$
m_{\theta, x}(z) = \exp_\omega\!\left\{\sum_{j=1}^\ell \frac{(z - x)^j}{j!}\,v_j\right\}, \qquad z \in \mathbb{R}.
$$
- $\omega = m_{\theta,x}(x)$ is the **local intercept** (the fitted value at $x$)
- $v_1$ is the local **slope** (velocity of the mean curve at $x$)
- $v_2$ is the acceleration.
**Estimation.**
For kernel $K$ and bandwidth $h > 0$, the population criterion is
$$
Q_{\ell, h}(\theta; x) = \mathbb{E}\Bigl[K_h(X - x)\, d^2(Y, m_{\theta, x}(X))\Bigr].
$$
A minimizer $\theta_{\ell, h}(x) = (\mu_{\ell, h}(x), \beta_{1,\ell,h}(x), \ldots)$ gives the **population target**; the fitted regression value at $x$ is the intercept $\mu_{\ell, h}(x)$.
The sample estimator is
$$
\hat{Q}_{\ell, h}(\theta; x) = \frac{1}{n}\sum_{i=1}^n K_h(X_i - x)\, d^2(Y_i, m_{\theta, x}(X_i)),
$$
with $\hat{\theta}_{\ell, h}(x) = (\hat{\mu}_{\ell, h}(x), \hat{\beta}_{1,\ell,h}(x), \ldots)$ being any minimizer.
**Special Cases.**
- **$\ell = 0$:** $m_{\theta, x}(z) \equiv \omega$, so $\hat{\mu}_{0,h}(x) = \arg\min_\omega \frac{1}{n}\sum K_h(X_i - x)d^2(Y_i, \omega)$ — this is exactly the **kernel Fréchet regression** estimator (Lecture 7).
**Contrast with Local Fréchet Regression.**
| Aspect | Local Fréchet Regression | Intrinsic Local Polynomial |
|---|---|---|
| **Structure needed** | Only metric | Smooth manifold + exp/log |
| **Linearization** | Predictor-side weights $s(z,x,h)$ | Tangent space at $\omega$ |
| **Weight interpretation** | Absorbed into weight function | Direct tangent-space polynomial |
| **Computation** | Weighted Fréchet mean | Optimization over $(\omega, v_1, \ldots, v_\ell)$ |
| **Generality** | Any metric space | Riemannian manifolds only |
| **Slope interpretation** | Implicit | Explicit: $v_1 \in T_\omega\mathcal{M}$ |
## Computation
### Optimization Landscape on the SPD Manifold
For responses $Y_i \in \mathcal{S}_{++}^m$ and the affine-invariant metric $d_{\mathrm{AI}}(A, B) = \|\log(A^{-1/2} B A^{-1/2})\|_F$, the objective is
$$
\hat{Q}_{\ell, h}(\theta; x)
= \frac{1}{n}\sum_{i=1}^n K_h(X_i - x)\,
d_{\mathrm{AI}}^2\!\left(Y_i,\;
\exp_\omega\!\left\{\sum_{j=1}^\ell \frac{(X_i - x)^j}{j!}\,v_j\right\}\right),
$$
where $\omega \in \mathcal{S}_{++}^m$ and $v_j \in \operatorname{Sym}(m)$. The parameters live in a space of dimension
$$
\dim(\theta) = \frac{m(m+1)}{2} \times (\ell + 1).
$$
For $m = 6$ (e.g., six assets) and $\ell = 1$, this is $21 \times 2 = 42$ parameters; for $\ell = 2$, it's 63 parameters. The objective is non-convex because the exponential map at $\omega$ depends nonlinearly on $\omega$, and the distance function involves matrix logarithms, inverses, and square roots.
For the trace/affine-invariant formulation studied by @YuanZhuLinMarron2012, the criterion for $\ell>0$ can have multiple local minima, so results from a local optimizer may depend on initialization.
### AESAMC: Annealing Evolutionary Stochastic Approximation Monte Carlo
@YuanZhuLinMarron2012 use the **annealing evolutionary stochastic approximation Monte Carlo** (AESAMC) algorithm of @Liang2011AESAMC for their non-convex trace-metric fit. AESAMC combines:
| Component | Mechanism | Role |
|---|---|---|
| **Stochastic approximation Monte Carlo** | Adaptive weights over energy subregions | Encourages exploration across objective levels |
| **Annealing** | Progressively restricts attention to lower-energy regions | Concentrates the search near the best objective values found |
| **Evolutionary Monte Carlo (EMC)** | Population-based crossover and mutation operators | Explores the parameter space efficiently through genetic operations on a population of candidate solutions |
At a high level, the method partitions objective values into energy regions, evolves a population with mutation/crossover proposals, adapts regional weights, and lowers the active energy ceiling as better solutions appear. The best point visited is returned.
### Computation on Log-Euclidean Metric
When the log-Euclidean (LE) metric is used instead of the affine-invariant metric, the SPD manifold becomes globally isometric to $\operatorname{Sym}(m)$ through the matrix logarithm [@ArsignyEtAl2007]. Write $a_j=D\log_\omega[v_j]$ for the log-coordinate representation of the tangent vector $v_j$. Then:
- The Riemannian exponential is $\operatorname{Exp}^{\mathrm{LE}}_\omega(v)=\exp\{\log\omega+D\log_\omega[v]\}$; equivalently, a log-coordinate increment $a$ maps back as $\exp(\log\omega+a)$.
- The distance is $d_{\mathrm{LE}}(A, B) = \|\log A - \log B\|_F$.
- The objective simplifies to **weighted least squares** in the log-domain:
$$
\hat{Q}_{\ell, h}(\theta; x)
= \frac{1}{n}\sum_{i=1}^n K_h(X_i - x)\,
\left\|\log Y_i - \log\omega - \sum_{j=1}^\ell \frac{(X_i - x)^j}{j!}\,a_j\right\|_F^2.
$$
This is **convex** in $(\log\omega,a_1,\ldots,a_\ell)$ and has the standard weighted-least-squares solution. The interactive SPD demo below uses these log coordinates.
## Interactive Exploration: Local-Constant vs. Tangent Local-Linear Approximation on $S^2$
```{ojs}
//| echo: false
//| output: false
s2_n_control = Inputs.range([20, 200], {step: 10, value: 80, label: "Sample size n"})
s2_sigma_control = Inputs.range([0.03, 0.40], {step: 0.01, value: 0.12, label: "Noise σ"})
s2_h_control = Inputs.range([0.06, 0.50], {step: 0.02, value: 0.18, label: "Bandwidth h"})
s2_curve_control = Inputs.select([
"Spiral (pole to equator)",
"Great circle",
"Sinusoidal oscillation"
], {value: "Spiral (pole to equator)", label: "True curve"})
s2_seed_control = Inputs.range([1, 20], {step: 1, value: 5, label: "Random seed"})
s2_n = Generators.input(s2_n_control)
s2_sigma = Generators.input(s2_sigma_control)
s2_h = Generators.input(s2_h_control)
s2_curve_label = Generators.input(s2_curve_control)
s2_curve = s2_curve_label === "Great circle"
? "great_circle"
: s2_curve_label === "Sinusoidal oscillation"
? "sinusoidal"
: "spiral"
s2_seed = Generators.input(s2_seed_control)
s2_controls_view = html`
<style>
.s2-slider-grid { display:flex; flex-wrap:wrap; gap:8px 20px; width:100%; max-width:840px; margin:0 auto 12px; font:0.85em system-ui,sans-serif; container-type:inline-size; }
.s2-slider-grid > * { min-width:0; margin:0; }
.s2-curve-control { flex:1 1 100%; }
.s2-numeric-control { flex:1 1 calc((100% - 20px)/2); }
.s2-slider-grid input[type="number"] { width:10rem !important; }
.s2-slider-grid select { width:100%; max-width:100%; }
@container (max-width:480px) { .s2-numeric-control { flex-basis:100%; } }
</style>
<div class="s2-slider-grid">
<div class="s2-curve-control">${s2_curve_control}</div>
<div class="s2-numeric-control">${s2_n_control}</div>
<div class="s2-numeric-control">${s2_sigma_control}</div>
<div class="s2-numeric-control">${s2_h_control}</div>
<div class="s2-numeric-control">${s2_seed_control}</div>
</div>`
// ===== S² Riemannian Geometry Utilities =====
function norm3(v) { return Math.sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]); }
function dot3(a, b) { return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; }
function cross3(a, b) { return [a[1]*b[2]-a[2]*b[1], a[2]*b[0]-a[0]*b[2], a[0]*b[1]-a[1]*b[0]]; }
function scale3(v, s) { return [v[0]*s, v[1]*s, v[2]*s]; }
function add3(a, b) { return [a[0]+b[0], a[1]+b[1], a[2]+b[2]]; }
function sub3(a, b) { return [a[0]-b[0], a[1]-b[1], a[2]-b[2]]; }
function normalize3(v) { const n = norm3(v); return n < 1e-12 ? [0,0,1] : scale3(v, 1/n); }
// Exponential map on S²: exp_p(v)
function s2Exp(p, v) {
const nv = norm3(v);
if (nv < 1e-12) return [...p];
const vDir = scale3(v, 1/nv);
return add3(scale3(p, Math.cos(nv)), scale3(vDir, Math.sin(nv)));
}
// Logarithm map on S²: log_p(q)
function s2Log(p, q) {
const dpq = Math.max(-1, Math.min(1, dot3(p, q)));
const d = Math.acos(dpq);
if (d < 1e-12) return [0, 0, 0];
const tangent = sub3(q, scale3(p, dpq));
const tangentNorm = norm3(tangent);
if (tangentNorm < 1e-8) {
if (dpq > 0) return [0, 0, 0];
// At or numerically too close to the antipode, log_p(q) is not unique.
const cand = Math.abs(p[0]) < 0.9 ? [1, 0, 0] : [0, 1, 0];
const proj = sub3(cand, scale3(p, dot3(cand, p)));
return scale3(normalize3(proj), Math.PI);
}
return scale3(tangent, d/tangentNorm);
}
// Orthonormal basis for T_p S²
function tangentBasis(p) {
const eps = 1e-8;
if (p[0]*p[0] + p[1]*p[1] < eps) {
// Near poles
const sign = p[2] > 0 ? 1 : -1;
return [[1, 0, 0], [0, sign, 0]];
}
const r = Math.sqrt(p[0]*p[0] + p[1]*p[1]);
const e1 = [-p[1]/r, p[0]/r, 0];
const e2 = cross3(p, e1);
return [e1, e2];
}
// Distance on S²
function s2Dist(p, q) { return Math.acos(Math.max(-1, Math.min(1, dot3(p, q)))); }
// ===== Rotation and Projection =====
function rotatePoint(p, yaw, pitch) {
// Rotate around Y axis (yaw)
const cy = Math.cos(yaw), sy = Math.sin(yaw);
const x1 = cy*p[0] + sy*p[2];
const y1 = p[1];
const z1 = -sy*p[0] + cy*p[2];
// Rotate around X axis (pitch)
const cp = Math.cos(pitch), sp = Math.sin(pitch);
return [x1, cp*y1 - sp*z1, sp*y1 + cp*z1];
}
function project(p, cx, cy, r) {
return {x: cx + r * p[0], y: cy - r * p[1], z: p[2]};
}
// ===== True Curves on S² =====
function trueCurveS2(t, type) {
// t in [0,1], returns point on S²
if (type === "great_circle") {
// Great circle tilted 25° from equator
const tilt = 0.436; // ~25°
const lon = 2*Math.PI*t;
const ct = Math.cos(tilt), st = Math.sin(tilt);
return [
Math.cos(lon)*ct,
Math.sin(lon),
-Math.cos(lon)*st
];
} else if (type === "sinusoidal") {
// Oscillating latitude
const lon = 2*Math.PI*t;
const lat = 0.9*Math.sin(5*t);
const cl = Math.cos(lat), sl = Math.sin(lat);
return [cl*Math.cos(lon), cl*Math.sin(lon), sl];
} else {
// Spiral: from pole toward equator
const colat = (Math.PI/2)*Math.pow(t, 1.6);
const lon = 4.5*Math.PI*t;
const sc = Math.sin(colat), cc = Math.cos(colat);
return [sc*Math.cos(lon), sc*Math.sin(lon), cc];
}
}
// ===== Data Generation =====
function generateS2Data(n, sigma, curveType, seed) {
const rng = (function(a) {
return function() {
a |= 0; a = a + 0x6D2B79F5 | 0;
var t = Math.imul(a ^ a >>> 15, 1 | a);
t = t + Math.imul(t ^ t >>> 7, 61 | t) ^ t;
return ((t ^ t >>> 14) >>> 0) / 4294967296;
};
})(seed || 42);
function rnorm() {
const u1 = rng(), u2 = rng();
return Math.sqrt(-2*Math.log(Math.max(u1,1e-10)))*Math.cos(2*Math.PI*u2);
}
// Predictor X in [0, 1]
const X = Array.from({length: n}, () => rng());
// Generate Y on S²: true curve + tangent noise
const Y = X.map(xi => {
const mu = trueCurveS2(xi, curveType);
const [e1, e2] = tangentBasis(mu);
// Generate 2D tangent noise scaled by sigma
const a = sigma * rnorm();
const b = sigma * rnorm();
const v = add3(scale3(e1, a), scale3(e2, b));
return s2Exp(mu, v);
});
return {X, Y};
}
// ===== Weighted Fréchet Mean (Local-Constant) on S² =====
function weightedFrechetMeanS2(points, weights) {
const n = points.length;
// Initialize as weighted Euclidean mean, normalized
let mu = [0, 0, 0];
let wSum = 0;
for (let i = 0; i < n; i++) {
mu = add3(mu, scale3(points[i], weights[i]));
wSum += weights[i];
}
if (wSum < 1e-12) return [0, 0, 1];
mu = normalize3(mu);
// Karcher iteration: the weighted average log vector is the descent direction.
function objective(p) {
let value = 0;
for (let i = 0; i < n; i++) {
const d = s2Dist(p, points[i]);
value += weights[i]*d*d;
}
return value;
}
for (let iter = 0; iter < 60; iter++) {
let direction = [0, 0, 0];
let wSum2 = 0;
for (let i = 0; i < n; i++) {
if (weights[i] < 1e-12) continue;
const logVec = s2Log(mu, points[i]);
direction = add3(direction, scale3(logVec, weights[i]));
wSum2 += weights[i];
}
if (wSum2 < 1e-12) break;
direction = scale3(direction, 1/wSum2);
const ng = norm3(direction);
if (ng < 1e-10) break;
const oldValue = objective(mu);
let step = 1;
let candidate = mu;
while (step >= 1/1024) {
candidate = normalize3(s2Exp(mu, scale3(direction, step)));
if (objective(candidate) < oldValue) break;
step *= 0.5;
}
if (step < 1/1024) break;
mu = candidate;
}
return mu;
}
// ===== Estimators =====
// Normalized Gaussian weights computed with a log-sum-exp shift.
// This preserves relative weights even when every unshifted kernel value is tiny.
function normalizedGaussianWeights(X, x, h) {
const logWeights = X.map(xi => -0.5*((xi-x)/h)**2);
const maxLogWeight = Math.max(...logWeights);
const raw = logWeights.map(lw => Math.exp(lw-maxLogWeight));
const total = raw.reduce((a,b) => a+b, 0);
return raw.map(w => w/total);
}
function computeFits(data, h, nGrid, curveType) {
const {X, Y} = data;
const n = X.length;
nGrid = nGrid || 80;
const xGrid = Array.from({length: nGrid + 1}, (_, j) => j/nGrid);
// Local-constant (ℓ=0)
const lcFit = xGrid.map(x => {
const nw = normalizedGaussianWeights(X, x, h);
return weightedFrechetMeanS2(Y, nw);
});
// Two-step tangent local-linear approximation (not the joint intrinsic minimizer)
const llFit = xGrid.map(x => {
const nw0 = normalizedGaussianWeights(X, x, h);
const effectiveN = 1/nw0.reduce((sum, w) => sum+w*w, 0);
// Step 1: local-constant estimate
const wHat = weightedFrechetMeanS2(Y, nw0);
const [e1, e2] = tangentBasis(wHat);
// Step 2: tangent-space WLS for each coordinate
function solveWLS(coordFn) {
let S0 = 0, S1 = 0, S2 = 0, Sy = 0, Sxy = 0;
for (let i = 0; i < n; i++) {
const wi = nw0[i];
const dx = (X[i] - x)/h;
S0 += wi;
S1 += wi*dx;
S2 += wi*dx*dx;
Sy += wi*coordFn(i);
Sxy += wi*dx*coordFn(i);
}
const det = S0*S2 - S1*S1;
if (effectiveN >= 3 && det > 1e-8) {
return {
intercept: (S2*Sy - S1*Sxy)/det,
slope: (S0*Sxy - S1*Sy)/(det*h)
};
}
return {intercept: Sy/Math.max(S0, 1e-15), slope: 0};
}
// Compute log-coordinates of each Y_i in T_{wHat}S²
const logCoords = Y.map(yi => {
const v = s2Log(wHat, yi);
return [dot3(v, e1), dot3(v, e2)];
});
const resA = solveWLS(i => logCoords[i][0]);
const resB = solveWLS(i => logCoords[i][1]);
// Reconstruct: β₀ = α_a·e₁ + α_b·e₂, β₁ = β_a·e₁ + β_b·e₂
const beta0 = add3(scale3(e1, resA.intercept), scale3(e2, resB.intercept));
const beta1 = add3(scale3(e1, resA.slope), scale3(e2, resB.slope));
// Fitted value at x: exp_{wHat}(β₀)
const muHat = s2Exp(wHat, beta0);
// For evaluation at other points z: exp_{wHat}(β₀ + β₁(z-x))
return {wHat, beta0, beta1, muHat};
});
// Evaluate fitted curves at grid points
const lcCurve = lcFit;
const llCurve = xGrid.map((x, idx) => {
// llFit[idx] gives the fitted value at xGrid[idx], which is muHat
return llFit[idx].muHat;
});
// True curve
const trueCurve = xGrid.map(x => trueCurveS2(x, curveType));
// Compute MSEs
let mseLC = 0, mseLL = 0;
for (let k = 0; k <= nGrid; k++) {
const dLC = s2Dist(lcCurve[k], trueCurve[k]);
const dLL = s2Dist(llCurve[k], trueCurve[k]);
mseLC += dLC*dLC;
mseLL += dLL*dLL;
}
mseLC /= (nGrid + 1);
mseLL /= (nGrid + 1);
return {xGrid, lcCurve, llCurve, trueCurve, mseLC, mseLL};
}
// ===== Generate data and fits =====
s2_data = generateS2Data(s2_n, s2_sigma, s2_curve, s2_seed);
s2_fits = computeFits(s2_data, s2_h, 100, s2_curve);
// ===== 3D Canvas Rendering =====
s2_plot = {
const W = 720, H = 500;
const CX = W/2, CY = H/2, R = 200;
const container = html`<div style="font-family:system-ui,sans-serif;max-width:780px;margin:0 auto;">
<canvas id="s2-globe-canvas" width="${W}" height="${H}" role="img" aria-label="Rotatable sphere showing the true curve, observations, and two local regression fits" style="cursor:grab;border:1px solid #dee2e6;border-radius:6px;display:block;width:100%;height:auto;margin:0 auto;background:#f0f4f8;touch-action:none;"></canvas>
<div style="display:flex;justify-content:space-between;align-items:flex-start;margin-top:8px;flex-wrap:wrap;gap:8px;">
<div style="font-size:0.85em;color:#495057;">
<span style="display:inline-block;width:12px;height:12px;border-radius:50%;background:#2b8a3e;margin-right:4px;vertical-align:middle;"></span> True μ(x)
<span style="display:inline-block;width:12px;height:12px;border-radius:50%;background:#f08c00;margin-right:4px;vertical-align:middle;"></span> LC (ℓ=0)
<span style="display:inline-block;width:12px;height:12px;border-radius:50%;background:#1971c2;margin-right:4px;vertical-align:middle;"></span> tangent LL
<span style="color:#adb5bd;">●</span> Data (Yᵢ)</span>
<div style="font-size:0.82em;color:#868e96;">
MSE<sub>LC</sub> = ${s2_fits.mseLC.toFixed(4)} |
MSE<sub>LL</sub> = ${s2_fits.mseLL.toFixed(4)} |
Δ = ${(s2_fits.mseLL - s2_fits.mseLC).toFixed(4)}
</div>
</div>
<p style="font-size:0.78em;color:#868e96;margin-top:4px;text-align:center;">
🖱 Drag to rotate globe. <i>Compare the two fits; the lower-MSE method depends on this sample, bandwidth, and noise level.</i>
</p>
</div>`;
const canvas = container.querySelector('#s2-globe-canvas');
const ctx = canvas.getContext('2d');
// Rotation state
let yaw = 0.55, pitch = -0.35;
let dragging = false, lastMouse = null;
// Precompute meridian and parallel points (on unit sphere)
const N_GRID = 120;
const meridians = [];
for (let k = 0; k < 12; k++) {
const lon = k*Math.PI/6;
const pts = [];
for (let j = 0; j <= N_GRID; j++) {
const phi = -Math.PI/2 + j*Math.PI/N_GRID;
pts.push([Math.cos(phi)*Math.cos(lon), Math.cos(phi)*Math.sin(lon), Math.sin(phi)]);
}
meridians.push(pts);
}
const parallels = [];
for (let k = 1; k < 6; k++) {
const phi = -Math.PI/2 + k*Math.PI/6;
const pts = [];
for (let j = 0; j <= N_GRID; j++) {
const lon = j*2*Math.PI/N_GRID;
pts.push([Math.cos(phi)*Math.cos(lon), Math.cos(phi)*Math.sin(lon), Math.sin(phi)]);
}
parallels.push(pts);
}
function projectPt(p) {
const rp = rotatePoint(p, yaw, pitch);
return project(rp, CX, CY, R);
}
function drawCurve(points3D, color, width, dashPattern) {
if (points3D.length < 2) return;
ctx.save();
ctx.strokeStyle = color;
ctx.lineWidth = width;
if (dashPattern) ctx.setLineDash(dashPattern);
ctx.beginPath();
const proj = points3D.map(projectPt);
let firstVisible = true;
for (let i = 0; i < proj.length - 1; i++) {
const p1 = proj[i], p2 = proj[i+1];
const z1 = p1.z, z2 = p2.z;
if (z1 > 0 && z2 > 0) {
// Both front: draw segment
if (firstVisible) { ctx.moveTo(p1.x, p1.y); firstVisible = false; }
ctx.lineTo(p2.x, p2.y);
} else if (z1 > 0 && z2 <= 0) {
// p1 front, p2 back: clip to silhouette
const t = z1 / (z1 - z2);
const intX = p1.x + t*(p2.x - p1.x);
const intY = p1.y + t*(p2.y - p1.y);
if (firstVisible) { ctx.moveTo(p1.x, p1.y); firstVisible = false; }
ctx.lineTo(intX, intY);
ctx.stroke();
ctx.beginPath();
firstVisible = true;
} else if (z1 <= 0 && z2 > 0) {
// p1 back, p2 front: clip from silhouette
const t = -z1 / (z2 - z1);
const intX = p1.x + t*(p2.x - p1.x);
const intY = p1.y + t*(p2.y - p1.y);
ctx.moveTo(intX, intY);
ctx.lineTo(p2.x, p2.y);
firstVisible = false;
}
// Both back: skip
}
ctx.stroke();
ctx.restore();
}
function render() {
ctx.clearRect(0, 0, W, H);
// === Draw sphere fill (radial gradient for 3D appearance) ===
const grad = ctx.createRadialGradient(CX - R*0.3, CY - R*0.35, R*0.05, CX, CY, R);
grad.addColorStop(0, '#f8f9fc');
grad.addColorStop(0.6, '#dce4ee');
grad.addColorStop(0.92, '#b0bec5');
grad.addColorStop(1, '#78909c');
ctx.fillStyle = grad;
ctx.beginPath();
ctx.arc(CX, CY, R, 0, 2*Math.PI);
ctx.fill();
// === Draw grid lines ===
ctx.save();
ctx.beginPath();
ctx.arc(CX, CY, R, 0, 2*Math.PI);
ctx.clip(); // Clip to sphere circle
// Meridians
for (const mer of meridians) {
const pts = mer.map(projectPt);
ctx.strokeStyle = 'rgba(120,144,156,0.25)';
ctx.lineWidth = 0.6;
ctx.beginPath();
let first = true;
for (let i = 0; i < pts.length; i++) {
if (pts[i].z > 0) {
if (first) { ctx.moveTo(pts[i].x, pts[i].y); first = false; }
else ctx.lineTo(pts[i].x, pts[i].y);
} else first = true;
}
ctx.stroke();
}
// Parallels
for (const par of parallels) {
const pts = par.map(projectPt);
ctx.strokeStyle = 'rgba(120,144,156,0.25)';
ctx.lineWidth = 0.6;
ctx.beginPath();
let first = true;
for (let i = 0; i < pts.length; i++) {
if (pts[i].z > 0) {
if (first) { ctx.moveTo(pts[i].x, pts[i].y); first = false; }
else ctx.lineTo(pts[i].x, pts[i].y);
} else first = true;
}
ctx.stroke();
}
ctx.restore();
// === Draw equator (slightly thicker) ===
const equatorPts = [];
for (let j = 0; j <= N_GRID; j++) {
const lon = j*2*Math.PI/N_GRID;
equatorPts.push([Math.cos(lon), Math.sin(lon), 0]);
}
drawCurve(equatorPts, 'rgba(84,110,122,0.45)', 1.0, []);
// === Draw curves ===
// True curve
drawCurve(s2_fits.trueCurve, '#2b8a3e', 2.8, [8, 3]);
// Local-constant fit
drawCurve(s2_fits.lcCurve, '#f08c00', 2.2, []);
// Local-linear fit
drawCurve(s2_fits.llCurve, '#1971c2', 2.6, []);
// === Draw sphere outline ===
ctx.strokeStyle = 'rgba(55,71,79,0.35)';
ctx.lineWidth = 1.2;
ctx.beginPath();
ctx.arc(CX, CY, R, 0, 2*Math.PI);
ctx.stroke();
// === Draw data points ===
for (const yi of s2_data.Y) {
const pp = projectPt(yi);
if (pp.z <= 0) continue; // back face
// Fade points near silhouette
const alpha = 0.18 + 0.42*Math.min(1, pp.z/0.5);
ctx.fillStyle = `rgba(73,80,87,${alpha.toFixed(2)})`;
ctx.beginPath();
ctx.arc(pp.x, pp.y, 2.5, 0, 2*Math.PI);
ctx.fill();
}
// === Draw north pole marker ===
const np = projectPt([0, 0, 1]);
if (np.z > 0) {
ctx.fillStyle = '#546e7a';
ctx.beginPath();
ctx.arc(np.x, np.y, 2.5, 0, 2*Math.PI);
ctx.fill();
ctx.fillStyle = '#37474f';
ctx.font = '10px system-ui';
ctx.fillText('N', np.x + 6, np.y + 4);
}
}
// === Mouse/Touch interaction ===
canvas.addEventListener('mousedown', (e) => {
dragging = true;
lastMouse = {x: e.offsetX, y: e.offsetY};
canvas.style.cursor = 'grabbing';
});
canvas.addEventListener('mousemove', (e) => {
if (!dragging) return;
const dx = e.offsetX - lastMouse.x;
const dy = e.offsetY - lastMouse.y;
yaw += dx * 0.008;
// Match the globe's vertical rotation to the mouse movement.
pitch += dy * 0.008;
pitch = Math.max(-Math.PI/2 + 0.05, Math.min(Math.PI/2 - 0.05, pitch));
lastMouse = {x: e.offsetX, y: e.offsetY};
render();
});
canvas.addEventListener('mouseup', () => { dragging = false; canvas.style.cursor = 'grab'; });
canvas.addEventListener('mouseleave', () => { dragging = false; canvas.style.cursor = 'grab'; });
canvas.addEventListener('touchstart', (e) => {
e.preventDefault();
dragging = true;
lastMouse = {x: e.touches[0].clientX, y: e.touches[0].clientY};
}, {passive: false});
canvas.addEventListener('touchmove', (e) => {
e.preventDefault();
if (!dragging) return;
const dx = e.touches[0].clientX - lastMouse.x;
const dy = e.touches[0].clientY - lastMouse.y;
yaw += dx * 0.008;
// Match the globe's vertical rotation to the touch movement.
pitch += dy * 0.008;
pitch = Math.max(-Math.PI/2 + 0.05, Math.min(Math.PI/2 - 0.05, pitch));
lastMouse = {x: e.touches[0].clientX, y: e.touches[0].clientY};
render();
}, {passive: false});
canvas.addEventListener('touchend', () => { dragging = false; });
canvas.addEventListener('touchcancel', () => { dragging = false; });
render();
return container
}
```
```{ojs}
//| label: fig-local-poly
//| fig-cap: "Interactive: local-constant regression and a two-step tangent local-linear approximation on S²"
//| echo: false
html`<div style="width:100%;max-width:840px;margin:0 auto;">
${s2_controls_view}
${s2_plot}
</div>`
```
The demo plots data and pointwise fits directly on $S^2$. **Rotate the globe by dragging**. Orange is a numerical local-constant fit found by a descent/Karcher iteration; on the positively curved sphere this need not be the unique global Fréchet minimizer. Blue is a computationally light **two-step approximation**: it first finds the local Fréchet mean, maps nearby observations to that tangent plane, performs Euclidean local-linear regression there, and maps the fitted intercept back. It is not the joint minimizer of $\hat Q_{1,h}(\omega,v_1;x)$ defined above.
- Local linear fitting often reduces leading boundary bias under the usual smoothness, design-density, and bandwidth conditions; a finite noisy sample need not show lower MSE.
- Larger $\sigma$ raises variability, and estimating a tangent slope adds uncertainty.
- The displayed blue line joins pointwise fitted intercepts. It is not one global geodesic, because its base point and local tangent fit are recomputed at every grid value.
## Application: From Metric-Space to Manifold — Local Polynomial Regression on $\mathcal{S}_{++}^m$ {#sec-l19-application}
### Motivation: What the Manifold Adds to Regression
Lectures 6–9 developed a progressive toolkit for regressing SPD-valued responses on scalar predictors:
| Lecture | Method | What it captures |
|---|---|---|---|
| 6 | Global Fréchet regression | Linear trend in log-covariance |
| 7 | Kernel Fréchet regression | Smooth nonlinear trend |
| 8 | Local-linear Fréchet regression | Nonlinear trend + boundary correction |
| 9 | TV-regularized Fréchet regression | Piecewise-constant jumps (regime changes) |
All four methods treat $\mathcal{S}_{++}^m$ as a **metric space** — they only need pairwise distances. The predictor-side weight functions $s(z, x, h)$ do the linearization. With the **Riemannian manifold** structure from Lectures 15–16, we now have intrinsic local polynomial regression: instead of predictor-side weights, we linearize in the **tangent space** at the fitted intercept.
This section revisits the portfolio stress-testing application from Lectures 6–9 and asks: **what does the manifold approach add, numerically and interpretively?**
### The SPD Manifold as a Regression Space
Suppose $r_t\in\mathbb{R}^p$ is the vector of asset returns at time $t$. Divide the return series into estimation windows and, for each window $i$, construct
$$
(X_i,\Sigma_i),\qquad i=1,\ldots,n,
$$
where $X_i$ is a scalar market-state variable---for example, realized market volatility, a volatility-index-like stress measure, or a macroeconomic score---and $\Sigma_i\in\mathcal{S}_{++}^p$ is the estimated covariance matrix of the $p$ asset returns in that window. The windows, covariance estimator, and stress measurement should all use information available at the same forecast origin; otherwise the fitted risk curve can contain look-ahead leakage.
Under the **log-Euclidean (LE) metric**, the logarithm map turns the SPD manifold into the vector space $\operatorname{Sym}(p)$:
$$
Z_i=\log\Sigma_i\in\operatorname{Sym}(p),
\qquad
d_{\mathrm{LE}}(\Sigma_i,\Sigma_j)=\|Z_i-Z_j\|_F.
$$
Thus fitting SPD-valued responses is equivalent to fitting the symmetric-matrix responses $Z_i$ by ordinary local polynomial regression and then applying the matrix exponential. For software that expects vectors, use the isometric half-vectorization
$$
\operatorname{svec}(Z)
=\bigl(Z_{11},\ldots,Z_{pp},\sqrt{2}Z_{12},\ldots,\sqrt{2}Z_{p-1,p}\bigr)^\top,
$$
with any fixed ordering of the off-diagonal entries. The factor $\sqrt{2}$ ensures $\|\operatorname{svec}(Z)\|_2=\|Z\|_F$; simply stacking the upper triangle without this factor changes the geometry.
Let $L(x)=\log\mu(x)$ denote the log-covariance curve at stress level $x$. Around a chosen target $x$, approximate it by
$$
L(x+u)\approx A_0(x)+uA_1(x)+\frac{u^2}{2}A_2(x)+\cdots+\frac{u^\ell}{\ell!}A_\ell(x),
$$
where each $A_j(x)$ is symmetric. $A_0(x)$ estimates the log-covariance at $x$, $A_1(x)$ estimates its stress derivative, and $A_2(x)$ estimates how that log-coordinate slope changes with stress.
#### Pointwise fitting algorithm
At target stress $x$, set
$$
u_i=X_i-x,
\qquad
r_i=\left(1,u_i,\frac{u_i^2}{2!},\ldots,\frac{u_i^\ell}{\ell!}\right)^\top,
\qquad
k_i=K_h(u_i).
$$
Let $R_x$ be the $n\times(\ell+1)$ matrix with rows $r_i^\top$, let $W_x=\operatorname{diag}(k_1,\ldots,k_n)$, and let $\mathbf Z$ be the $n\times q$ matrix with rows $\operatorname{svec}(Z_i)^\top$, where $q=p(p+1)/2$. The log-Euclidean local polynomial estimator is the weighted least-squares solution
$$
\widehat C(x)
=\left(R_x^\top W_xR_x\right)^{-1}R_x^\top W_x\mathbf Z,
$$
provided $R_x^\top W_xR_x$ is well-conditioned. Applying $\operatorname{svec}^{-1}$ to row $j+1$ of $\widehat C(x)$ gives $\widehat A_j(x)$. The estimates on the SPD manifold are then
$$
\widehat\mu(x)=\exp\{\widehat A_0(x)\}
$$
For the common local-linear case ($\ell=1$), the calculation can also be written directly in matrix form. Define
$$
S_r=\sum_{i=1}^n k_i u_i^r,
\qquad
T_r=\sum_{i=1}^n k_i u_i^r Z_i,
\qquad r=0,1,2,
$$
where $T_r$ is a symmetric matrix, and let $D=S_0S_2-S_1^2$. Then
$$
\widehat A_0(x)=\frac{S_2T_0-S_1T_1}{D},
\qquad
\widehat A_1(x)=\frac{S_0T_1-S_1T_0}{D}.
$$
This formula makes the boundary correction transparent: near an edge of the observed stress range, the local-linear fit compensates for the one-sided predictor design instead of merely averaging the nearby log-covariances.
#### From the fitted covariance curve to portfolio risk
Let $a\in\mathbb{R}^p$ be a **fixed** vector of portfolio weights. Once $\widehat\Sigma(x)=\widehat\mu(x)$ has been obtained, the fitted conditional portfolio variance and volatility are
$$
\widehat V_a(x)=a^\top\widehat\Sigma(x)a,
\qquad
\widehat\sigma_a(x)=\sqrt{\widehat V_a(x)}.
$$
These scalar curves answer the operational question: *what risk would the current portfolio carry if the market were at stress level $x$?* The local slope also quantifies the rate at which that risk changes. For $\ell\geq1$,
$$
\widehat{\dot\Sigma}(x)
:=\left.\frac{d}{dz}\widehat\mu_x(z)\right|_{z=x}
=D\exp_{\widehat A_0(x)}[\widehat A_1(x)],
$$
and therefore
$$
\widehat{\dot V}_a(x)
=a^\top\widehat{\dot\Sigma}(x)a.
$$
::: {.callout-note title="Why is the covariance slope $D\exp_{\widehat A_0(x)}[\widehat A_1(x)]$ rather than $\widehat A_1(x)$?" collapse="true"}
$\widehat A_1(x)$ is the slope of the **log-covariance**, not the covariance itself. The locally fitted curve is
$$
\widehat\mu_x(z)
=\exp\!\left\{\widehat A_0(x)+(z-x)\widehat A_1(x)\right\},
$$
so the chain rule gives
$$
\left.\frac{d}{dz}\widehat\mu_x(z)\right|_{z=x}
=D\exp_{\widehat A_0(x)}[\widehat A_1(x)].
$$
The scalar analogue is $\Sigma(x)=e^{L(x)}$, for which $\Sigma'(x)=e^{L(x)}L'(x)$ rather than merely $L'(x)$. For matrices the conversion is more involved because matrix multiplication need not commute:
$$
D\exp_A[H]
=\int_0^1 e^{(1-s)A}H e^{sA}\,ds.
$$
If $A$ and $H$ commute, this reduces to $D\exp_A[H]=e^A H$; it equals $H$ only in special cases such as $A=0$.
For example, if
$$
\widehat A_0=
\begin{pmatrix}\log 4&0\\0&\log 9\end{pmatrix},
\qquad
\widehat A_1=
\begin{pmatrix}0.1&0\\0&0.2\end{pmatrix},
$$
then the two matrices commute and
$$
\widehat{\dot\Sigma}
=D\exp_{\widehat A_0}[\widehat A_1]
=\begin{pmatrix}4(0.1)&0\\0&9(0.2)\end{pmatrix}
=\begin{pmatrix}0.4&0\\0&1.8\end{pmatrix},
$$
not $\widehat A_1$. Thus $\widehat A_1$ describes motion in global log-covariance coordinates, while $D\exp_{\widehat A_0}[\widehat A_1]$ expresses the same motion on the original covariance scale used in portfolio-risk sensitivity $a^\top\widehat{\dot\Sigma}(x)a$.
:::
#### Choosing the bandwidth and validating the risk curve
Choose $h$ using leave-one-out or blocked cross-validation with the intrinsic prediction loss
$$
\operatorname{CV}(h)
=\sum_i d_{\mathrm{LE}}^2\!\left(
\Sigma_i,\widehat\mu_{-i,h}(X_i)
\right).
$$
For overlapping rolling covariance windows, ordinary random folds leak many of the same returns into training and validation sets. Time-ordered or blocked folds, preferably with a gap at least as long as the overlap, give a more honest assessment. The bandwidth selected to predict covariance matrices or portfolio variance need not be optimal for estimating $\widehat A_1(x)$; derivative estimation is noisier and usually benefits from a separate, smoother bandwidth choice.
### Interactive Demo: Local-Constant vs. Local-Linear on Simulated SPD Data
The following demo simulates $2 \times 2$ SPD covariance matrices whose log-entries depend nonlinearly on a stress predictor, mimicking the portfolio stress-testing setup. We compare:
1. **Local-constant** ($\ell = 0$, kernel Fréchet) — the metric-space approach from Lecture 7.
2. **Local-linear** ($\ell = 1$, intrinsic) — the manifold approach from this lecture.
3. **Boundary behavior** — local-linear should reduce bias near the edges ($x \approx 0$ and $x \approx 1$).
4. **Estimated slopes** — the manifold method recovers the true tangent-space derivative.
```{ojs}
//| echo: false
//| output: false
spd_n19_control = Inputs.range([30, 200], {step: 10, value: 80, label: "Sample size n"})
spd_h19_control = Inputs.range([0.08, 0.5], {step: 0.02, value: 0.18, label: "Bandwidth h"})
spd_noise19_control = Inputs.range([0.1, 1.2], {step: 0.05, value: 0.4, label: "Noise level σ"})
spd_seed19_control = Inputs.range([1, 20], {step: 1, value: 7, label: "Random seed"})
spd_n19 = Generators.input(spd_n19_control)
spd_h19 = Generators.input(spd_h19_control)
spd_noise19 = Generators.input(spd_noise19_control)
spd_seed19 = Generators.input(spd_seed19_control)
spdlp_controls_view = {
for (const control of [spd_n19_control, spd_h19_control, spd_noise19_control, spd_seed19_control]) {
const numberInput = control.querySelector('input[type="number"]');
if (numberInput) numberInput.style.setProperty('width', '10rem', 'important');
}
return html`
<style>
.spdlp-demo { width:100%; max-width:none; font-family:system-ui,sans-serif; container-type:inline-size; }
.spdlp-slider-grid { display:flex; flex-wrap:wrap; gap:8px 20px; width:100%; max-width:840px; margin:0 auto 12px; font:0.85em system-ui,sans-serif; container-type:inline-size; }
.spdlp-slider-grid > * { min-width:0; margin:0; }
.spdlp-numeric-control { flex:1 1 calc((100% - 20px)/2); }
.spdlp-slider-grid input[type="number"] { width:10rem !important; }
.spdlp-panels { display:grid; grid-template-columns:minmax(0, 1fr) minmax(0, 1fr); gap:16px; align-items:start; width:100%; }
.spdlp-panel { min-width:0; width:100%; }
.spdlp-plot-svg { display:block; width:100%; max-width:100%; height:auto; }
@container (max-width:520px) { .spdlp-panels { grid-template-columns:minmax(0, 1fr); } }
@container (max-width:480px) { .spdlp-numeric-control { flex-basis:100%; } }
</style>
<div class="spdlp-slider-grid">
<div class="spdlp-numeric-control">${spd_n19_control}</div>
<div class="spdlp-numeric-control">${spd_h19_control}</div>
<div class="spdlp-numeric-control">${spd_noise19_control}</div>
<div class="spdlp-numeric-control">${spd_seed19_control}</div>
</div>`;
}
function runSPDLocalPoly(n, h, sigma, seed) {
const rng = (function(a) { return function() { a |= 0; a = a + 0x6D2B79F5 | 0; var t = Math.imul(a ^ a >>> 15, 1 | a); t = t + Math.imul(t ^ t >>> 7, 61 | t) ^ t; return ((t ^ t >>> 14) >>> 0) / 4294967296; }; })(seed || 137);
function rnorm() {
const u1 = rng(), u2 = rng();
return Math.sqrt(-2*Math.log(Math.max(u1,1e-10)))*Math.cos(2*Math.PI*u2);
}
// Dimension of Sym(2) = 3: matrix-log entries [a, b, c] for [[a,b],[b,c]]
// True mean curve in log-domain (nonlinear)
// a(x) = 0.3 + 0.8*x + 0.4*sin(3*x) (log-var1)
// b(x) = 0.1 + 0.6*x + 0.3*sin(2.5*x - 0.5) (log-cov)
// c(x) = 0.2 + 1.0*x + 0.5*sin(3.5*x + 0.3) (log-var2)
function trueLogMu(x) {
return [
0.3 + 0.8*x + 0.4*Math.sin(3*x),
0.1 + 0.6*x + 0.3*Math.sin(2.5*x - 0.5),
0.2 + 1.0*x + 0.5*Math.sin(3.5*x + 0.3)
];
}
// True derivative (tangent slope) at x
function trueDeriv(x) {
return [
0.8 + 1.2*Math.cos(3*x),
0.6 + 0.75*Math.cos(2.5*x - 0.5),
1.0 + 1.75*Math.cos(3.5*x + 0.3)
];
}
// Predictor: uniform on [0, 1]
const X = Array.from({length: n}, () => rng());
// Generate log(Sigma) with noise isotropic in the Frobenius metric.
// The off-diagonal coordinate appears twice in a symmetric matrix.
const Y = X.map(xi => {
const mu = trueLogMu(xi);
return mu.map((mj, j) => mj + sigma*rnorm()/(j === 1 ? Math.sqrt(2) : 1));
});
// Gaussian kernel
function gKern(u) { return Math.exp(-0.5*u*u)/Math.sqrt(2*Math.PI); }
// Evaluation grid
const nGrid = 60;
const xGrid = Array.from({length: nGrid+1}, (_, j) => j/nGrid);
// ---- Local-constant (ℓ=0, kernel Fréchet) ----
const lcFit = xGrid.map(x => {
const rw = X.map(xi => gKern((xi-x)/h));
const wS = rw.reduce((a,b) => a+b, 0);
const nw = wS > 1e-10 ? rw.map(w => w/wS) : Array(n).fill(1/n);
const m = [0, 0, 0];
for (let i = 0; i < n; i++) {
for (let j = 0; j < 3; j++) m[j] += nw[i]*Y[i][j];
}
return m;
});
// ---- Local-linear (ℓ=1, intrinsic tangent-space) ----
const llFit = xGrid.map(x => {
// Step 1: local-constant as initial guess
const rw = X.map(xi => gKern((xi-x)/h));
const wS = rw.reduce((a,b) => a+b, 0);
const nw0 = wS > 1e-10 ? rw.map(w => w/wS) : Array(n).fill(1/n);
const w0 = [0, 0, 0];
for (let i = 0; i < n; i++) {
for (let j = 0; j < 3; j++) w0[j] += nw0[i]*Y[i][j];
}
// Step 2: solve weighted least squares for intercept and slope
// For each coordinate j independently (LE metric → Euclidean in log-domain):
// minimize Σ w_i * (Y_ij - α_j - (X_i-x)*β_j)²
const slopes = [0, 0, 0];
const intercepts = [0, 0, 0];
for (let j = 0; j < 3; j++) {
let S0 = 0, S1 = 0, S2 = 0, Sy = 0, Sxy = 0;
for (let i = 0; i < n; i++) {
const wi = rw[i];
const dx = X[i] - x;
S0 += wi;
S1 += wi*dx;
S2 += wi*dx*dx;
Sy += wi*Y[i][j];
Sxy += wi*dx*Y[i][j];
}
const det = S0*S2 - S1*S1;
if (Math.abs(det) > 1e-12) {
intercepts[j] = (S2*Sy - S1*Sxy)/det;
slopes[j] = (S0*Sxy - S1*Sy)/det;
} else {
intercepts[j] = Sy/Math.max(S0, 1e-12);
slopes[j] = 0;
}
}
// The fitted value at z=x is the intercept (evaluation at u=0)
return {fit: intercepts, slope: slopes};
});
// Compute MSE, boundary MSE, and slope error
function sqDist(a, b) {
return (a[0]-b[0])**2 + 2*(a[1]-b[1])**2 + (a[2]-b[2])**2;
}
let mseLC = 0, mseLL = 0;
let mseLCbound = 0, mseLLbound = 0;
let nBound = 0;
let slopeErr = 0;
for (let k = 0; k <= nGrid; k++) {
const x = xGrid[k];
const truth = trueLogMu(x);
mseLC += sqDist(lcFit[k], truth);
const llF = llFit[k].fit;
mseLL += sqDist(llF, truth);
slopeErr += sqDist(llFit[k].slope, trueDeriv(x));
if (x < 0.15 || x > 0.85) {
mseLCbound += sqDist(lcFit[k], truth);
mseLLbound += sqDist(llF, truth);
nBound++;
}
}
mseLC /= (nGrid+1);
mseLL /= (nGrid+1);
mseLCbound /= nBound;
mseLLbound /= nBound;
slopeErr /= (nGrid+1);
return {
X, Y, xGrid, lcFit,
llFit: llFit.map(d => d.fit),
llSlopes: llFit.map(d => d.slope),
trueCurve: xGrid.map(x => trueLogMu(x)),
trueDeriv: xGrid.map(x => trueDeriv(x)),
mseLC, mseLL, mseLCbound, mseLLbound, slopeErr,
h, sigma, n
};
}
spdLpRes = runSPDLocalPoly(spd_n19, spd_h19, spd_noise19, spd_seed19);
spd_plot = html`
<div class="spdlp-panels">
<!-- Plot of fits -->
<div class="spdlp-panel">
<h4>Fitted Curves: Local-Constant vs. Local-Linear on 𝓢<sub>++</sub><sup>2</sup></h4>
<svg class="spdlp-plot-svg" viewBox="0 0 460 280" width="460" height="280" preserveAspectRatio="xMidYMid meet" style="border: 1px solid #dee2e6; border-radius: 4px;">
${(() => {
const mg = {top:15, right:110, bottom:35, left:55};
const pW = 460-mg.left-mg.right, pH = 280-mg.top-mg.bottom;
const xS = x => mg.left + x*pW;
// Pick entry to display — show all 3 entries using different rows
// For compactness, show entry 0 (log-var1) as representative
const yMin = -1.2, yMax = 2.5;
const yS = y => mg.top + pH - (y-yMin)/(yMax-yMin)*pH;
function path(curve, entryIdx) {
return curve.map((pt,i) => `${i===0?'M':'L'} ${xS(spdLpRes.xGrid[i])} ${yS(pt[entryIdx])}`).join(' ');
}
const entry = 0; // log-var1
const entryLabel = '(log Σ)₁₁';
const scat = spdLpRes.X.map((xi,i) =>
`<circle cx="${xS(xi)}" cy="${yS(spdLpRes.Y[i][entry])}" r="1.8" fill="#adb5bd" opacity="0.45"/>`
).join('');
return `
<line x1="${mg.left}" y1="${mg.top}" x2="${mg.left}" y2="${mg.top+pH}" stroke="#adb5bd"/>
<line x1="${mg.left}" y1="${mg.top+pH}" x2="${mg.left+pW}" y2="${mg.top+pH}" stroke="#adb5bd"/>
<text x="${mg.left+pW/2}" y="${mg.top+pH+25}" text-anchor="middle" font-size="10">Stress predictor X (e.g., SPY volatility)</text>
<text x="${mg.left-45}" y="${mg.top+pH/2}" text-anchor="middle" font-size="10" transform="rotate(-90,${mg.left-45},${mg.top+pH/2})">${entryLabel}</text>
<path d="${path(spdLpRes.lcFit, entry)}" fill="none" stroke="#f08c00" stroke-width="2"/>
<path d="${path(spdLpRes.llFit, entry)}" fill="none" stroke="#1971c2" stroke-width="2.5"/>
<path d="${path(spdLpRes.trueCurve, entry)}" fill="none" stroke="#2b8a3e" stroke-width="1.5" stroke-dasharray="5,3"/>
${scat}
<rect x="${mg.left}" y="${mg.top}" width="${pW}" height="${pH}" fill="none" stroke="#dee2e6"/>
<g transform="translate(${mg.left+pW+5},${mg.top+5})">
<line x1="0" y1="5" x2="15" y2="5" stroke="#f08c00" stroke-width="2"/><text x="18" y="10" font-size="9">LC (ℓ=0)</text>
<line x1="0" y1="20" x2="15" y2="20" stroke="#1971c2" stroke-width="2.5"/><text x="18" y="25" font-size="9">LL (ℓ=1)</text>
<line x1="0" y1="35" x2="15" y2="35" stroke="#2b8a3e" stroke-width="1.5" stroke-dasharray="5,3"/><text x="18" y="40" font-size="9">True μ(x)</text>
</g>`;
})()}
</svg>
<p style="font-size:0.8em;color:#868e96;margin-top:2px;">
Displaying the (1,1) entry of the matrix logarithm. Compare the reported MSEs; either fit can win in a finite sample.
</p>
</div>
<!-- Diagnostics -->
<div class="spdlp-panel">
<div style="padding: 12px; background: #f8f9fa; border-radius: 6px; margin-bottom: 12px;">
<h4 style="margin-top: 0;">Performance Comparison</h4>
<table style="width: 100%; border-collapse: collapse; font-size: 0.88em;">
<tr style="border-bottom: 2px solid #dee2e6;">
<th style="padding: 3px 8px; text-align: left;">Metric</th>
<th style="padding: 3px 8px; text-align: center;">LC (ℓ=0)</th>
<th style="padding: 3px 8px; text-align: center;">LL (ℓ=1)</th>
<th style="padding: 3px 8px; text-align: center;">Δ</th>
</tr>
<tr style="border-bottom: 1px solid #e9ecef;">
<td style="padding: 3px 8px;">Overall MSE</td>
<td style="padding: 3px 8px; text-align: center; font-family: monospace;">${spdLpRes.mseLC.toFixed(4)}</td>
<td style="padding: 3px 8px; text-align: center; font-family: monospace; font-weight: bold;">${spdLpRes.mseLL.toFixed(4)}</td>
<td style="padding: 3px 8px; text-align: center; font-family: monospace; color: ${spdLpRes.mseLL<spdLpRes.mseLC?'#2b8a3e':'#c92a2a'};">${(spdLpRes.mseLL-spdLpRes.mseLC).toFixed(4)}</td>
</tr>
<tr style="border-bottom: 1px solid #e9ecef;">
<td style="padding: 3px 8px;">Boundary MSE</td>
<td style="padding: 3px 8px; text-align: center; font-family: monospace;">${spdLpRes.mseLCbound.toFixed(4)}</td>
<td style="padding: 3px 8px; text-align: center; font-family: monospace; font-weight: bold;">${spdLpRes.mseLLbound.toFixed(4)}</td>
<td style="padding: 3px 8px; text-align: center; font-family: monospace; color: ${spdLpRes.mseLLbound<spdLpRes.mseLCbound?'#2b8a3e':'#c92a2a'};">${(spdLpRes.mseLLbound-spdLpRes.mseLCbound).toFixed(4)}</td>
</tr>
<tr>
<td style="padding: 3px 8px;">Slope MSE</td>
<td style="padding: 3px 8px; text-align: center; color: #868e96; font-size: 0.83em;">N/A</td>
<td style="padding: 3px 8px; text-align: center; font-family: monospace; font-weight: bold;">${spdLpRes.slopeErr.toFixed(4)}</td>
<td style="padding: 3px 8px; text-align: center; font-size: 0.83em; color: #2b8a3e;">No fitted LC slope</td>
</tr>
</table>
</div>
<div style="padding: 12px; background: #f1f3f5; border-radius: 6px; margin-bottom: 12px;">
<h4 style="margin-top: 0;">Estimated Slopes at Selected Stress Levels</h4>
${(() => {
const idxs = [0, Math.floor(spdLpRes.xGrid.length/2), spdLpRes.xGrid.length-1];
const labels = ['Low stress (x≈0)', 'Mid stress (x≈0.5)', 'High stress (x≈1)'];
return idxs.map((idx, k) => {
const s = spdLpRes.llSlopes[idx];
const t = spdLpRes.trueDeriv[idx];
const xv = spdLpRes.xGrid[idx];
return `<div style="margin-bottom: 4px; font-size: 0.83em;">
<b>${labels[k]}:</b> Â₁ = [${s[0].toFixed(2)}, ${s[1].toFixed(2)}, ${s[2].toFixed(2)}]
<span style="color: #868e96; font-size: 0.88em;"> (true: [${t[0].toFixed(2)}, ${t[1].toFixed(2)}, ${t[2].toFixed(2)}])</span>
</div>`;
}).join('');
})()}
<p style="font-size: 0.78em; color: #868e96; margin-top: 6px;">
Entries: [∂(log Σ)₁₁/∂x, ∂(log Σ)₁₂/∂x, ∂(log Σ)₂₂/∂x], where log is the matrix logarithm.
</p>
</div>
</div>
</div>
`
```
```{ojs}
//| label: fig-spd-local-poly
//| fig-cap: "Interactive: Local-constant vs. local-linear regression on the SPD manifold (log-Euclidean metric)"
//| echo: false
html`<div class="spdlp-demo">${spdlp_controls_view}${spd_plot}</div>`
```
::: {.callout-note title="How to read this"}
- **Boundary MSE** averages squared LE distance for $x < 0.15$ or $x > 0.85$; local linear removes the leading design-asymmetry term under standard smoothness conditions.
- The **slope MSE** assesses the directly fitted derivative of the matrix-log curve; a local-constant fit has no slope parameter.
- Increasing **bandwidth $h$** smooths both fits and generally trades variance for bias.
- Increase **noise $\sigma$**: advantage of LL diminishes — with very noisy data, the extra parameter (slope) is harder to estimate.
:::
::: {.callout-tip title="Key observations from the demo"}
- **Boundary behavior:** Under standard smoothness and design conditions, local linear fitting removes the leading $O(h)$ boundary term of local constant fitting. The displayed finite-sample MSEs can still reverse because variance also matters.
- **Slope estimation:** The fitted slope estimates $L'(x)$ for $L(x)=\log\mu(x)$. Convert it with $D\exp_{L(x)}$ before interpreting changes in the entries of $\mu(x)$ itself.
- **Bias-variance tradeoff:** With small $h$ or large $\sigma$, slope estimates become noisy. Distance-loss cross-validation can select $h$ for prediction, while derivative estimation may favor a different bandwidth.
:::
### Discussion: What the Manifold Approach Adds to Regression
**1. Explicit slopes.** A basic local Fréchet fit produces $\hat\mu(x)$ without a response-space slope parameter. Intrinsic local polynomial regression estimates the intercept and slope jointly. Under LE geometry, $\hat A_1(x)$ is the derivative of the matrix-log curve; $D\exp_{\hat A_0(x)}[\hat A_1(x)]$ is the corresponding covariance sensitivity.
**2. Boundary correction without a boundary-specific estimator.** Under the regularity conditions used in local polynomial theory, local-constant bias is $O(h)$ at a boundary while local-linear bias is $O(h^2)$ [@YuanZhuLinMarron2012]. This concerns the leading asymptotic bias; it does not guarantee lower finite-sample MSE.
**3. Higher-order fits.** The manifold framework supports $\ell = 2$ (local quadratic), which estimates a second normal-coordinate coefficient. This describes bending of the regression curve; it is distinct from curvature of the ambient manifold. Metric-only methods can still represent nonlinear mean curves, but do not parameterize this intrinsic second-order term directly.
## Key Takeaways
- **Intrinsic local polynomial regression** uses the manifold's smooth structure (exp/log) to fit polynomials directly in tangent spaces.
- The **intrinsic error model** $Y = \exp_{\mu(X)}(\varepsilon)$, $\mathbb{E}(\varepsilon \mid X) = 0$ generalizes the Euclidean mean-zero error.
- **Order $\ell = 0$** recovers local-constant (kernel Fréchet); **$\ell = 1$** gives local-linear with explicit tangent slope.
- Contrast with **local Fréchet regression**: manifold method uses tangent-space linearization; metric method uses predictor-side weights.
- Under the LE metric, computation is weighted least squares. For affine-invariant local polynomials, use a carefully derived Riemannian optimizer or a stochastic search such as AESAMC.
- **On $\mathcal{S}_{++}^m$ with the LE metric**, local-linear regression reduces to weighted least squares in $\mathbb{R}^{m(m+1)/2}$ — no iterative manifold optimization required.
- The fitted slope provides first-order sensitivity directly; a basic metric-space Fréchet fit has no such direct parameter.
- Local-linear fitting reduces leading boundary bias under regularity conditions.
## Exercises
1. **Recovering kernel Fréchet.** Show that for $\ell = 0$, the intrinsic local polynomial estimator is exactly the kernel Fréchet regression estimator from Lecture 7. <a href="javascript:void(0)" onclick="showSolution('l19-sol-1')" class="solution-link">📝 Show Solution</a>
2. **Tangent slope interpretation.** For $\ell = 1$ on $S^1$, interpret $v_1$ geometrically. How does it relate to the derivative of the unwrapped mean curve? <a href="javascript:void(0)" onclick="showSolution('l19-sol-2')" class="solution-link">📝 Show Solution</a>
3. **Metric vs. manifold.** In what scenarios would you prefer local Fréchet regression over intrinsic local polynomial regression? <a href="javascript:void(0)" onclick="showSolution('l19-sol-3')" class="solution-link">📝 Show Solution</a>
4. **Affine-invariant local-linear.** Implement the affine-invariant local-linear estimator from @YuanZhuLinMarron2012.
<style>
.solution-link { font-size: 0.9em; text-decoration: none; white-space: nowrap; margin-left: 0.3em; }
.solution-link:hover { text-decoration: underline; }
.solution-dialog { padding: 0; max-width: 720px; }
.solution-dialog-header { display: flex; justify-content: space-between; align-items: flex-start; border-bottom: 1px solid #dee2e6; padding: 1.25rem 1.5rem 1rem; background: #f8f9fa; border-radius: 8px 8px 0 0; }
.solution-dialog-header h4 { margin: 0; font-size: 1.15rem; }
.solution-dialog-close { background: none; border: 1px solid #adb5bd; border-radius: 4px; padding: 0.2rem 0.75rem; cursor: pointer; font-size: 0.9rem; color: #495057; white-space: nowrap; flex-shrink: 0; }
.solution-dialog-close:hover { background: #e9ecef; }
.solution-original { padding: 1rem 1.5rem; background: #f1f3f5; border-left: 4px solid #868e96; margin: 1rem 1.5rem; border-radius: 4px; font-size: 0.95rem; }
.solution-answer { padding: 0.5rem 1.5rem 1.5rem; }
.solution-answer strong { color: #2b8a3e; }
dialog { border: none; border-radius: 8px; box-shadow: 0 8px 32px rgba(0,0,0,0.22); padding: 0; max-width: 750px; width: 90vw; }
dialog::backdrop { background: rgba(0,0,0,0.45); }
</style>
<dialog id="l19-sol-1"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 1</h4><button onclick="closeSolution('l19-sol-1')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Show ℓ=0 recovers kernel Fréchet regression.</div><div class="solution-answer"><strong>Solution:</strong> For $\ell=0$, $m_{\theta,x}(z) = \exp_\omega(0) = \omega$ for all $z$. The parameter $\theta = (\omega)$ with no $v_j$. The objective is $\hat{Q}_{0,h}(\omega; x) = \frac{1}{n}\sum K_h(X_i-x)d^2(Y_i, \omega)$, and the minimizer is $\hat{\mu}_{0,h}(x) = \arg\min_\omega \frac{1}{n}\sum K_h(X_i-x)d^2(Y_i, \omega)$. This is exactly the kernel Fréchet regression estimator with unnormalized kernel weights. After normalizing by $\sum K_h(X_j-x)$, it's identical to the estimator in Lecture 7.</div></div></dialog>
<dialog id="l19-sol-2"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 2</h4><button onclick="closeSolution('l19-sol-2')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Interpret v₁ geometrically on S¹.</div><div class="solution-answer"><strong>Solution:</strong> On $S^1$, the tangent space at $\omega$ is $\mathbb{R}$ (identified with the real line). The local linear curve is $m(z) = (\omega + (z-x)v_1) \bmod 2\pi$ (for small $(z-x)v_1$). $v_1$ is the **angular derivative** of the mean curve at $x$: $v_1 = \frac{d}{dz}\mu(z)|_{z=x}$, measuring how fast (in radians per unit $x$) the conditional Fréchet mean moves along the circle. Positive $v_1$ = counterclockwise movement; negative = clockwise. Unlike Euclidean local-linear regression where slopes are vectors in $\mathbb{R}^p$, here the slope $v_1$ is a single scalar — the rate of angular change.</div></div></dialog>
<dialog id="l19-sol-3"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 3</h4><button onclick="closeSolution('l19-sol-3')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> When prefer metric vs. manifold local regression?</div><div class="solution-answer"><strong>Solution:</strong> **Prefer local Fréchet regression (metric-only) when:** (1) The response space is not a smooth manifold (e.g., BHV tree space with singularities, or a general metric space). (2) You do not have exp/log maps—only distances. (3) You need only the fitted mean curve. **Prefer intrinsic local polynomial regression when:** (1) the response is a smooth Riemannian manifold with computable exp/log maps; (2) you want an explicit tangent slope; or (3) you want explicit higher-order normal-coordinate coefficients. Both approaches can fit nonlinear mean curves; the distinction is whether response-space derivatives are parameters of the fit.</div></div></dialog>
<dialog id="l19-sol-4"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 4</h4><button onclick="closeSolution('l19-sol-4')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Bandwidth selection for manifold local polynomial regression.</div><div class="solution-answer"><strong>Solution:</strong> Bandwidth still controls a bias–variance tradeoff, but response geometry enters the expansion through log/exp derivatives and the Hessian of the population criterion. There is no rate that follows merely from saying the manifold is Hadamard; it also depends on predictor dimension, polynomial order, design density, smoothness, and the estimator. A practical choice is cross-validation with the same squared geodesic loss used for prediction. If the goal is slope estimation rather than prediction, tune for derivative loss when a validation proxy or justified plug-in rule is available. Under the LE metric the problem is Euclidean in matrix-log coordinates, so standard local-polynomial bandwidth theory applies in those coordinates.</div></div></dialog>
<dialog id="l19-sol-5"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 5</h4><button onclick="closeSolution('l19-sol-5')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Interpret log-coordinate slopes for portfolio stress-testing.</div><div class="solution-answer"><strong>Solution:</strong> $\hat A_1(x)$ estimates $\frac{d}{dx}\log\mu(x)$, so comparing $\hat A_1(0.2)$ with $\hat A_1(0.8)$ compares the direction and speed of the fitted curve in the common LE coordinate space. Its off-diagonal entry is the derivative of an off-diagonal entry of the <em>matrix logarithm</em>; it is not automatically the derivative of covariance or correlation. The covariance-matrix sensitivity is
$$
\frac{d}{dx}\mu(x)=D\exp_{\log\mu(x)}[A_1(x)].
$$
Only after computing this Fréchet derivative (and, for correlation, applying the variance normalization and product rule) should one make an entrywise financial interpretation. A basic local-constant fit has no directly estimated slope parameter, although a derivative could be approximated from a separately smoothed curve.</div></div></dialog>
<dialog id="l19-sol-6"><div class="solution-dialog"><div class="solution-dialog-header"><h4>Exercise 6</h4><button onclick="closeSolution('l19-sol-6')" class="solution-dialog-close">✕ Close</button></div><div class="solution-original"><strong>Exercise:</strong> Derive boundary bias of ℓ=0 and show how ℓ=1 eliminates it.</div><div class="solution-answer"><strong>Solution:</strong> **Local-constant boundary bias (ℓ=0):** At a boundary point $x$ (say $x=0$), data only exists on one side ($X_i \ge 0$). The tangent-space expansion gives $\log_{\mu(x)}(\mu(X_i)) = (X_i - x)\dot{\mu}(x) + O((X_i-x)^2)$. The local-constant estimator solves $\min_\omega \sum K_h(X_i-x) d^2(Y_i, \omega)$. In the tangent space:
$$
\mathbb{E}\!\left[\log_{\mu(x)}\hat{\mu}_{0,h}(x)\right] \approx h\,\dot{\mu}(x)\frac{\int_0^\infty u K(u)\,du}{\int_0^\infty K(u)\,du} + O(h^2),
$$
which is $O(h)$ because $\int_0^\infty u K(u) du \neq 0$ (the kernel is asymmetric at the boundary). At interior points, the integral is symmetric and the $O(h)$ term vanishes, leaving $O(h^2)$ bias.
**Local-linear correction (ℓ=1):** The estimator fits $\min_{\omega, v_1} \sum K_h(X_i-x) d^2(Y_i, \exp_\omega((X_i-x)v_1))$. In the tangent space, this is weighted least squares with an intercept and slope. The design matrix includes the $(X_i-x)$ term, which **absorbs the first-order bias**:
$$
\mathbb{E}\!\left[\log_{\mu(x)}\hat{\mu}_{1,h}(x)\right] = O(h^2)
$$
at an interior or boundary target, under the usual smoothness, positive design-density, kernel-moment, and bandwidth conditions. This is a bias statement, not a guarantee about finite-sample MSE. **On the SPD manifold with the LE metric**, the calculation is exactly Euclidean in matrix-log coordinates because $(\mathcal{S}_{++}^m,g^{\mathrm{LE}})$ is isometric to a vector space.</div></div></dialog>
<script>
function showSolution(id) { const d = document.getElementById(id); if (d) { d.showModal(); if (window.MathJax && MathJax.typesetPromise) MathJax.typesetPromise([d]).catch(function(e) { console.log('MathJax error:', e); }); } }
function closeSolution(id) { const d = document.getElementById(id); if (d) d.close(); }
document.addEventListener('click', function(e) { if (e.target.tagName === 'DIALOG') e.target.close(); });
</script>
## Further Reading
- @YuanZhuLinMarron2012 — Local polynomial regression on SPD matrices.
- @Fletcher2013 — Geodesic regression (the parametric counterpart).
- @PetersenMueller2019 — Local Fréchet regression (the metric-only counterpart).
- @Liang2011AESAMC — AESAMC algorithm used for optimization.
- @BhattacharyaPatrangenaru2005 — Large-sample theory for ordinary intrinsic/extrinsic sample means; useful background, but not a local-polynomial confidence-band theorem.
- @ArsignyEtAl2007 — The log-Euclidean metric and its global matrix-log isometry.
- @PennecFillardAyache2006 — Affine-invariant SPD geometry and tensor-computing operations.
- @Moakher2005GeometricMean — Affine-invariant geometric mean on SPD matrices; alternative to LE metric for regression.
## Self-Assessment Quiz
Test your understanding of this lecture with the interactive MCQ quiz:
👉 **[Lecture 19 Quiz — 10 Multiple-Choice Questions](../quizzes/lecture-19-quiz.qmd)**