---
title: "Lecture 7: Fréchet Regression — Kernel Regression"
subtitle: "Nadaraya–Watson smoothing for metric-space 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:
- Explain why global Fréchet regression may be misspecified for nonlinear predictor–response relationships.
- Formulate the kernel Fréchet regression target (local-constant) as a Nadaraya–Watson extension to metric spaces.
- Distinguish between the kernel target $\mu_h^{\mathrm{NW}}(x)$ and the true conditional Fréchet mean $\mu(x)$.
- Derive the normalized-kernel-weight representation of the sample estimator.
- Apply the split proximal-point algorithms (cyclic and randomized) from Lecture 5 to compute the kernel Fréchet regression estimator in a Hadamard space, using the closed-form component resolvent $t_i = 2\lambda w_{in}^{\mathrm{NW}} / (1 + 2\lambda w_{in}^{\mathrm{NW}})$.
- Discuss bandwidth selection and the bias–variance tradeoff in the metric-space context.
## From Global to Local: Why Kernel Smoothing?
Global Fréchet regression (Lecture 6) is the analogue of a **linear model**. Like any linear model, it may be misspecified when the relation between predictor and response is nonlinear. To address this, one can generalize nonparametric regression — specifically kernel regression — to general metric spaces.
The most direct local smoother for the conditional Fréchet mean is the metric-space analogue of the classical **Nadaraya–Watson** kernel estimator. This local-constant point of view is historically important: for manifold-valued and shape-valued responses under random design, the regression estimator of @DavisFletcherBullittJoshi2007 is obtained by replacing the ordinary weighted Euclidean average in kernel regression with a weighted Fréchet mean. The extension to general metric spaces is straightforward.
## Kernel Setup
Let $K : \mathbb{R}^p \to [0, \infty)$ be a kernel function satisfying
$$
\int_{\mathbb{R}^p} K(u)\, du = 1,
$$
and for bandwidth $h > 0$ define the rescaled kernel
$$
K_h(u) = h^{-p} K(u / h), \qquad u \in \mathbb{R}^p.
$$
We assume $\mathbb{E}\{K_h(X - x)\} > 0$ for the predictor value $x$ under consideration.
## The Kernel Fréchet Regression Target
First recall the ordinary Euclidean case $Y \in \mathbb{R}$. The population Nadaraya–Watson smoother is
$$
m_h(x) = \frac{\mathbb{E}\{K_h(X - x)Y\}}{\mathbb{E}\{K_h(X - x)\}}.
$$
Since the denominator does not depend on the optimization variable, this same quantity can be written in weighted least-squares form as
$$
m_h(x) = \arg\min_{y \in \mathbb{R}} \mathbb{E}\bigl[K_h(X - x)(Y - y)^2\bigr].
$$
This reformulation is the **key to the extension** beyond Euclidean responses: once the response enters only through squared distance, the linear averaging operation can be replaced by a Fréchet mean.
::: {#def-kernel-frechet .definition title="Kernel Fréchet regression (local-constant Fréchet regression)"}
For fixed bandwidth $h > 0$ and predictor value $x \in \mathbb{R}^p$, define the localized Fréchet objective
$$
M_h^{\mathrm{NW}}(\omega, x) = \mathbb{E}\bigl\{K_h(X - x) d^2(Y, \omega)\bigr\}, \qquad \omega \in \mathcal{M}.
$$
Any minimizer of $M_h^{\mathrm{NW}}(\cdot, x)$ is called the **kernel Fréchet regression target** or **local-constant Fréchet regression target** at $x$; when unique, we denote it by
$$
\mu_h^{\mathrm{NW}}(x) = \arg\min_{\omega \in \mathcal{M}} M_h^{\mathrm{NW}}(\omega, x).
$$
Equivalently, since $\mathbb{E}\{K_h(X - x)\}$ does not depend on $\omega$,
$$
\mu_h^{\mathrm{NW}}(x) = \arg\min_{\omega \in \mathcal{M}} \frac{\mathbb{E}\bigl\{K_h(X - x) d^2(Y, \omega)\bigr\}}{\mathbb{E}\{K_h(X - x)\}}.
$$
:::
Thus, for each $x$, kernel Fréchet regression is simply the **Fréchet mean** of the response distribution localized around $x$ by kernel weights. If $\mathcal{M} = \mathbb{R}^q$ with the Euclidean metric, the minimizer is the ordinary Nadaraya–Watson estimator componentwise — so this construction is a genuine extension rather than a new method in disguise.
## Bias–Variance Perspective
It is important to distinguish $\mu_h^{\mathrm{NW}}(x)$ from the true conditional Fréchet mean
$$
\mu(x) = \arg\min_{\omega \in \mathcal{M}} \mathbb{E}\{d^2(Y, \omega) \mid X = x\}.
$$
For fixed $h > 0$, the kernel target $\mu_h^{\mathrm{NW}}(x)$ is a **smoothed version** of $\mu(x)$, because it averages nearby conditional laws according to $K_h(X - x)$. Under the usual continuity assumptions on the joint law of $(X, Y)$, one has $\mu_h^{\mathrm{NW}}(x) \to \mu(x)$ as $h \to 0$, but at finite bandwidth the two are not identical.
::: {.callout-note title="Bias of the kernel target"}
This is exactly the same bias phenomenon that occurs for ordinary local-constant kernel regression. The bias depends on the curvature of $\mu(x)$ and the bandwidth $h$. For Euclidean responses and twice-differentiable $\mu$, the bias is $O(h^2)$. In a general metric space, the bias rate depends on the local geometry of $\mathcal{M}$ near $\mu(x)$.
:::
## Sample Estimator
Given data $(X_1, Y_1), \ldots, (X_n, Y_n)$, the sample analogue is immediate. Define the empirical local objective
$$
\hat{M}_{n,h}^{\mathrm{NW}}(\omega, x) = \frac{1}{n}\sum_{i=1}^n K_h(X_i - x) d^2(Y_i, \omega), \qquad \omega \in \mathcal{M},
$$
and, provided $\sum_{j=1}^n K_h(X_j - x) > 0$, the normalized kernel weights
$$
w_{in}^{\mathrm{NW}}(x, h) = \frac{K_h(X_i - x)}{\sum_{j=1}^n K_h(X_j - x)}, \qquad i = 1, \ldots, n.
$$
Then the sample objective can be written equivalently as
$$
\hat{M}_{n,h}^{\mathrm{NW}}(\omega, x) \propto \sum_{i=1}^n w_{in}^{\mathrm{NW}}(x, h) d^2(Y_i, \omega),
$$
so the **kernel Fréchet regression estimator** is
$$
\hat{\mu}_h^{\mathrm{NW}}(x) = \arg\min_{\omega \in \mathcal{M}} \hat{M}_{n,h}^{\mathrm{NW}}(\omega, x)
= \arg\min_{\omega \in \mathcal{M}} \sum_{i=1}^n w_{in}^{\mathrm{NW}}(x, h) d^2(Y_i, \omega).
$$
Hence estimation at each target point $x$ reduces to a **weighted sample Fréchet mean** problem.
## Computation via the Proximal Point Algorithm {#sec-computation-ppa}
The computation is conceptually simple:
1. Choose a kernel $K$ and bandwidth $h$.
2. For the desired predictor value $x$, compute the local weights $w_{in}^{\mathrm{NW}}(x, h)$.
3. Minimize the weighted Fréchet objective $\sum_i w_{in}^{\mathrm{NW}}(x, h) d^2(Y_i, \omega)$ over $\omega \in \mathcal{M}$.
The third step depends on the geometry of the response space. We focus on **Hadamard spaces** (complete, simply connected metric spaces of nonpositive curvature). In such spaces, each summand $\omega \mapsto d^2(Y_i, \omega)$ is geodesically convex, and the kernel-weighted objective
$$
Q_x(\omega) = \sum_{i=1}^n w_{in}^{\mathrm{NW}}(x, h)\, d^2(Y_i, \omega)
\tag{1}
$$
inherits this convexity — the weights $w_{in}^{\mathrm{NW}}(x, h)$ are nonnegative and satisfy $\sum_{i=1}^n w_{in}^{\mathrm{NW}} = 1$, so $Q_x$ is a convex combination of convex functions. Consequently, $Q_x$ admits a unique minimizer, namely the kernel Fréchet regression estimator $\hat{\mu}_h^{\mathrm{NW}}(x)$.
Equation (1) is exactly a **weighted sample Fréchet mean** problem — identical in form to the objective studied in Lecture 5, with the kernel weights $w_{in}^{\mathrm{NW}}(x, h)$ playing the role of the Fréchet weights. The proximal-point algorithms developed there therefore apply **directly**. We recall the key constructions and show how they specialize to the NW setting.
### Proximal Maps and the Full PPA
As introduced in Lecture 5, the **proximal map** (or resolvent) of a convex function $f$ on a Hadamard space is
$$
J_{\lambda}^{f}(\omega) = \arg\min_{z \in \mathcal{M}} \left\{ f(z) + \frac{1}{2\lambda} d^2(z, \omega) \right\}, \qquad \lambda > 0.
$$
The added term $\frac{1}{2\lambda} d^2(z, \omega)$ penalizes movement away from the current point, rendering the subproblem strongly convex and its solution unique. The **full proximal-point algorithm (PPA)** for minimizing $Q_x$ iterates
$$
\omega_{k+1} = J_{\lambda_k}^{Q_x}(\omega_k), \qquad k = 0, 1, 2, \ldots,
\tag{2}
$$
starting from an arbitrary $\omega_0 \in \mathcal{M}$. Each step solves the **entire** regularized weighted Fréchet mean problem — conceptually clean, but often as hard as the original minimization. For practical computation, especially with the data-dependent NW weights, a **split** approach is preferred.
### Split PPA for the NW Estimator
Following the template of Lecture 5, split the NW objective (1) into its summands:
$$
Q_x = \sum_{i=1}^n f_i, \qquad
f_i(\omega) = w_{in}^{\mathrm{NW}}(x, h)\, d^2(Y_i, \omega).
$$
A proximal step is then applied to **one summand at a time**, avoiding the full minimization. Two standard schedules are available:
**Cyclic split PPA.** Cycle through the data in a fixed order. For cycle index $k = 0, 1, \ldots$:
$$
\begin{aligned}
\omega_{kn+1} &= J_{\lambda_k}^{f_1}(\omega_{kn}), \\
\omega_{kn+2} &= J_{\lambda_k}^{f_2}(\omega_{kn+1}), \\
&\;\;\vdots \\
\omega_{kn+n} &= J_{\lambda_k}^{f_n}(\omega_{kn+n-1}),
\end{aligned}
\tag{3}
$$
with the step-size parameter $\lambda_k$ held **fixed** within the $k$-th cycle.
**Randomized split PPA.** At each iteration $k$, pick an index $r_k \in \{1, \ldots, n\}$ (e.g., uniformly) and update
$$
\omega_{k+1} = J_{\lambda_k}^{f_{r_k}}(\omega_k).
\tag{4}
$$
Both variants are particularly natural for kernel regression: when $h$ is small, only a few observations carry substantial weight $w_{in}^{\mathrm{NW}}(x, h)$, so most $f_i$ are nearly zero and contribute negligible updates. The randomized scheme automatically samples the influential points more often (since they carry larger weights, one may also sample proportionally to $w_{in}^{\mathrm{NW}}$), while the cyclic scheme guarantees every data point is visited.
### Closed Form of the Component Resolvent
For the NW summand $f_i(\omega) = w_{in}^{\mathrm{NW}}(x, h)\, d^2(Y_i, \omega)$, the component resolvent has a **closed form** — a key result from Lecture 5. The subproblem
$$
J_{\lambda}^{f_i}(\omega) = \arg\min_{z \in \mathcal{M}} \left\{ w_{in}^{\mathrm{NW}}(x, h)\, d^2(Y_i, z) + \frac{1}{2\lambda} d^2(z, \omega) \right\}
$$
is solved by a point on the geodesic $[\omega, Y_i]$ connecting the current iterate to the $i$-th response:
$$
J_{\lambda}^{f_i}(\omega) = \gamma_{\omega \to Y_i}(t_i), \qquad
t_i = \frac{2\lambda\, w_{in}^{\mathrm{NW}}(x, h)}{1 + 2\lambda\, w_{in}^{\mathrm{NW}}(x, h)},
\tag{5}
$$
where $\gamma_{\omega \to Y_i}(t)$ denotes the point at fraction $t$ along the geodesic from $\omega$ to $Y_i$. The interpolation parameter $t_i$ lies in $(0, 1)$ for any $\lambda > 0$ and $w_{in}^{\mathrm{NW}} > 0$, guaranteeing that each split proximal step is a **geodesic interpolation** — move a fraction $t_i$ of the way toward data point $Y_i$.
::: {.callout-note title="NW-specific behavior of the component step"}
The step length $t_i$ depends on the product $\lambda \cdot w_{in}^{\mathrm{NW}}(x, h)$:
- **Large kernel weight** ($Y_i$ near $x$ in predictor space): $w_{in}^{\mathrm{NW}}$ is large, so $t_i$ is close to 1 — the algorithm takes a substantial step toward $Y_i$.
- **Small kernel weight** ($Y_i$ far from $x$): $w_{in}^{\mathrm{NW}}$ is small, so $t_i \approx 2\lambda w_{in}^{\mathrm{NW}}$ is tiny — the algorithm barely moves toward $Y_i$.
- **As $h \to 0$:** only the few observations with $X_i \approx x$ receive non-negligible weight; the split PPA automatically focuses on the relevant local data.
- **As $h \to \infty$:** all $w_{in}^{\mathrm{NW}} \to 1/n$, recovering the uniform-weight Fréchet mean computation of Lecture 5.
:::
### Algorithm Summary: Randomized Split PPA for Kernel Fréchet Regression
For a fixed evaluation point $x$ and bandwidth $h$:
**Step 1.** Compute the kernel weights $w_i = w_{in}^{\mathrm{NW}}(x, h)$ for $i = 1, \ldots, n$.
**Step 2.** Choose an initial point $\omega_0 \in \mathcal{M}$ (e.g., the unconditional Fréchet mean of $\{Y_i\}$, or any data point $Y_i$ with large $w_i$).
**Step 3.** Choose a step-size sequence $(\lambda_k)$ satisfying
$$
\sum_{k=0}^{\infty} \lambda_k = \infty, \qquad \sum_{k=0}^{\infty} \lambda_k^2 < \infty.
\tag{6}
$$
A standard choice is $\lambda_k = c / (k+1)$ for some $c > 0$.
**Step 4.** At iteration $k = 0, 1, 2, \ldots$:
- Select an index $r_k \in \{1, \ldots, n\}$ (uniformly at random, or with probability proportional to $w_i$).
- Compute $t_k = \dfrac{2\lambda_k\, w_{r_k}}{1 + 2\lambda_k\, w_{r_k}}$.
- Update $\omega_{k+1} = \gamma_{\omega_k \to Y_{r_k}}(t_k)$.
**Step 5.** Stop when $d(\omega_{k+1}, \omega_k)$ falls below a chosen tolerance.
For the **cyclic** variant, replace Step 4 with cycling through $i = 1, \ldots, n$ in order, keeping $\lambda_k$ fixed within each full cycle.
The practical consequence is that **no new optimization machinery** is needed for kernel Fréchet regression beyond what was already developed for Fréchet means in Lecture 5. The only change is that the fixed weights $w_i = 1/n$ are replaced by the data-dependent kernel weights $w_{in}^{\mathrm{NW}}(x, h)$, and the algorithm is re-run for each evaluation point $x$ of interest.
::: {.callout-tip title="Efficiency for multiple evaluation points"}
When computing $\hat{\mu}_h^{\mathrm{NW}}(x)$ on a grid of evaluation points $x_1, \ldots, x_m$ (as in the interactive demos below), a natural warm-start strategy is to use the converged estimate at $x_j$ as the initial point $\omega_0$ for $x_{j+1}$. Since $\hat{\mu}_h^{\mathrm{NW}}(x)$ is typically a smooth function of $x$, the warm-started PPA converges in far fewer iterations than a cold start.
:::
## Interactive Exploration: Kernel vs. Global Fréchet Regression on the Sphere
The following demo compares **global Fréchet regression** (Lecture 6) and **kernel (local-constant) Fréchet regression** on the positive quadrant (octant) of the unit sphere $\mathbb{S}^2$. The response space is the octant $\{x, y, z \ge 0\}$ — the same setting as Lecture 6 — but the true regression curve is now **nonlinear**: the spherical angles follow sinusoidal patterns. This makes the global linear fit misspecified, while the kernel fit adapts to local structure.
**Data generation.** $X_i \sim \text{Uniform}(-1, 1)$. The true regression function $\mu(x)$ is a curve on the octant whose spherical angles $\theta(x)$ and $\phi(x)$ are sinusoidal (or constant, depending on the regression type). Each response $Y_i$ is generated by moving a random fraction along the geodesic from a random octant point toward $\mu(X_i)$:
$$Y_i = \gamma_{U_i \to \mu(X_i)}(1 - \sigma\varepsilon_i), \qquad U_i \sim \text{Uniform(octant)},\; \varepsilon_i \sim \text{Uniform}(0,1).$$
**Two estimators compared:**
- **Global Fréchet** (red): $\hat{\mu}_{\text{global}}(x) = \arg\min_\omega \frac{1}{n}\sum_i s_{in}(x)\, d^2(Y_i, \omega)$ with $s_{in}(x) = 1 + \frac{(X_i - \bar{X})(x - \bar{X})}{\hat{\sigma}_X^2}$.
- **Kernel Fréchet** (blue): $\hat{\mu}_h^{\mathrm{NW}}(x) = \arg\min_\omega \sum_i w_{in}^{\mathrm{NW}}(x, h)\, d^2(Y_i, \omega)$ with $w_{in}^{\mathrm{NW}} \propto K_h(X_i - x)$.
Both are computed via **Riemannian gradient descent** on $\mathbb{S}^2$ (the sphere is positively curved, so the PPA of @sec-computation-ppa does not directly apply; an intrinsic gradient scheme with geodesic updates is used instead).
**Visual guide:**
- **Colored dots** = observed responses $Y_i$ (color intensity shows $X_i$)
- **Green dashed curve** = true regression function $\mu(x)$
- **Red curve** = global Fréchet fit $\hat{\mu}_{\text{global}}(x)$
- **Blue curve** = kernel Fréchet fit $\hat{\mu}_h^{\mathrm{NW}}(x)$
- **Gold diamond** = unconditional Fréchet mean of all $Y_i$
```{ojs}
//| label: fig-kernel-frechet-sphere
//| fig-cap: "Interactive: Kernel vs. Global Fréchet regression on the sphere octant"
// ---- Controls ----
n_pts_control = Inputs.range([20, 250], {step: 10, value: 80, label: "Sample size n"})
bandwidth_control = Inputs.range([0.15, 1.2], {step: 0.05, value: 0.4, label: "Bandwidth h"})
reg_type_control = Inputs.select(["Sinusoidal θ & φ", "Varying θ only", "Varying φ only"], {value: "Sinusoidal θ & φ", label: "Regression type"})
noise_level_control = Inputs.range([0.03, 0.4], {step: 0.01, value: 0.12, label: "Noise level σ"})
data_seed_control = Inputs.range([1, 100], {step: 1, value: 42, label: "Random seed"})
eval_x_control = Inputs.range([-1, 1], {step: 0.05, value: 0.5, label: "Evaluation point x"})
n_pts = Generators.input(n_pts_control)
bandwidth = Generators.input(bandwidth_control)
reg_type = Generators.input(reg_type_control)
noise_level = Generators.input(noise_level_control)
data_seed = Generators.input(data_seed_control)
eval_x = Generators.input(eval_x_control)
controls_view = html`
<style>
.kfr-control-panel {
width: 100%;
margin: 0 0 12px;
font-family: system-ui, sans-serif;
font-size: 0.85em;
container-type: inline-size;
}
.kfr-control-panel > *,
.kfr-slider-grid > * {
min-width: 0;
margin: 0;
}
.kfr-slider-grid {
display: flex;
flex-wrap: wrap;
gap: 6px 20px;
margin-top: 6px;
}
.kfr-slider-grid > * {
flex: 1 1 calc((100% - 40px) / 3);
}
.kfr-slider-grid input[type="number"] {
width: 7.5rem !important;
}
@container (max-width: 700px) {
.kfr-slider-grid > * {
flex-basis: calc((100% - 20px) / 2);
}
}
@container (max-width: 480px) {
.kfr-slider-grid > * {
flex-basis: 100%;
}
}
</style>
<div class="kfr-control-panel">
<div>${reg_type_control}</div>
<div class="kfr-slider-grid">
<div>${n_pts_control}</div>
<div>${bandwidth_control}</div>
<div>${noise_level_control}</div>
<div>${data_seed_control}</div>
<div>${eval_x_control}</div>
</div>
</div>
`
// ---- Seeded 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;
};
}
// ---- Sphere geometry ----
function sphereGeodesic(p, q, t) {
var dot = p[0]*q[0] + p[1]*q[1] + p[2]*q[2];
dot = Math.max(-1, Math.min(1, dot));
var theta = Math.acos(dot);
if (theta < 1e-12) {
var rx = (1-t)*p[0] + t*q[0], ry = (1-t)*p[1] + t*q[1], rz = (1-t)*p[2] + t*q[2];
var nr = Math.sqrt(rx*rx + ry*ry + rz*rz);
return [rx/nr, ry/nr, rz/nr];
}
var s = Math.sin(theta);
var w0 = Math.sin((1-t)*theta) / s, w1 = Math.sin(t*theta) / s;
return [w0*p[0] + w1*q[0], w0*p[1] + w1*q[1], w0*p[2] + w1*q[2]];
}
function sphereDist(p, q) {
var dot = p[0]*q[0] + p[1]*q[1] + p[2]*q[2];
return Math.acos(Math.max(-1, Math.min(1, dot)));
}
function sphereLog(p, q) {
var dot = p[0]*q[0] + p[1]*q[1] + p[2]*q[2];
dot = Math.max(-1, Math.min(1, dot));
var theta = Math.acos(dot);
if (theta < 1e-14) return [0, 0, 0];
var scale = theta / Math.sin(theta);
return [scale*(q[0] - dot*p[0]), scale*(q[1] - dot*p[1]), scale*(q[2] - dot*p[2])];
}
function sphereExp(p, v) {
var vNorm = Math.sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
if (vNorm < 1e-14) return [p[0], p[1], p[2]];
var cv = Math.cos(vNorm), sv = Math.sin(vNorm);
return [cv*p[0] + sv/vNorm*v[0], cv*p[1] + sv/vNorm*v[1], cv*p[2] + sv/vNorm*v[2]];
}
// Weighted Fréchet mean on S^2 via Riemannian GD
function weightedSphereFrechetMean(dataPts, weights, maxIter, initLr) {
var n = dataPts.length;
var w = [0, 0, 0], wSum = 0;
for (var i = 0; i < n; i++) {
var aw = Math.abs(weights[i]);
w[0] += aw * dataPts[i][0]; w[1] += aw * dataPts[i][1]; w[2] += aw * dataPts[i][2];
wSum += aw;
}
var norm = Math.sqrt(w[0]*w[0] + w[1]*w[1] + w[2]*w[2]);
w = [w[0]/norm, w[1]/norm, w[2]/norm];
var lr = initLr || 0.3;
for (var iter = 0; iter < maxIter; iter++) {
var grad = [0, 0, 0];
for (var i = 0; i < n; i++) {
var logVal = sphereLog(w, dataPts[i]);
grad[0] -= 2 * weights[i] * logVal[0];
grad[1] -= 2 * weights[i] * logVal[1];
grad[2] -= 2 * weights[i] * logVal[2];
}
grad[0] /= n; grad[1] /= n; grad[2] /= n;
var gNorm = Math.sqrt(grad[0]*grad[0] + grad[1]*grad[1] + grad[2]*grad[2]);
if (gNorm < 1e-10) break;
var step = -lr;
w = sphereExp(w, [step*grad[0], step*grad[1], step*grad[2]]);
lr = initLr / (1 + iter / 80);
}
return w;
}
function weightedSphereFrechetMeanWarm(dataPts, weights, init, maxIter, initLr) {
var n = dataPts.length;
var w = [init[0], init[1], init[2]];
var lr = initLr;
for (var iter = 0; iter < maxIter; iter++) {
var grad = [0, 0, 0];
for (var i = 0; i < n; i++) {
var logVal = sphereLog(w, dataPts[i]);
grad[0] -= 2 * weights[i] * logVal[0];
grad[1] -= 2 * weights[i] * logVal[1];
grad[2] -= 2 * weights[i] * logVal[2];
}
grad[0] /= n; grad[1] /= n; grad[2] /= n;
var gNorm = Math.sqrt(grad[0]*grad[0] + grad[1]*grad[1] + grad[2]*grad[2]);
if (gNorm < 1e-10) break;
var step = -lr;
w = sphereExp(w, [step*grad[0], step*grad[1], step*grad[2]]);
lr = initLr / (1 + iter / 60);
}
return w;
}
// ---- Orthographic projection (matching lecture 6) ----
nview = [1/Math.sqrt(3), 1/Math.sqrt(3), 1/Math.sqrt(3)];
u_scr = [-1/Math.sqrt(2), 1/Math.sqrt(2), 0];
v_scr = [-1/Math.sqrt(6), -1/Math.sqrt(6), 2/Math.sqrt(6)];
function project(p) {
var sx = p[0]*u_scr[0] + p[1]*u_scr[1] + p[2]*u_scr[2];
var sy = p[0]*v_scr[0] + p[1]*v_scr[1] + p[2]*v_scr[2];
return [250 + 210*sx, 200 - 190*sy];
}
// ---- Data generation ----
function randomOctantPoint(rng) {
var z = rng();
var phi = rng() * Math.PI / 2;
var r = Math.sqrt(1 - z*z);
return [r * Math.cos(phi), r * Math.sin(phi), z];
}
// True regression function — NONLINEAR (sinusoidal) unlike lecture 6
function trueRegFunc(x, regType) {
var theta, phi;
if (regType === "Sinusoidal θ & φ") {
theta = Math.PI/4 + 0.35 * Math.sin(Math.PI * x);
phi = Math.PI/4 + 0.35 * Math.cos(Math.PI * x);
} else if (regType === "Varying θ only") {
theta = Math.PI/4 + 0.4 * Math.sin(Math.PI * x);
phi = Math.PI/4;
} else { // Varying φ only
theta = Math.PI/4;
phi = Math.PI/4 + 0.4 * Math.cos(Math.PI * x);
}
theta = Math.max(0.05, Math.min(Math.PI/2 - 0.05, theta));
phi = Math.max(0.05, Math.min(Math.PI/2 - 0.05, phi));
return [Math.sin(theta)*Math.cos(phi), Math.sin(theta)*Math.sin(phi), Math.cos(theta)];
}
function generateSphereRegressionData(n, regType, noise, seed) {
var rng = mulberry32(seed);
var X = [], Y = [];
for (var i = 0; i < n; i++) {
var xi = (rng() - 0.5) * 2;
X.push(xi);
var mu_i = trueRegFunc(xi, regType);
var randPt = randomOctantPoint(rng);
var eps = rng();
var w = 1 - noise * eps;
var yi = sphereGeodesic(randPt, mu_i, w);
Y.push(yi);
}
return { X: X, Y: Y };
}
// ---- Global Fréchet regression (Lecture 6) ----
function computeGlobalFrechetCurve(X, Y, xGrid) {
var n = X.length;
var meanX = 0;
for (var i = 0; i < n; i++) meanX += X[i];
meanX /= n;
var varX = 0;
for (var i = 0; i < n; i++) { var d = X[i] - meanX; varX += d * d; }
varX /= n;
var fitted = [], prevFit = null;
for (var j = 0; j < xGrid.length; j++) {
var x = xGrid[j];
var weights = [];
for (var i = 0; i < n; i++) {
weights.push(1 + (X[i] - meanX) * (x - meanX) / varX);
}
var fit;
if (prevFit) {
fit = weightedSphereFrechetMeanWarm(Y, weights, prevFit, 80, 0.2);
} else {
fit = weightedSphereFrechetMean(Y, weights, 150, 0.25);
}
fitted.push(fit);
prevFit = fit;
}
return fitted;
}
// ---- Kernel Fréchet regression (Lecture 7) ----
function gaussKernel(u) { return Math.exp(-0.5 * u * u); }
function computeKernelFrechetCurve(X, Y, xGrid, h) {
var n = X.length;
var fitted = [], prevFit = null;
for (var j = 0; j < xGrid.length; j++) {
var x = xGrid[j];
var rawW = [], wSum = 0;
for (var i = 0; i < n; i++) {
var kw = gaussKernel((X[i] - x) / h);
rawW.push(kw);
wSum += kw;
}
var weights;
if (wSum > 1e-10) {
weights = rawW.map(function(w) { return w / wSum; });
} else {
weights = rawW.map(function() { return 1/n; });
}
var fit;
if (prevFit) {
fit = weightedSphereFrechetMeanWarm(Y, weights, prevFit, 80, 0.2);
} else {
fit = weightedSphereFrechetMean(Y, weights, 150, 0.25);
}
fitted.push(fit);
prevFit = fit;
}
return fitted;
}
// ---- Run everything ----
regData = generateSphereRegressionData(n_pts, reg_type, noise_level, data_seed);
nGrid = 50;
xGridEval = Array.from({length: nGrid + 1}, function(_, j) { return -1 + j * 2 / nGrid; });
globalCurve = computeGlobalFrechetCurve(regData.X, regData.Y, xGridEval);
kernelCurve = computeKernelFrechetCurve(regData.X, regData.Y, xGridEval, bandwidth);
trueCurve = xGridEval.map(function(x) { return trueRegFunc(x, reg_type); });
// ---- Evaluation at chosen x ----
meanX_val = regData.X.reduce(function(a,b) { return a+b; }, 0) / regData.X.length;
varX_val = regData.X.reduce(function(s, xi) { var d = xi - meanX_val; return s + d*d; }, 0) / regData.X.length;
// Global at eval
gWeightsEval = regData.X.map(function(xi) { return 1 + (xi - meanX_val) * (eval_x - meanX_val) / varX_val; });
globalAtEval = weightedSphereFrechetMean(regData.Y, gWeightsEval, 200, 0.3);
// Kernel at eval
rawEvalW = regData.X.map(function(xi) { return gaussKernel((xi - eval_x) / bandwidth); });
wSumEval = rawEvalW.reduce(function(a,b) { return a+b; }, 0);
nWeightsEval = wSumEval > 1e-10 ? rawEvalW.map(function(w) { return w / wSumEval; }) : rawEvalW.map(function() { return 1/regData.X.length; });
kernelAtEval = weightedSphereFrechetMean(regData.Y, nWeightsEval, 200, 0.3);
trueAtEval = trueRegFunc(eval_x, reg_type);
globalErr = sphereDist(globalAtEval, trueAtEval);
kernelErr = sphereDist(kernelAtEval, trueAtEval);
// Overall Fréchet mean
unifWeights = regData.Y.map(function() { return 1; });
overallMean = weightedSphereFrechetMean(regData.Y, unifWeights, 200, 0.3);
// Kernel diagnostics
effN = rawEvalW.filter(function(w) { return w > 1e-6; }).length;
gNegCount = gWeightsEval.filter(function(w) { return w < 0; }).length;
gWeightSum = gWeightsEval.reduce(function(a,b) { return a+b; }, 0);
// ---- 3D Visualization of the sphere octant ----
function build3DView(regData, trueCurve, globalCurve, kernelCurve, xGridEval, globalAtEval, kernelAtEval, trueAtEval, overallMean, evalX) {
function projectPolyline(pts3d) {
return pts3d.map(function(p) {
var proj = project(p);
return proj[0].toFixed(1) + "," + proj[1].toFixed(1);
}).join(" ");
}
function sampleArc(p, q, nPts) {
var pts = [];
for (var i = 0; i <= nPts; i++) pts.push(sphereGeodesic(p, q, i / nPts));
return pts;
}
var vx = [1,0,0], vy = [0,1,0], vz = [0,0,1];
var arc_xy = sampleArc(vx, vy, 40);
var arc_yz = sampleArc(vy, vz, 40);
var arc_zx = sampleArc(vz, vx, 40);
var boundaryPts = [];
for (var i = 0; i < arc_xy.length; i++) boundaryPts.push(arc_xy[i]);
for (var i = 1; i < arc_yz.length; i++) boundaryPts.push(arc_yz[i]);
for (var i = 1; i < arc_zx.length - 1; i++) boundaryPts.push(arc_zx[i]);
function constantZArc(zVal, nPts) {
var r = Math.sqrt(1 - zVal*zVal);
var pts = [];
for (var i = 0; i <= nPts; i++) {
var alpha = (i / nPts) * Math.PI / 2;
pts.push([r * Math.cos(alpha), r * Math.sin(alpha), zVal]);
}
return pts;
}
var gridLines = [];
[0.2, 0.4, 0.6, 0.8].forEach(function(zv) { gridLines.push(constantZArc(zv, 30)); });
var parts = [];
// Octant surface
parts.push('<path d="M ' + projectPolyline(boundaryPts) + ' Z" fill="#e3f2fd" stroke="none" opacity="0.45"/>');
// Grid lines
gridLines.forEach(function(arc) {
parts.push('<polyline points="' + projectPolyline(arc) + '" fill="none" stroke="#bbdefb" stroke-width="0.5"/>');
});
// Boundary arcs
[arc_xy, arc_yz, arc_zx].forEach(function(arc) {
parts.push('<polyline points="' + projectPolyline(arc) + '" fill="none" stroke="#64b5f6" stroke-width="2.0" stroke-linecap="round"/>');
});
// True regression curve (green, dashed, thicker)
var truePtsStr = trueCurve.map(function(p) { var proj = project(p); return proj[0].toFixed(1) + "," + proj[1].toFixed(1); }).join(" ");
parts.push('<polyline points="' + truePtsStr + '" fill="none" stroke="#2e7d32" stroke-width="3.0" stroke-dasharray="8,4" stroke-linecap="round" stroke-linejoin="round" opacity="0.85"/>');
// Global Fréchet curve (red, solid)
var globalPtsStr = globalCurve.map(function(p) { var proj = project(p); return proj[0].toFixed(1) + "," + proj[1].toFixed(1); }).join(" ");
parts.push('<polyline points="' + globalPtsStr + '" fill="none" stroke="#e53935" stroke-width="2.8" stroke-linecap="round" stroke-linejoin="round" opacity="0.8"/>');
// Kernel Fréchet curve (blue, solid, thicker)
var kernelPtsStr = kernelCurve.map(function(p) { var proj = project(p); return proj[0].toFixed(1) + "," + proj[1].toFixed(1); }).join(" ");
parts.push('<polyline points="' + kernelPtsStr + '" fill="none" stroke="#1971c2" stroke-width="3.0" stroke-linecap="round" stroke-linejoin="round" opacity="0.85"/>');
// Data points (colored by X value)
var xMin = -1, xMax = 1;
regData.Y.forEach(function(yi, idx) {
var xi = regData.X[idx];
var t = (xi - xMin) / (xMax - xMin);
t = Math.max(0, Math.min(1, t));
var r = Math.round(30 + 200 * t);
var g = Math.round(80 + 60 * (1 - Math.abs(t - 0.5) * 2));
var b = Math.round(220 - 200 * t);
var color = "rgb(" + r + "," + g + "," + b + ")";
var proj = project(yi);
parts.push('<circle cx="' + proj[0].toFixed(1) + '" cy="' + proj[1].toFixed(1) +
'" r="3.5" fill="' + color + '" stroke="rgba(0,0,0,0.3)" stroke-width="0.6" opacity="0.8"/>');
parts.push('<title>X=' + xi.toFixed(3) + '</title>');
});
// Evaluation point markers
var gEvalProj = project(globalAtEval);
var kEvalProj = project(kernelAtEval);
var tEvalProj = project(trueAtEval);
// Dashed connectors from both fits to truth
parts.push('<line x1="' + gEvalProj[0].toFixed(1) + '" y1="' + gEvalProj[1].toFixed(1) +
'" x2="' + tEvalProj[0].toFixed(1) + '" y2="' + tEvalProj[1].toFixed(1) +
'" stroke="#795548" stroke-width="1.2" stroke-dasharray="3,2" opacity="0.55"/>');
parts.push('<line x1="' + kEvalProj[0].toFixed(1) + '" y1="' + kEvalProj[1].toFixed(1) +
'" x2="' + tEvalProj[0].toFixed(1) + '" y2="' + tEvalProj[1].toFixed(1) +
'" stroke="#795548" stroke-width="1.2" stroke-dasharray="3,2" opacity="0.55"/>');
// Rings at eval point
parts.push('<circle cx="' + gEvalProj[0].toFixed(1) + '" cy="' + gEvalProj[1].toFixed(1) +
'" r="7" fill="none" stroke="#e53935" stroke-width="3"/>');
parts.push('<circle cx="' + kEvalProj[0].toFixed(1) + '" cy="' + kEvalProj[1].toFixed(1) +
'" r="7" fill="none" stroke="#1971c2" stroke-width="3"/>');
parts.push('<circle cx="' + tEvalProj[0].toFixed(1) + '" cy="' + tEvalProj[1].toFixed(1) +
'" r="7" fill="none" stroke="#2e7d32" stroke-width="3" stroke-dasharray="5,3"/>');
// Overall Fréchet mean (gold diamond)
var omProj = project(overallMean);
parts.push('<polygon points="' +
(omProj[0]).toFixed(1) + ',' + (omProj[1]-7).toFixed(1) + ' ' +
(omProj[0]+5.5).toFixed(1) + ',' + omProj[1].toFixed(1) + ' ' +
omProj[0].toFixed(1) + ',' + (omProj[1]+7).toFixed(1) + ' ' +
(omProj[0]-5.5).toFixed(1) + ',' + omProj[1].toFixed(1) +
'" fill="#f9a825" stroke="#f57f17" stroke-width="1.5"/>');
// Axis labels
var lvx = project([1.15, 0, 0]), lvy = project([0, 1.15, 0]), lvz = project([0, 0, 1.15]);
parts.push('<text x="' + lvx[0].toFixed(1) + '" y="' + (lvx[1]+4).toFixed(1) +
'" text-anchor="middle" font-size="13" fill="#546e7a" font-style="italic">x</text>');
parts.push('<text x="' + lvy[0].toFixed(1) + '" y="' + (lvy[1]+4).toFixed(1) +
'" text-anchor="middle" font-size="13" fill="#546e7a" font-style="italic">y</text>');
parts.push('<text x="' + lvz[0].toFixed(1) + '" y="' + (lvz[1]-6).toFixed(1) +
'" text-anchor="middle" font-size="13" fill="#546e7a" font-style="italic">z</text>');
// Legend
var lx = 370, ly0 = 40;
parts.push('<circle cx="' + lx + '" cy="' + ly0 + '" r="4" fill="rgb(160,100,100)"/>');
parts.push('<text x="' + (lx+12) + '" y="' + (ly0+4) + '" font-size="9" fill="#37474f">Data Yᵢ (color = Xᵢ)</text>');
parts.push('<line x1="' + (lx-6) + '" y1="' + (ly0+17) + '" x2="' + (lx+14) + '" y2="' + (ly0+17) +
'" stroke="#2e7d32" stroke-width="2.5" stroke-dasharray="5,3"/>');
parts.push('<text x="' + (lx+22) + '" y="' + (ly0+21) + '" font-size="9" fill="#37474f">True μ(x)</text>');
parts.push('<line x1="' + (lx-6) + '" y1="' + (ly0+34) + '" x2="' + (lx+14) + '" y2="' + (ly0+34) +
'" stroke="#e53935" stroke-width="2.5"/>');
parts.push('<text x="' + (lx+22) + '" y="' + (ly0+38) + '" font-size="9" fill="#37474f">Global μ̂(x)</text>');
parts.push('<line x1="' + (lx-6) + '" y1="' + (ly0+51) + '" x2="' + (lx+14) + '" y2="' + (ly0+51) +
'" stroke="#1971c2" stroke-width="2.5"/>');
parts.push('<text x="' + (lx+22) + '" y="' + (ly0+55) + '" font-size="9" fill="#37474f">Kernel μ̂(x)</text>');
parts.push('<polygon points="' + lx + ',' + (ly0+65) + ' ' + (lx+4) + ',' + (ly0+70) +
' ' + lx + ',' + (ly0+75) + ' ' + (lx-4) + ',' + (ly0+70) + '" fill="#f9a825"/>');
parts.push('<text x="' + (lx+12) + '" y="' + (ly0+74) + '" font-size="9" fill="#37474f">Fréchet mean</text>');
return '<svg viewBox="0 0 520 420" xmlns="http://www.w3.org/2000/svg" style="width:100%;max-width:540px;">\n' +
parts.join('\n') + '\n</svg>';
}
sphere_view = html`<div style="text-align:center;">${build3DView(regData, trueCurve, globalCurve, kernelCurve, xGridEval, globalAtEval, kernelAtEval, trueAtEval, overallMean, eval_x)}</div>`
// ---- 2D Unwrapped View: θ and φ vs X ----
function sphericalAngles(p) {
var theta = Math.acos(Math.max(-1, Math.min(1, p[2])));
var phi = Math.atan2(p[1], p[0]);
if (phi < 0) phi += 2 * Math.PI;
return { theta: theta, phi: phi };
}
trueAngles = trueCurve.map(function(p) { return sphericalAngles(p); });
globalAngles = globalCurve.map(function(p) { return sphericalAngles(p); });
kernelAngles = kernelCurve.map(function(p) { return sphericalAngles(p); });
dataAngles = regData.Y.map(function(yi) { return sphericalAngles(yi); });
evalTrueAng = sphericalAngles(trueAtEval);
evalGlobalAng = sphericalAngles(globalAtEval);
evalKernelAng = sphericalAngles(kernelAtEval);
function buildUnwrappedPlot(title, yLabel, trueVals, globalVals, kernelVals, dataVals, xVals, xGrid, yMin, yMax) {
var W = 520, H = 200;
var margin = {top: 18, right: 15, bottom: 32, left: 42};
var pw = W - margin.left - margin.right;
var ph = H - margin.top - margin.bottom;
function xToSvg(x) { return margin.left + (x - (-1)) / 2 * pw; }
function yToSvg(y) { return margin.top + ph - (y - yMin) / (yMax - yMin) * ph; }
var parts = [];
parts.push('<rect x="' + margin.left + '" y="' + margin.top + '" width="' + pw +
'" height="' + ph + '" fill="#fafafa" rx="2"/>');
// Grid lines
for (var gy = 0; gy <= 4; gy++) {
var yv = yMin + (gy/4) * (yMax - yMin);
var yp = yToSvg(yv);
parts.push('<line x1="' + margin.left + '" y1="' + yp + '" x2="' + (margin.left+pw) +
'" y2="' + yp + '" stroke="#e0e0e0" stroke-width="0.5"/>');
parts.push('<text x="' + (margin.left-4) + '" y="' + (yp+4) +
'" text-anchor="end" font-size="8" fill="#9e9e9e">' + yv.toFixed(2) + '</text>');
}
// Data points
dataVals.forEach(function(dv, idx) {
var sx = xToSvg(xVals[idx]), sy = yToSvg(dv);
parts.push('<circle cx="' + sx + '" cy="' + sy + '" r="2.2" fill="#90caf9" opacity="0.45"/>');
});
// True curve (green, dashed)
var trueLine = "";
for (var j = 0; j < xGrid.length; j++) {
trueLine += (j === 0 ? 'M' : 'L') + " " + xToSvg(xGrid[j]).toFixed(1) + " " + yToSvg(trueVals[j]).toFixed(1);
}
parts.push('<path d="' + trueLine + '" fill="none" stroke="#2e7d32" stroke-width="2.2" stroke-dasharray="7,4"/>');
// Global fit (red)
var globalLine = "";
for (var j = 0; j < xGrid.length; j++) {
globalLine += (j === 0 ? 'M' : 'L') + " " + xToSvg(xGrid[j]).toFixed(1) + " " + yToSvg(globalVals[j]).toFixed(1);
}
parts.push('<path d="' + globalLine + '" fill="none" stroke="#e53935" stroke-width="2.0"/>');
// Kernel fit (blue)
var kernelLine = "";
for (var j = 0; j < xGrid.length; j++) {
kernelLine += (j === 0 ? 'M' : 'L') + " " + xToSvg(xGrid[j]).toFixed(1) + " " + yToSvg(kernelVals[j]).toFixed(1);
}
parts.push('<path d="' + kernelLine + '" fill="none" stroke="#1971c2" stroke-width="2.2"/>');
// Evaluation point line
var exSvg = xToSvg(eval_x);
parts.push('<line x1="' + exSvg + '" y1="' + margin.top + '" x2="' + exSvg +
'" y2="' + (margin.top+ph) + '" stroke="#795548" stroke-width="1.2" stroke-dasharray="3,3" opacity="0.5"/>');
// Axes
parts.push('<line x1="' + margin.left + '" y1="' + (margin.top+ph) + '" x2="' +
(margin.left+pw) + '" y2="' + (margin.top+ph) + '" stroke="#424242" stroke-width="1"/>');
parts.push('<line x1="' + margin.left + '" y1="' + margin.top + '" x2="' +
margin.left + '" y2="' + (margin.top+ph) + '" stroke="#424242" stroke-width="1"/>');
// Axis labels
parts.push('<text x="' + (margin.left+pw/2) + '" y="' + (H-4) +
'" text-anchor="middle" font-size="10" fill="#616161">Predictor X</text>');
parts.push('<text x="' + (margin.left-34) + '" y="' + (margin.top+ph/2) +
'" text-anchor="middle" font-size="10" fill="#616161" ' +
'transform="rotate(-90,' + (margin.left-34) + ',' + (margin.top+ph/2) + ')">' + yLabel + '</text>');
// Title
parts.push('<text x="' + (margin.left+pw/2) + '" y="' + (margin.top-2) +
'" text-anchor="middle" font-size="11" font-weight="bold" fill="#37474f">' + title + '</text>');
// Legend
var lx = margin.left + pw - 155, ly = margin.top + 4;
parts.push('<line x1="' + lx + '" y1="' + ly + '" x2="' + (lx+14) + '" y2="' + ly +
'" stroke="#2e7d32" stroke-width="1.8" stroke-dasharray="4,3"/>');
parts.push('<text x="' + (lx+18) + '" y="' + (ly+4) + '" font-size="8" fill="#616161">True</text>');
parts.push('<line x1="' + (lx+48) + '" y1="' + ly + '" x2="' + (lx+62) + '" y2="' + ly +
'" stroke="#e53935" stroke-width="1.8"/>');
parts.push('<text x="' + (lx+66) + '" y="' + (ly+4) + '" font-size="8" fill="#616161">Global</text>');
parts.push('<line x1="' + (lx+100) + '" y1="' + ly + '" x2="' + (lx+114) + '" y2="' + ly +
'" stroke="#1971c2" stroke-width="2"/>');
parts.push('<text x="' + (lx+118) + '" y="' + (ly+4) + '" font-size="8" fill="#616161">Kernel</text>');
parts.push('<rect x="' + margin.left + '" y="' + margin.top + '" width="' + pw +
'" height="' + ph + '" fill="none" stroke="#dee2e6"/>');
return '<svg viewBox="0 0 ' + W + ' ' + H + '" xmlns="http://www.w3.org/2000/svg" style="width:100%;max-width:540px;">\n' +
parts.join('\n') + '\n</svg>';
}
unwrapped_view = html`
<div style="max-width:540px; margin:8px auto; font-family:system-ui,sans-serif;">
${buildUnwrappedPlot("Polar angle θ vs X", "θ (rad)",
trueAngles.map(function(a) { return a.theta; }),
globalAngles.map(function(a) { return a.theta; }),
kernelAngles.map(function(a) { return a.theta; }),
dataAngles.map(function(a) { return a.theta; }),
regData.X, xGridEval, 0, Math.PI/2)}
<div style="margin-top:10px;">
${buildUnwrappedPlot("Azimuthal angle φ vs X", "φ (rad)",
trueAngles.map(function(a) { return a.phi; }),
globalAngles.map(function(a) { return a.phi; }),
kernelAngles.map(function(a) { return a.phi; }),
dataAngles.map(function(a) { return a.phi; }),
regData.X, xGridEval, 0, Math.PI/2)}
</div>
</div>
`
// ---- Statistics panel ----
statistics_view = html`
<div style="font-family:system-ui,sans-serif; max-width:540px; margin:10px auto; font-size:0.88em;">
<div style="background:#f8f9fa; border-radius:6px; padding:12px 16px;">
<table style="width:100%; border-collapse:collapse;">
<tr style="border-bottom:1px solid #dee2e6;">
<td style="padding:5px 8px;" colspan="2"><b>Evaluation at x = ${eval_x.toFixed(2)}</b> (X̄ = ${meanX_val.toFixed(3)})</td>
</tr>
<tr style="border-bottom:1px solid #eee;">
<td style="padding:4px 8px;" colspan="2"><b>Kernel Fréchet</b> (h = ${bandwidth.toFixed(2)})</td>
</tr>
<tr>
<td style="padding:4px 8px;">Kernel μ̂(x)</td>
<td style="padding:4px 8px; text-align:right; color:#1971c2; font-weight:bold;">θ=${evalKernelAng.theta.toFixed(3)}, φ=${evalKernelAng.phi.toFixed(3)} rad</td>
</tr>
<tr>
<td style="padding:4px 8px;">Kernel error</td>
<td style="padding:4px 8px; text-align:right; color:#1971c2; font-weight:bold;">${kernelErr.toFixed(5)} rad (${(kernelErr*180/Math.PI).toFixed(3)}°)</td>
</tr>
<tr>
<td style="padding:4px 8px;">Effective local n</td>
<td style="padding:4px 8px; text-align:right;">~${effN} of ${regData.X.length}</td>
</tr>
<tr style="border-bottom:1px solid #eee;">
<td style="padding:4px 8px;" colspan="2"><b>Global Fréchet</b></td>
</tr>
<tr>
<td style="padding:4px 8px;">Global μ̂(x)</td>
<td style="padding:4px 8px; text-align:right; color:#e53935;">θ=${evalGlobalAng.theta.toFixed(3)}, φ=${evalGlobalAng.phi.toFixed(3)} rad</td>
</tr>
<tr>
<td style="padding:4px 8px;">Global error</td>
<td style="padding:4px 8px; text-align:right; color:#e53935; font-weight:bold;">${globalErr.toFixed(5)} rad (${(globalErr*180/Math.PI).toFixed(3)}°)</td>
</tr>
<tr>
<td style="padding:4px 8px;">Negative weights</td>
<td style="padding:4px 8px; text-align:right;">${gNegCount} of ${regData.X.length} (${(gNegCount/regData.X.length*100).toFixed(1)}%)</td>
</tr>
<tr style="border-bottom:1px solid #eee;">
<td style="padding:4px 8px;" colspan="2"><b>True value</b></td>
</tr>
<tr>
<td style="padding:4px 8px;">True μ(x)</td>
<td style="padding:4px 8px; text-align:right;">θ=${evalTrueAng.theta.toFixed(3)}, φ=${evalTrueAng.phi.toFixed(3)} rad</td>
</tr>
<tr>
<td style="padding:4px 8px;">Global Fréchet mean μ̂(X̄)</td>
<td style="padding:4px 8px; text-align:right;">(${overallMean[0].toFixed(3)}, ${overallMean[1].toFixed(3)}, ${overallMean[2].toFixed(3)})</td>
</tr>
</table>
</div>
</div>
`
html`<div>${controls_view}${sphere_view}${unwrapped_view}${statistics_view}</div>`
```
::: {.callout-tip title="Try these experiments"}
- **Toggle regression types:** "Sinusoidal θ & φ" makes both angles oscillate — the global fit (a straight-line analogue on the sphere) cannot capture the curvature. The kernel fit tracks it closely.
- **Vary the bandwidth $h$:** Small $h$ (0.15–0.25) gives a wiggly kernel fit (high variance, low bias); large $h$ (0.8–1.2) shrinks the kernel fit toward the unconditional Fréchet mean (gold diamond), approaching the global fit's behavior at $\bar{X}$.
- **Increase noise $\sigma$ (0.3–0.4):** Both fits shrink toward the Fréchet mean, but the kernel fit adapts locally — it shrinks more where data are sparse.
- **Evaluate near edges ($x = \pm 1$):** The global fit uses negative weights (reported in the stats panel). The kernel weights remain nonnegative but the effective local $n$ drops — watch the error grow.
- **Increase $n$ (200+):** The kernel fit converges to the true curve everywhere; the global fit remains biased in the sinusoidal regions no matter how large $n$ is.
- **Study the 2D unwrapped views:** The kernel fit tracks the sinusoidal $\theta(x)$ and $\phi(x)$ while the global fit approximates them with straight lines — this is exactly the bias–variance tradeoff.
:::
## Application: Capturing Nonlinear Risk Dynamics with Kernel Smoothing
Lecture 6 modeled how the typical market covariance structure varies with a continuous stress measure using **global Fréchet regression** — the metric analogue of a straight-line fit. But what if the relationship between market stress and risk structure is **nonlinear**? For instance, covariance matrices might change slowly at low stress levels, then shift abruptly as stress crosses a critical threshold, and finally saturate. A global linear model would smooth over these features; kernel Fréchet regression can capture them.
### Why Kernel? Nonlinear Stress–Covariance Relationships
In the market risk application (Lecture 4: Fréchet ANOVA for market risk regimes; Lecture 6: global Fréchet regression for stress-dependent covariance), we used the same six-asset universe (SPY, TLT, GLD, XLF, XLE, QQQ) with rolling covariance matrices $\Sigma_t \in \mathcal{S}_{++}^p$ as responses and realized SPY volatility $\hat{\sigma}_t$ as the predictor. Under global Fréchet regression, the fitted covariance at stress level $x$ is
$$
\log \hat{\mu}_{\text{global}}(x) = \frac{1}{n}\sum_{i=1}^n s_{in}(x)\, \log \Sigma_i,
\qquad
s_{in}(x) = 1 + \frac{(X_i - \bar{X})(x - \bar{X})}{\hat{\sigma}_X^2}.
$$
The weights $s_{in}(x)$ are **linear in $x$** — the model assumes that each $\log \Sigma_i$ contributes to the fit with an influence that grows linearly with the predictor value. This is the key limitation: if the true relationship has curvature, thresholds, or saturation, a linear weight function cannot adapt.
Kernel Fréchet regression replaces $s_{in}(x)$ with **localized** weights:
$$
\log \hat{\mu}_h^{\mathrm{NW}}(x) = \sum_{i=1}^n w_{in}^{\mathrm{NW}}(x, h)\, \log \Sigma_i,
\qquad
w_{in}^{\mathrm{NW}}(x, h) = \frac{K_h(X_i - x)}{\sum_{j=1}^n K_h(X_j - x)}.
$$
These weights concentrate around observations with $X_i \approx x$ and decay smoothly with distance — no global linearity assumption is imposed.
::: {.callout-important title="From linear to local: the same computational primitive"}
Under the log-Euclidean metric, **all three methods** (global, kernel, local-linear) share the same computational core: a weighted average of matrix logarithms followed by matrix exponentiation. The only difference is **how the weights are computed**. This means switching from global to kernel regression requires no new geometry code — only a different weight function.
:::
### What the Kernel Fit Reveals
The kernel Fréchet regression curve $\hat{\mu}_h^{\mathrm{NW}}(x)$ traces a **data-adaptive** path through the SPD manifold:
- **Near stress level $x$**, only rolling windows with similar volatility matter — the fit reflects the *local* average covariance structure.
- **In flat regions** of the relationship, the kernel fit behaves like the global fit (both approximate a local constant).
- **Near transitions or thresholds**, the kernel fit can bend to follow the data, while the global fit is constrained to a straight-line path in the log domain.
- **At the edges** of the predictor range, the kernel fit may exhibit boundary bias (see Lecture 8 for the fix).
The difference between the global and kernel fits at a given stress level is a **diagnostic for nonlinearity**: large discrepancies suggest that the linear model is misspecified and that local structure matters.
### Portfolio-Risk Interpretation
For a portfolio weight vector $w$, the stress-dependent variance estimate is
$$
\sigma_P^2(x) = w^\top \hat{\mu}_h^{\mathrm{NW}}(x) w.
$$
If the kernel fit differs substantially from the global fit at certain stress levels, then a risk manager using the global model would **misestimate** portfolio variance at those levels. Kernel regression provides a more honest uncertainty assessment when the true stress–risk relationship is unknown.
### Interactive Exploration: Kernel vs. Global on Simulated SPD Data
The following demo illustrates the difference between global and kernel Fréchet regression on simulated $2 \times 2$ SPD matrices with a **nonlinear** true relationship. We work directly in the log domain (symmetric matrices) so that the weighted Fréchet mean is simply a weighted average.
**Setup:** $X_i \sim \text{Uniform}(0, 1)$. The true log-matrix $\log \mu(x)$ has entries that are nonlinear functions of $x$ (sinusoidal). Noisy log-matrices are generated by adding Gaussian noise. Both global and kernel regression compute weighted averages of the noisy log-matrices; the fitted SPD matrices are obtained by matrix exponentiation.
**Visual guide:**
- **Gray dots** = noisy log-matrix entries of the observed data
- **Green dashed curve** = true function (sinusoidal)
- **Red curve** = global Fréchet fit (straight-line analogue)
- **Blue curve** = kernel Fréchet fit
```{ojs}
//| label: fig-kernel-app-spd
//| fig-cap: "Interactive: Kernel vs. Global Fréchet regression on simulated 2×2 SPD matrices"
// ---- Controls ----
n_app_control = Inputs.range([20, 200], {step: 10, value: 80, label: "Sample size n"})
h_app_control = Inputs.range([0.05, 0.6], {step: 0.02, value: 0.15, label: "Bandwidth h"})
noise_app_control = Inputs.range([0.05, 0.5], {step: 0.02, value: 0.2, label: "Noise level σ"})
seed_app_control = Inputs.range([1, 100], {step: 1, value: 42, label: "Random seed"})
n_app = Generators.input(n_app_control)
h_app = Generators.input(h_app_control)
noise_app = Generators.input(noise_app_control)
seed_app = Generators.input(seed_app_control)
controls_app_view = html`
<style>
.spd-kfr-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;
}
.spd-kfr-slider-grid > * {
min-width: 0;
margin: 0;
}
@media (max-width: 500px) {
.spd-kfr-slider-grid {
grid-template-columns: 1fr;
}
}
</style>
<div class="spd-kfr-slider-grid">
<div>${n_app_control}</div>
<div>${h_app_control}</div>
<div>${noise_app_control}</div>
<div>${seed_app_control}</div>
</div>
`
// ---- 2×2 symmetric matrix utilities ----
// Matrix represented as [a, b, c] for [[a, b], [b, c]]
function symTrace(m) { return m[0] + m[2]; }
function symDet(m) { return m[0]*m[2] - m[1]*m[1]; }
// Eigenvalues of 2×2 symmetric matrix
function symEigenvalues(m) {
var tr = symTrace(m), det = symDet(m);
var disc = Math.sqrt(Math.max(0, tr*tr - 4*det));
return [(tr - disc)/2, (tr + disc)/2];
}
// Matrix exponential of 2×2 symmetric matrix (log-domain → SPD)
function symExp(m) {
var evals = symEigenvalues(m);
var e1 = Math.exp(evals[0]), e2 = Math.exp(evals[1]);
var a = m[0], b = m[1], c = m[2];
// Eigenvector for eval[0]: [b, evals[0]-a] normalized
var vx = b, vy = evals[0] - a;
var vn = Math.sqrt(vx*vx + vy*vy);
if (vn < 1e-14) {
// Diagonal case: m is already diagonal
return [Math.exp(a), 0, Math.exp(c)];
}
vx /= vn; vy /= vn;
// Reconstruct: V * diag(e1,e2) * V^T
return [
e1*vx*vx + e2*vy*vy,
e1*vx*(-vy) + e2*vy*vx, // second eigenvector is [-vy, vx]
e1*vy*vy + e2*vx*vx
];
}
// ---- Data generation ----
// True log-matrix entries as nonlinear functions of x
function trueLogMatrix(x) {
var a = 0.5 + 1.5*x + 0.6*Math.sin(2.5*Math.PI*x); // diagonal entry 1
var b = -0.3 + 0.8*x + 0.4*Math.cos(2.0*Math.PI*x); // off-diagonal
var c = 0.8 + 1.0*x + 0.5*Math.sin(1.8*Math.PI*x + 0.5); // diagonal entry 2
return [a, b, c];
}
function generateSPDData(n, noise, seed) {
var rng = mulberry32(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);
}
var X = [], logs = [];
for (var i = 0; i < n; i++) {
var xi = rng(); // Uniform(0, 1)
X.push(xi);
var truth = trueLogMatrix(xi);
logs.push([truth[0] + noise*randn(), truth[1] + noise*randn(), truth[2] + noise*randn()]);
}
return { X: X, logs: logs };
}
// Weighted average of symmetric matrices
function symWeightedMean(mats, weights) {
var n = mats.length;
var wSum = 0, a = 0, b = 0, c = 0;
for (var i = 0; i < n; i++) {
a += weights[i] * mats[i][0];
b += weights[i] * mats[i][1];
c += weights[i] * mats[i][2];
wSum += weights[i];
}
return [a/wSum, b/wSum, c/wSum];
}
// ---- Generate data ----
appData = generateSPDData(n_app, noise_app, seed_app);
// ---- Evaluation grid ----
nGrid_app = 60;
xGrid_app = Array.from({length: nGrid_app + 1}, function(_, j) { return j / nGrid_app; });
// True curve
trueLogs_app = xGrid_app.map(function(x) { return trueLogMatrix(x); });
trueTraces = trueLogs_app.map(symTrace);
trueDets = trueLogs_app.map(symDet);
// Global Fréchet regression
meanX_app = appData.X.reduce(function(a,b) { return a+b; }, 0) / n_app;
varX_app = appData.X.reduce(function(s, xi) { var d = xi - meanX_app; return s + d*d; }, 0) / n_app;
globalLogs = xGrid_app.map(function(x) {
var weights = appData.X.map(function(xi) { return 1 + (xi - meanX_app)*(x - meanX_app)/varX_app; });
return symWeightedMean(appData.logs, weights);
});
globalTraces = globalLogs.map(symTrace);
globalDets = globalLogs.map(symDet);
// Kernel Fréchet regression
kernelLogs = xGrid_app.map(function(x) {
var rawW = appData.X.map(function(xi) { return gaussKernel((xi - x)/h_app); });
var wSum = rawW.reduce(function(a,b) { return a+b; }, 0);
var weights = wSum > 1e-10 ? rawW.map(function(w) { return w/wSum; }) : rawW.map(function() { return 1/n_app; });
return symWeightedMean(appData.logs, weights);
});
kernelTraces = kernelLogs.map(symTrace);
kernelDets = kernelLogs.map(symDet);
// Data scalar summaries
dataTraces = appData.logs.map(symTrace);
dataDets = appData.logs.map(symDet);
// Global Fréchet mean (unconditional, at mean X)
unifW = appData.logs.map(function() { return 1; });
meanLog = symWeightedMean(appData.logs, unifW);
meanSPD = symExp(meanLog);
// ---- Visualization ----
function buildAppPlot(title, yLabel, trueVals, globalVals, kernelVals, dataVals, xVals, xGrid, yMin, yMax) {
var W = 500, H = 200;
var margin = {top: 18, right: 15, bottom: 32, left: 45};
var pw = W - margin.left - margin.right;
var ph = H - margin.top - margin.bottom;
function xS(x) { return margin.left + x * pw; }
function yS(y) { return margin.top + ph - (y - yMin)/(yMax - yMin)*ph; }
var parts = [];
parts.push('<rect x="' + margin.left + '" y="' + margin.top + '" width="' + pw +
'" height="' + ph + '" fill="#fafafa" rx="2"/>');
// Grid lines
for (var gy = 0; gy <= 4; gy++) {
var yv = yMin + (gy/4)*(yMax - yMin);
var yp = yS(yv);
parts.push('<line x1="' + margin.left + '" y1="' + yp + '" x2="' + (margin.left+pw) +
'" y2="' + yp + '" stroke="#e0e0e0" stroke-width="0.5"/>');
parts.push('<text x="' + (margin.left-4) + '" y="' + (yp+4) +
'" text-anchor="end" font-size="8" fill="#9e9e9e">' + yv.toFixed(2) + '</text>');
}
// Data points
dataVals.forEach(function(dv, idx) {
parts.push('<circle cx="' + xS(xVals[idx]) + '" cy="' + yS(dv) +
'" r="2.5" fill="#adb5bd" opacity="0.5"/>');
});
// True curve (green dashed)
var d = "";
for (var j = 0; j < xGrid.length; j++) d += (j===0?'M':'L') + " " + xS(xGrid[j]).toFixed(1) + " " + yS(trueVals[j]).toFixed(1);
parts.push('<path d="' + d + '" fill="none" stroke="#2e7d32" stroke-width="2.2" stroke-dasharray="7,4"/>');
// Global fit (red)
d = "";
for (var j = 0; j < xGrid.length; j++) d += (j===0?'M':'L') + " " + xS(xGrid[j]).toFixed(1) + " " + yS(globalVals[j]).toFixed(1);
parts.push('<path d="' + d + '" fill="none" stroke="#e53935" stroke-width="2.0"/>');
// Kernel fit (blue)
d = "";
for (var j = 0; j < xGrid.length; j++) d += (j===0?'M':'L') + " " + xS(xGrid[j]).toFixed(1) + " " + yS(kernelVals[j]).toFixed(1);
parts.push('<path d="' + d + '" fill="none" stroke="#1971c2" stroke-width="2.4"/>');
// Axes
parts.push('<line x1="' + margin.left + '" y1="' + (margin.top+ph) + '" x2="' + (margin.left+pw) +
'" y2="' + (margin.top+ph) + '" stroke="#424242" stroke-width="1"/>');
parts.push('<line x1="' + margin.left + '" y1="' + margin.top + '" x2="' +
margin.left + '" y2="' + (margin.top+ph) + '" stroke="#424242" stroke-width="1"/>');
parts.push('<text x="' + (margin.left+pw/2) + '" y="' + (H-4) +
'" text-anchor="middle" font-size="10" fill="#616161">Stress level X</text>');
parts.push('<text x="' + (margin.left-38) + '" y="' + (margin.top+ph/2) +
'" text-anchor="middle" font-size="10" fill="#616161" ' +
'transform="rotate(-90,' + (margin.left-38) + ',' + (margin.top+ph/2) + ')">' + yLabel + '</text>');
parts.push('<text x="' + (margin.left+pw/2) + '" y="' + (margin.top-2) +
'" text-anchor="middle" font-size="11" font-weight="bold" fill="#37474f">' + title + '</text>');
// Legend
var lx = margin.left + pw - 140, ly = margin.top + 3;
parts.push('<line x1="' + lx + '" y1="' + ly + '" x2="' + (lx+14) + '" y2="' + ly +
'" stroke="#2e7d32" stroke-width="1.8" stroke-dasharray="4,3"/>');
parts.push('<text x="' + (lx+18) + '" y="' + (ly+4) + '" font-size="8" fill="#616161">True</text>');
parts.push('<line x1="' + (lx+46) + '" y1="' + ly + '" x2="' + (lx+60) + '" y2="' + ly +
'" stroke="#e53935" stroke-width="1.8"/>');
parts.push('<text x="' + (lx+64) + '" y="' + (ly+4) + '" font-size="8" fill="#616161">Global</text>');
parts.push('<line x1="' + (lx+96) + '" y1="' + ly + '" x2="' + (lx+110) + '" y2="' + ly +
'" stroke="#1971c2" stroke-width="2"/>');
parts.push('<text x="' + (lx+114) + '" y="' + (ly+4) + '" font-size="8" fill="#616161">Kernel</text>');
parts.push('<rect x="' + margin.left + '" y="' + margin.top + '" width="' + pw +
'" height="' + ph + '" fill="none" stroke="#dee2e6"/>');
return '<svg viewBox="0 0 ' + W + ' ' + H + '" xmlns="http://www.w3.org/2000/svg" style="width:100%;max-width:520px;">\n' +
parts.join('\n') + '\n</svg>';
}
// Compute y ranges
traceMin = Math.min.apply(null, dataTraces.concat(trueTraces, globalTraces, kernelTraces));
traceMax = Math.max.apply(null, dataTraces.concat(trueTraces, globalTraces, kernelTraces));
tracePad = (traceMax - traceMin) * 0.12;
traceMinAdj = traceMin - tracePad; traceMaxAdj = traceMax + tracePad;
detMin = Math.min.apply(null, dataDets.concat(trueDets, globalDets, kernelDets));
detMax = Math.max.apply(null, dataDets.concat(trueDets, globalDets, kernelDets));
detPad = (detMax - detMin) * 0.12;
detMinAdj = detMin - detPad; detMaxAdj = detMax + detPad;
plot_app_view = html`
<div style="font-family:system-ui,sans-serif; max-width:540px; margin:0 auto;">
${buildAppPlot("Trace of log Σ(x) vs. Stress", "tr(log Σ)", trueTraces, globalTraces, kernelTraces, dataTraces, appData.X, xGrid_app, traceMinAdj, traceMaxAdj)}
<div style="margin-top:10px;">
${buildAppPlot("Determinant of log Σ(x) vs. Stress", "det(log Σ)", trueDets, globalDets, kernelDets, dataDets, appData.X, xGrid_app, detMinAdj, detMaxAdj)}
</div>
<p style="margin-top:12px; font-size:0.85em; color:#495057;">
<b>Observe:</b> The global fit (red) approximates the sinusoidal truth with a straight line — fundamentally missing the oscillation.
The kernel fit (blue, h = ${h_app.toFixed(2)}) adapts to local structure. The trace and determinant are scalar summaries of the
2×2 log-covariance matrix; the actual regression operates on the full matrix, but these summaries reveal where the global
model is most severely misspecified.
</p>
<div style="background:#f8f9fa; border-radius:6px; padding:10px 14px; margin-top:8px; font-size:0.85em;">
<b>Fitted SPD matrices at X = 0.5 (mid-stress):</b>
<table style="width:100%; border-collapse:collapse; margin-top:6px;">
<tr><td style="padding:2px 8px;">Global μ̂(0.5)</td>
<td style="padding:2px 8px; text-align:right; font-family:monospace;">
[[${(Math.exp(globalLogs[30][0])).toFixed(2)}, ${(Math.exp(globalLogs[30][1])).toFixed(2)}],
[${(Math.exp(globalLogs[30][1])).toFixed(2)}, ${(Math.exp(globalLogs[30][2])).toFixed(2)}]]</td></tr>
<tr><td style="padding:2px 8px;">Kernel μ̂(0.5)</td>
<td style="padding:2px 8px; text-align:right; font-family:monospace;">
[[${(Math.exp(kernelLogs[30][0])).toFixed(2)}, ${(Math.exp(kernelLogs[30][1])).toFixed(2)}],
[${(Math.exp(kernelLogs[30][1])).toFixed(2)}, ${(Math.exp(kernelLogs[30][2])).toFixed(2)}]]</td></tr>
<tr><td style="padding:2px 8px;">Unconditional Fréchet mean</td>
<td style="padding:2px 8px; text-align:right; font-family:monospace;">
[[${meanSPD[0].toFixed(2)}, ${meanSPD[1].toFixed(2)}],
[${meanSPD[1].toFixed(2)}, ${meanSPD[2].toFixed(2)}]]</td></tr>
</table>
</div>
</div>
`
html`<div>${controls_app_view}${plot_app_view}</div>`
```
::: {.callout-tip title="Try these experiments"}
- **Decrease $h$ (0.05–0.10):** The kernel fit becomes wiggly, chasing individual data points — high variance, low bias.
- **Increase $h$ (0.4–0.6):** The kernel fit flattens toward the unconditional Fréchet mean — low variance, high bias. At very large $h$, it nearly coincides with the global fit at $\bar{X}$.
- **Increase noise $\sigma$ (0.35–0.5):** Both fits become less certain. The kernel fit's local adaptivity means it can still track the general shape, but the global fit's bias dominates.
- **Compare trace and determinant:** Mismatches between the global and kernel fits in different scalar summaries reveal which aspects of the covariance structure are most nonlinear in stress.
:::
### Bandwidth Selection in Practice
For SPD-valued kernel regression, bandwidth selection follows the same principles as in the Euclidean case but must respect the geometry:
- **Cross-validation on the manifold:** Define a prediction error using the metric distance, $\mathrm{CV}(h) = \sum_i d_{\mathrm{LE}}^2(\Sigma_i, \hat{\mu}_{h,-i}^{\mathrm{NW}}(X_i))$, where $\hat{\mu}_{h,-i}^{\mathrm{NW}}$ is the leave-one-out kernel estimator.
- **Rule-of-thumb bandwidths:** For the Gaussian kernel and a scalar predictor, $h \approx 0.2\, \mathrm{sd}(X)$ often provides a reasonable starting point in the SPD setting (analogous to Silverman's rule for density estimation).
- **Adaptive bandwidths:** In regions with sparse data (e.g., extreme stress levels), a variable bandwidth $h(x)$ that expands where data are scarce can stabilize the estimator — this is the metric analogue of nearest-neighbor bandwidth selection.
### From Global to Kernel: A Practical Guide
| Criterion | Prefer Global (Lecture 6) | Prefer Kernel (Lecture 7) |
|---|---|---|
| **True relationship** | Approximately linear in log-domain | Unknown or known to be nonlinear |
| **Interpretability** | Single "slope" interpretation | Local average — harder to summarize |
| **Extrapolation** | Possible (linear extrapolation) | Not recommended (no data → no fit) |
| **Boundary behavior** | Consistent everywhere | Boundary bias (see Lecture 8 for fix) |
| **Computational cost** | Same ($n$ weights per $x$) | Same ($n$ kernel evaluations per $x$) |
| **Dimension of $X$** | Degrades with $p$ (curse of dimensionality for $\hat{\Sigma}^{-1}$) | Degrades with $p$ (kernel sparsity) |
### Takeaway
Kernel Fréchet regression applies the same computational primitive as global Fréchet regression — weighted Fréchet means on the SPD manifold — to the same market risk data, but with a fundamentally different assumption: the relationship between stress and covariance structure is **local** rather than **global-linear**. This makes kernel regression the preferred tool when the true stress–risk relationship may involve thresholds, saturation, or other nonlinearities that a straight-line model cannot capture.
## Key Takeaways
- **Kernel Fréchet regression** replaces the global weight function $s(z, x)$ with localized kernel weights $K_h(X - x)$, yielding a local-constant smoother in the response space.
- The **target** $\mu_h^{\mathrm{NW}}(x)$ is a smoothed version of the true conditional Fréchet mean; they coincide only in the limit $h \to 0$.
- The **sample estimator** reduces to a weighted Fréchet mean with normalized kernel weights — computationally identical to the Fréchet mean problem, but with data-dependent weights.
- **Bias–variance tradeoff:** Small $h$ gives low bias but high variance (few effective data points); large $h$ gives low variance but high bias (oversmoothing).
- In a **Hadamard space**, the split proximal-point algorithms from Lecture 5 apply directly to the NW estimator: the component resolvent has the closed form $t_i = 2\lambda w_{in}^{\mathrm{NW}} / (1 + 2\lambda w_{in}^{\mathrm{NW}})$, and each split step is simply a geodesic interpolation toward one response $Y_i$.
- Kernel regression is **genuinely nonparametric** — it makes no global linearity assumption and can capture arbitrary smooth relationships.
## Exercises
1. **Bandwidth extremes.** What happens to $\hat{\mu}_h^{\mathrm{NW}}(x)$ as $h \to 0$ and as $h \to \infty$? Interpret both limits geometrically. <a href="javascript:void(0)" onclick="showSolution('l7-sol-1')" class="solution-link">📝 Show Solution</a>
2. **Kernel choice.** The Gaussian kernel has infinite support, while the Epanechnikov kernel has bounded support. Discuss the practical tradeoffs for Fréchet regression, particularly regarding computational cost and the existence of the weighted Fréchet mean. <a href="javascript:void(0)" onclick="showSolution('l7-sol-2')" class="solution-link">📝 Show Solution</a>
3. **Nadaraya–Watson as a weighted Fréchet mean.** Show directly that for $\mathcal{M} = \mathbb{R}^q$, the kernel Fréchet regression estimator reduces exactly to the componentwise Nadaraya–Watson estimator. <a href="javascript:void(0)" onclick="showSolution('l7-sol-3')" class="solution-link">📝 Show Solution</a>
4. **Local-constant bias at boundaries.** Local-constant kernel regression is known to have $O(h)$ (rather than $O(h^2)$) bias near the boundary of the predictor support. Does this phenomenon carry over to metric-space responses? Explain why or why not. <a href="javascript:void(0)" onclick="showSolution('l7-sol-4')" class="solution-link">📝 Show Solution</a>
5. **PPA implementation.** Implement the split proximal-point algorithm for the Nadaraya–Watson estimator in the SPD matrix space endowed with the log-Euclidean metric.
<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="l7-sol-1">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 1: Bandwidth Extremes</h4>
<button onclick="closeSolution('l7-sol-1')" class="solution-dialog-close">✕ Close</button>
</div>
<div class="solution-original"><strong>Exercise:</strong> What happens to $\hat{\mu}_h^{\mathrm{NW}}(x)$ as $h \to 0$ and as $h \to \infty$?</div>
<div class="solution-answer">
<strong>Solution:</strong>
**As $h \to 0$:** The kernel becomes a Dirac delta at $x$. Only observations with $X_i = x$ (exactly) receive nonzero weight. If there are multiple such observations, $\hat{\mu}_h^{\mathrm{NW}}(x)$ is their Fréchet mean. If there are none, the estimator is undefined. In practice, for very small $h$, only a few neighbors contribute, leading to extremely high variance.
**As $h \to \infty$:** The kernel becomes flat: $K_h(X_i - x) \to K(0)/h^p \to 0$, but the normalized weights $w_{in}$ approach $1/n$ for all $i$. Therefore $\hat{\mu}_h^{\mathrm{NW}}(x)$ converges to the **unconditional sample Fréchet mean** of all $Y_i$ — a single constant, independent of $x$. This is maximum bias (no localization) but minimum variance.
The optimal bandwidth balances these extremes.
</div>
</div>
</dialog>
<dialog id="l7-sol-2">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 2: Kernel Choice</h4>
<button onclick="closeSolution('l7-sol-2')" class="solution-dialog-close">✕ Close</button>
</div>
<div class="solution-original"><strong>Exercise:</strong> Discuss practical tradeoffs between Gaussian and Epanechnikov kernels for Fréchet regression.</div>
<div class="solution-answer">
<strong>Solution:</strong>
**Gaussian kernel (infinite support):**
- *Pro:* Smooth, all observations contribute, numerically stable.
- *Con:* Every evaluation requires $O(n)$ computation of distances and weights — no sparse representation.
- *Metric-space consideration:* The weighted Fréchet mean always uses all data points; some may have negligible weight but still contribute computationally.
**Epanechnikov kernel (bounded support):**
- *Pro:* Sparse — only observations within $h$ of $x$ contribute → faster computation for large $n$.
- *Con:* Discontinuous derivative at the boundary; may cause numerical issues in gradient-based optimization.
- *Metric-space consideration:* The effective sample size is smaller, which may lead to non-unique or poorly defined weighted Fréchet means when too few points are local.
</div>
</div>
</dialog>
<dialog id="l7-sol-3">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 3: Recovery of Nadaraya–Watson</h4>
<button onclick="closeSolution('l7-sol-3')" class="solution-dialog-close">✕ Close</button>
</div>
<div class="solution-original"><strong>Exercise:</strong> Show that for $\mathcal{M} = \mathbb{R}^q$, the kernel Fréchet regression estimator reduces to componentwise Nadaraya–Watson.</div>
<div class="solution-answer">
<strong>Solution:</strong>
For $\mathcal{M} = \mathbb{R}^q$ with $d(y, \omega) = \|y - \omega\|_2$, the weighted Fréchet objective separates componentwise:
$$
\sum_{i=1}^n w_i \|Y_i - \omega\|_2^2 = \sum_{j=1}^q \sum_{i=1}^n w_i (Y_{ij} - \omega_j)^2.
$$
For each component $j$, the minimizer of $\sum_i w_i (Y_{ij} - \omega_j)^2$ is the weighted average:
$$
\hat{\omega}_j = \frac{\sum_i w_i Y_{ij}}{\sum_i w_i}.
$$
With kernel weights $w_i = K_h(X_i - x)$, this is exactly the Nadaraya–Watson estimator for the $j$-th coordinate. Since the components are independent in the Euclidean metric, the vector of componentwise minimizers is the joint minimizer:
$$
\hat{\mu}_h^{\mathrm{NW}}(x) = \frac{\sum_i K_h(X_i - x) Y_i}{\sum_i K_h(X_i - x)}.
$$
</div>
</div>
</dialog>
<dialog id="l7-sol-4">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 4: Boundary Bias</h4>
<button onclick="closeSolution('l7-sol-4')" class="solution-dialog-close">✕ Close</button>
</div>
<div class="solution-original"><strong>Exercise:</strong> Does the $O(h)$ boundary bias of local-constant regression carry over to metric-space responses?</div>
<div class="solution-answer">
<strong>Solution:</strong>
**Yes**, the boundary bias phenomenon carries over. The reason is geometric:
- Near a boundary of the predictor support, the kernel weights are **asymmetric** — more weight falls on one side of $x$ than the other.
- The weighted Fréchet mean is then pulled toward the more heavily weighted side, introducing first-order ($O(h)$) bias even if the true regression function is locally linear.
- This is a consequence of the **local-constant** nature of the estimator: it approximates $\mu(x)$ by a constant in a neighborhood, and at boundaries that approximation is first-order wrong.
**Solution:** Local Fréchet regression (Lecture 8) corrects the boundary bias by fitting a local **linear** model, reducing bias from $O(h)$ to $O(h^2)$ everywhere, including at boundaries.
</div>
</div>
</dialog>
<dialog id="l7-sol-5">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 5: PPA Convergence</h4>
<button onclick="closeSolution('l7-sol-5')" class="solution-dialog-close">✕ Close</button>
</div>
<div class="solution-original"><strong>Exercise:</strong> Explain why $Q_x(\omega)$ is convex in a Hadamard space and why the PPA converges for any $\lambda_k \ge \lambda_{\min} > 0$.</div>
<div class="solution-answer">
<strong>Solution:</strong>
**Convexity:** In a Hadamard space, for any fixed $y$, the function $\omega \mapsto d^2(y, \omega)$ is **1-strongly convex** along geodesics. The NW weights satisfy $w_{in}^{\mathrm{NW}}(x, h) \ge 0$ and sum to 1, so
$$
Q_x(\omega) = \sum_{i=1}^n w_{in}^{\mathrm{NW}}(x, h)\, d^2(Y_i, \omega)
$$
is a convex combination of convex functions, hence itself 1-strongly convex. It therefore admits a unique minimizer $\hat{\mu}_h^{\mathrm{NW}}(x)$ in any Hadamard space.
**PPA convergence:** The PPA iteration is
$$
\omega_{k+1} = \arg\min_{\omega \in \mathcal{M}} \Bigl\{ Q_x(\omega) + \frac{1}{2\lambda_k} d^2(\omega_k, \omega) \Bigr\}.
$$
Each subproblem is $(1 + 1/\lambda_k)$-strongly convex, guaranteeing a unique solution at every step. In a Hadamard space, the resolvent (proximal) operator
$$
R_{\lambda_k}(\omega_k) = \arg\min_{\omega \in \mathcal{M}} \Bigl\{ Q_x(\omega) + \frac{1}{2\lambda_k} d^2(\omega_k, \omega) \Bigr\}
$$
is **firmly nonexpansive**. This means that for any two points $u, v \in \mathcal{M}$, the distance between their images under $R_{\lambda_k}$ is no larger than the distance between $u$ and $v$ — with a strict contraction unless both are already fixed points.
Because the fixed points of $R_{\lambda_k}$ are exactly the minimizers of $Q_x$ (for any $\lambda_k > 0$), firm nonexpansivity together with the existence of a unique fixed point implies that the iterates $\omega_{k+1} = R_{\lambda_k}(\omega_k)$ converge to $\hat{\mu}_h^{\mathrm{NW}}(x)$ for any starting point $\omega_0$.
**Why $\lambda_k \ge \lambda_{\min} > 0$ matters:** If $\lambda_k \to 0$, the proximal term $\frac{1}{2\lambda_k}d^2(\omega_k, \omega)$ dominates and the iterates may stall ($\omega_{k+1} \approx \omega_k$) before reaching the minimizer. Bounding $\lambda_k$ away from zero ensures that each step makes sufficient progress. In practice, a constant $\lambda_k \equiv 1$ or a slowly diminishing schedule such as $\lambda_k = 1/\sqrt{k+1}$ (with a positive floor) works well.
</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 7 Quiz — 10 Multiple-Choice Questions](../quizzes/lecture-07-quiz.qmd)**
## Further Reading
- @PetersenMueller2019 — Sections 2 and 5: the global-to-local progression of Fréchet regression.
- @DavisFletcherBullittJoshi2007 — Early work on manifold-valued kernel regression for shape and diffusion-tensor responses.
- @LinMuller2021TV — Total-variation regularized Fréchet regression, an alternative to kernel smoothing for piecewise-constant signals.