Lecture 8: Fréchet Regression — Local Fréchet Regression

Local linear smoothing in metric spaces

1 Learning Goals

By the end of this lecture, learners should be able to:

Explain why local-constant kernel regression suffers from boundary bias and how local linear regression corrects it.
Define the local moments $\nu_j(x,h)$ and the local weight function $s(z, x, h)$.
Derive the weight function $s(z, x, h)$ from the normal equations of Euclidean local linear regression.
Explain why $\mathbb{E}[s(X, x, h)] = 1$ and $\mathbb{E}[s(X, x, h)(X - x)] = 0$ are the two key properties enabling the metric-space generalization.
State the local Fréchet regression target and sample estimator.
Interpret the bias, stochastic error, and overall convergence rate of the local Fréchet regression estimator.
Derive the optimal bandwidth $h \asymp n^{-1/5}$ for the quadratic-curvature case $\beta_1 = \beta_2 = 2$.
Contrast local-constant (kernel) and local-linear Fréchet regression.

2 Why Local Linear?

Local-constant (kernel) Fréchet regression, discussed in Lecture 7, approximates the regression function as constant within each kernel window. This leads to:

$O(h^2)$ bias in the interior (acceptable).
$O(h)$ bias near boundaries of the predictor support (problematic).

Local linear regression fixes this by fitting a line within each window, reducing the bias to $O(h^2)$ everywhere — including at boundaries. Petersen and Müller (2019) extend this idea to metric spaces. For exposition, we consider the case $p = 1$ (scalar predictor).

3 Local Moments and Weight Functions

Let $K$ be a kernel density and $K_h(u) = h^{-1}K(u/h)$ for bandwidth $h > 0$. Define the local moments

\[ \nu_j(x, h) = \mathbb{E}\bigl\{K_h(X - x)(X - x)^j\bigr\}, \qquad j = 0, 1, 2, \]

and

\[ \sigma_0^2(x, h) = \nu_0(x, h) \nu_2(x, h) - \nu_1(x, h)^2. \]

First consider the ordinary Euclidean case $Y \in \mathbb{R}$. The local linear estimator at $x$ is the fitted intercept $\hat{\mu}(x) = \hat{\beta}_0$, where $(\hat{\beta}_0, \hat{\beta}_1)$ minimizes

\[ \frac{1}{n}\sum_{i=1}^n K_h(X_i - x)\bigl\{Y_i - \beta_0 - \beta_1(X_i - x)\bigr\}^2. \]

At the population level, the corresponding target $(\bar{\beta}_0, \bar{\beta}_1)$ minimizes

\[ \mathbb{E}\Bigl[K_h(X - x)\bigl\{Y - \beta_0 - \beta_1(X - x)\bigr\}^2\Bigr]. \]

4 From Euclidean Local Linear to the Weight Function

The key algebraic step — and the one that makes the metric-space generalization possible — is to rewrite the local linear intercept as a weighted expectation with a data-dependent weight function $s(z, x, h)$. We now derive this transformation in detail.

4.1 Step 1: The Normal Equations

Define the cross-moments

\[ r_j(x, h) = \mathbb{E}\bigl\{K_h(X - x)(X - x)^j Y\bigr\}, \qquad j = 0, 1. \]

At the population level, $(\bar{\beta}_0, \bar{\beta}_1)$ minimizes

\[ Q(\beta_0, \beta_1) = \mathbb{E}\Bigl[K_h(X - x)\bigl\{Y - \beta_0 - \beta_1(X - x)\bigr\}^2\Bigr]. \]

Differentiating with respect to $\beta_0$ and $\beta_1$ and setting to zero yields the normal equations:

\[ \begin{aligned} \mathbb{E}\bigl[K_h(X - x)\{Y - \bar{\beta}_0 - \bar{\beta}_1(X - x)\}\bigr] &= 0, \\ \mathbb{E}\bigl[K_h(X - x)(X - x)\{Y - \bar{\beta}_0 - \bar{\beta}_1(X - x)\}\bigr] &= 0. \end{aligned} \]

Expanding and using the definitions of $\nu_j$ and $r_j$:

\[ \boxed{ \begin{aligned} \nu_0(x, h)\,\bar{\beta}_0 + \nu_1(x, h)\,\bar{\beta}_1 &= r_0(x, h), \\ \nu_1(x, h)\,\bar{\beta}_0 + \nu_2(x, h)\,\bar{\beta}_1 &= r_1(x, h). \end{aligned}} \tag{1} \]

4.2 Step 2: Solving for the Intercept

The system (1) is a $2 \times 2$ linear system whose coefficient matrix has determinant $\sigma_0^2 = \nu_0\nu_2 - \nu_1^2 > 0$ (strictly positive since $K \ge 0$ and $X$ has a continuous density). By Cramer’s rule, the intercept is

\[ \tilde{\mu}(x) = \bar{\beta}_0 = \frac{\det\begin{pmatrix} r_0 & \nu_1 \\ r_1 & \nu_2 \end{pmatrix}}{\det\begin{pmatrix} \nu_0 & \nu_1 \\ \nu_1 & \nu_2 \end{pmatrix}} = \frac{\nu_2(x, h)\, r_0(x, h) - \nu_1(x, h)\, r_1(x, h)}{\sigma_0^2(x, h)}. \tag{2} \]

4.3 Step 3: The Weight Function Emerges

Now substitute the definition $r_j(x, h) = \mathbb{E}[K_h(X - x)(X - x)^j Y]$ into (2):

\[ \begin{aligned} \tilde{\mu}(x) &= \frac{\nu_2\,\mathbb{E}[K_h(X - x)Y] - \nu_1\,\mathbb{E}[K_h(X - x)(X - x)Y]}{\sigma_0^2} \\[4pt] &= \mathbb{E}\!\left[\frac{K_h(X - x)\bigl\{\nu_2(x, h) - \nu_1(x, h)(X - x)\bigr\}}{\sigma_0^2(x, h)}\, Y\right]. \end{aligned} \]

This is of the form $\tilde{\mu}(x) = \mathbb{E}[s(X, x, h)\,Y]$, where the weight function is

\[ \boxed{s(z, x, h) = \frac{K_h(z - x)\bigl\{\nu_2(x, h) - \nu_1(x, h)(z - x)\bigr\}}{\sigma_0^2(x, h)}}. \tag{3} \]

4.4 Step 4: Properties of the Weight Function

Two properties of $s(z, x, h)$ are essential for what follows.

Sum-to-one. Take $Y \equiv 1$ in the derivation above. Then $r_0 = \nu_0$ and $r_1 = \nu_1$, so (2) gives $\tilde{\mu}(x) = 1$. But $\tilde{\mu}(x) = \mathbb{E}[s(X, x, h) \cdot 1]$, hence

\[ \boxed{\mathbb{E}[s(X, x, h)] = 1}. \tag{4} \]

First-moment annihilation. Take $Y = X - x$. Then $r_0 = \nu_1$ and $r_1 = \nu_2$, so (2) gives $\tilde{\mu}(x) = 0$. But $\tilde{\mu}(x) = \mathbb{E}[s(X, x, h)(X - x)]$, hence

\[ \boxed{\mathbb{E}[s(X, x, h)(X - x)] = 0}. \tag{5} \]

Equation (5) is the mathematical expression of the bias-correction mechanism: the weights are orthogonal to linear trends, so the local-linear estimator is exactly unbiased for linear conditional mean functions — a property the local-constant (kernel) estimator lacks.

4.5 Step 5: Reformulation as Weighted Least Squares

Because $\mathbb{E}[s(X, x, h)] = 1$, we can add a constant term (with respect to $y$) without changing the minimizer:

\[ \begin{aligned} \tilde{\mu}(x) &= \mathbb{E}[s(X, x, h) Y] \\ &= \arg\min_{y \in \mathbb{R}} \Bigl\{\mathbb{E}[s(X, x, h) Y^2] - 2y\,\mathbb{E}[s(X, x, h) Y] + y^2\Bigr\} \\ &= \arg\min_{y \in \mathbb{R}} \mathbb{E}\bigl[s(X, x, h)(Y - y)^2\bigr]. \end{aligned} \]

Thus

\[ \boxed{\tilde{\mu}(x) = \arg\min_{y \in \mathbb{R}} \mathbb{E}\bigl\{s(X, x, h)(Y - y)^2\bigr\}}. \tag{6} \]

Why this form matters

Equation (6) expresses the local linear target purely in terms of a weighted expectation of a squared distance. This is the form that generalizes: the weight function $s(z, x, h)$ is built entirely from the predictor $X$ and depends on the response $Y$ only through the squared-error criterion. To extend to a metric-space-valued response, we simply:

Replace the squared Euclidean distance $(Y - y)^2$ by the squared metric $d^2(Y, \omega)$,
Minimize over $\omega$ in the metric space $\mathcal{M}$ rather than over $y \in \mathbb{R}$.

The weight function $s(z, x, h)$ remains unchanged — it is computed from the scalar predictors exactly as in the Euclidean case.

5 The Local Fréchet Regression Target and Estimator

Now consider data $(X_i, Y_i) \in \mathbb{R} \times \mathcal{M}$, $i = 1, \dots, n$, where $\mathcal{M}$ is a metric space with distance $d$.

Definition 1 Let

\[ \widetilde{L}_h(\omega, x) = \mathbb{E}\{s(X, x, h) d^2(Y, \omega)\}, \]

The local Fréchet regression target is

\[ \tilde{\mu}(x) = \arg\min_{\omega \in \mathcal{M}} \widetilde{L}_h(\omega, x). \]

This target is not exactly the conditional Fréchet mean $\mu(x)$; it is the metric-space analogue of the usual local linear population target. The key fact proved by Petersen and Müller is that the induced weighted measure differs from the true conditional law by $O(h^2)$, so minimizing $\widetilde{L}_h(\cdot, x)$ approximates minimizing the conditional Fréchet function.

Given data, define empirical local moments

\[ \hat{\nu}_j(x, h) = \frac{1}{n}\sum_{i=1}^n K_h(X_i - x)(X_i - x)^j, \qquad j = 0, 1, 2, \]

\[ \hat{\sigma}_0^2(x, h) = \hat{\nu}_0(x, h) \hat{\nu}_2(x, h) - \hat{\nu}_1(x, h)^2, \]

and weights

\[ s_{in}(x, h) = \frac{K_h(X_i - x)\{\hat{\nu}_2(x, h) - \hat{\nu}_1(x, h)(X_i - x)\}}{\hat{\sigma}_0^2(x, h)}. \]

The sample local objective is

\[ \hat{L}_n(\omega, x) = \frac{1}{n}\sum_{i=1}^n s_{in}(x, h) d^2(Y_i, \omega), \]

and the local Fréchet regression estimator is

\[ \hat{\mu}(x) = \arg\min_{\omega \in \mathcal{M}} \hat{L}_n(\omega, x). \]

Important distinction

Unlike kernel regression weights, the local linear weights $s_{in}(x, h)$ can be negative. This is the price for the bias correction — the weights automatically adjust for the local slope. The empirical weights satisfy $\sum_i s_{in}(x, h) = 1$ and $\sum_i s_{in}(x, h)(X_i - x) = 0$.

6 Bias, Stochastic Error, and Convergence Rates

The following regularity assumptions are from Petersen and Müller (2019), Section 4:

Assumptions (K0), (L0)–(L3), and (P1)

(K0) Kernel regularity. $K$ is a symmetric probability density with finite moments $|K_{14}| < \infty$ and $|K_{26}| < \infty$, where $K_{kj} = \int_{\mathbb{R}} K(u)^k u^j\, du$.

(L0) Existence and separation. $\mu(x)$ and $\tilde{\mu}(x)$ exist and are unique. Also, $\hat{\mu}(x)$ exists and is unique, almost surely. For any $\epsilon > 0$,

\[ \inf_{d(\omega, \mu(x)) > \epsilon}\{F(\omega, x) - F(\mu(x), x)\} > 0, \qquad \liminf_{n \to \infty} \inf_{d(\omega, \tilde{\mu}(x)) > \epsilon}\{\widetilde{L}_n(\omega, x) - \widetilde{L}_n(\tilde{\mu}(x), x)\} > 0. \]

(L1) Smooth predictor law. The marginal density $f$ of $X$ and the conditional densities $g_y$ of $X \mid Y = y$ exist and are twice continuously differentiable, with $\sup_{x, y} |g_y''(x)| < \infty$.

(L2) Bias curvature. There exist $\eta_1 > 0$, $C_1 > 0$, and $\beta_1 > 1$ such that

\[ F(\omega, x) - F(\mu(x), x) \ge C_1 d(\omega, \mu(x))^{\beta_1} \quad \text{whenever } d(\omega, \mu(x)) < \eta_1. \]

(L3) Stochastic curvature. There exist $\eta_2 > 0$, $C_2 > 0$, and $\beta_2 > 1$ such that

\[ \liminf_{n \to \infty} \bigl[\widetilde{L}_n(\omega, x) - \widetilde{L}_n(\tilde{\mu}(x), x)\bigr] \ge C_2 d(\omega, \tilde{\mu}(x))^{\beta_2} \quad \text{whenever } d(\omega, \tilde{\mu}(x)) < \eta_2. \]

(P1) Local entropy. Let $B_\delta(\mu(x)) \subset \mathcal{M}$ be the ball of radius $\delta$ centered at $\mu(x)$, and let $N(\varepsilon, B_\delta(\mu(x)), d)$ denote its covering number by balls of radius $\varepsilon$. Then

\[ \int_0^1 \bigl\{1 + \log N(\delta\varepsilon, B_\delta(\mu(x)), d)\bigr\}^{1/2}\, d\varepsilon = O(1) \qquad \text{as } \delta \downarrow 0. \]

This controls the local complexity of the metric space near the target, ensuring that the stochastic fluctuation of the empirical objective on shrinking neighborhoods of $\mu(x)$ remains bounded. For finite-dimensional manifolds (e.g., $\mathbb{S}^d$, $\mathcal{S}_{++}^m$), the covering number grows as $\varepsilon^{-\dim}$, and the integral condition is automatically satisfied. For infinite-dimensional function spaces, it may fail and requires verification.

Theorem 1 Suppose $(\mathcal{M}, d)$ is totally bounded.

Bias: If (K0), (L0), (L1), (L2), and (P1) hold, then as $h = h_n \to 0$, \[ d(\mu(x), \tilde{\mu}(x)) = O\bigl(h^{2/(\beta_1 - 1)}\bigr). \]
Stochastic error: If (K0), (L0), (L3), and (P1) hold, and $h \to 0$ with $nh \to \infty$, then \[ d(\tilde{\mu}(x), \hat{\mu}(x)) = O_p\bigl((nh)^{-1/[2(\beta_2 - 1)]}\bigr). \]
Overall rate: Under all assumptions, \[ d(\mu(x), \hat{\mu}(x)) = O_p\bigl(h^{2/(\beta_1 - 1)} + (nh)^{-1/[2(\beta_2 - 1)]}\bigr). \]

Choosing $h = n^{-\gamma}$ with $\gamma^* = (\beta_1 - 1)/(4\beta_2 + \beta_1 - 5)$ yields the optimal rate

\[ d(\mu(x), \hat{\mu}(x)) = O_p\bigl(n^{-2/(\beta_1 + 4\beta_2 - 5)}\bigr). \]

The quadratic-curvature case $\beta_1 = \beta_2 = 2$

In the important special case where both curvature exponents are 2 (which holds for SPD spaces, BHV space, and Euclidean spaces):

\[ d(\mu(x), \hat{\mu}(x)) = O_p\bigl(h^2 + (nh)^{-1/2}\bigr), \]

and the optimal bandwidth $h \asymp n^{-1/5}$ yields the familiar nonparametric rate $n^{-2/5}$.

This shows that Fréchet regression inherits the same bias–variance tradeoff as classical local linear regression, but with the exponents $\beta_1$ and $\beta_2$ determined by the local geometry of the response space.

7 Comparison: Local-Constant vs. Local-Linear

Property	Local-Constant (Kernel)	Local-Linear
Bias (interior)	$O(h^2)$	$O(h^2)$
Bias (boundary)	$O(h)$	$O(h^2)$
Weights	Always nonnegative	Can be negative
Finite-sample	Simpler, more stable	Better at boundaries
Computational	Weighted Fréchet mean	Weighted Fréchet mean (same complexity)

8 Interactive Exploration: Global vs. Kernel vs. Local-Linear on the Sphere

The following demo compares all three Fréchet regression methods on the positive quadrant (octant) of the unit sphere $\mathbb{S}^2$:

Global Fréchet regression (Lecture 6) — the metric analogue of a global linear model.
Kernel (local-constant) Fréchet regression (Lecture 7) — the metric analogue of Nadaraya–Watson smoothing.
Local-linear Fréchet regression (this lecture) — the metric analogue of local linear regression.

The true regression curve is nonlinear (sinusoidal), so the global fit is misspecified. The kernel fit captures the nonlinearity but suffers from boundary bias. The local-linear fit corrects the boundary bias by adjusting weights to account for the local slope — some weights become negative as a result.

Data generation. Same as Lectures 6–7: $X_i \sim \text{Uniform}(-1, 1)$, true $\mu(x)$ is a sinusoidal curve on the sphere octant, and $Y_i$ is generated with geodesic noise.

Three weight functions compared:

Method	Weights $w_i(x)$	Can be negative?
Global (L6)	$1 + \frac{(X_i - \bar{X})(x - \bar{X})}{\hat{\sigma}_X^2}$	Yes
Kernel (L7)	$\frac{K_h(X_i - x)}{\sum_j K_h(X_j - x)}$	No
Local-linear (L8)	$\frac{K_h(X_i - x)\{\hat{\nu}_2 - \hat{\nu}_1(X_i - x)\}}{\hat{\nu}_0\hat{\nu}_2 - \hat{\nu}_1^2}$	Yes

All three then compute the weighted Fréchet mean on $\mathbb{S}^2$ via Riemannian gradient descent.

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 (L6)
Orange curve = kernel / local-constant fit (L7)
Blue curve = local-linear fit (L8)
Gold diamond = unconditional Fréchet mean of all $Y_i$

Code

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_sphere_view = html`
<style>
  .lfr-control-panel {
    width: 100%;
    margin: 0 0 12px;
    font-family: system-ui, sans-serif;
    font-size: 0.85em;
    container-type: inline-size;
  }
  .lfr-control-panel > *,
  .lfr-slider-grid > * {
    min-width: 0;
    margin: 0;
  }
  .lfr-slider-grid {
    display: flex;
    flex-wrap: wrap;
    gap: 6px 20px;
    margin-top: 6px;
  }
  .lfr-slider-grid > * {
    flex: 1 1 calc((100% - 40px) / 3);
  }
  .lfr-slider-grid input[type="number"] {
    width: 7.5rem !important;
  }
  @container (max-width: 700px) {
    .lfr-slider-grid > * {
      flex-basis: calc((100% - 20px) / 2);
    }
  }
  @container (max-width: 480px) {
    .lfr-slider-grid > * {
      flex-basis: 100%;
    }
  }
</style>
<div class="lfr-control-panel">
  <div>${reg_type_control}</div>
  <div class="lfr-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 lectures 6–7) ----
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)
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 };
}

// ---- Gaussian kernel ----
function gaussKernel(u) { return Math.exp(-0.5 * u * u); }

// ---- 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 (local-constant) Fréchet regression (Lecture 7) ----
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;
}

// ---- Local-linear Fréchet regression (Lecture 8) ----
function computeLocalLinearCurve(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 = [];
    for (var i = 0; i < n; i++) {
      rawW.push(gaussKernel((X[i] - x) / h));
    }
    // Compute local moments
    var nu0 = 0, nu1 = 0, nu2 = 0;
    for (var i = 0; i < n; i++) {
      var di = X[i] - x;
      nu0 += rawW[i];
      nu1 += rawW[i] * di;
      nu2 += rawW[i] * di * di;
    }
    nu0 /= n; nu1 /= n; nu2 /= n;
    var sigma02 = nu0 * nu2 - nu1 * nu1;

    var weights;
    if (Math.abs(sigma02) > 1e-10) {
      // Local-linear weights (can be negative — the bias-correction mechanism)
      weights = X.map(function(xi, i) {
        return rawW[i] * (nu2 - nu1 * (xi - x)) / sigma02;
      });
    } else {
      // Fallback to kernel (local-constant)
      var wSum = rawW.reduce(function(a,b) { return a+b; }, 0);
      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);
localLinCurve = computeLocalLinearCurve(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);

// Local-linear at eval
nu0e = rawEvalW.reduce(function(a,b) { return a+b; }, 0) / regData.X.length;
nu1e = rawEvalW.reduce(function(s, w, i) { return s + w * (regData.X[i] - eval_x); }, 0) / regData.X.length;
nu2e = rawEvalW.reduce(function(s, w, i) { return s + w * Math.pow(regData.X[i] - eval_x, 2); }, 0) / regData.X.length;
sigma02e = nu0e * nu2e - nu1e * nu1e;
llWeightsEval = (Math.abs(sigma02e) > 1e-10)
  ? regData.X.map(function(xi, i) { return rawEvalW[i] * (nu2e - nu1e * (xi - eval_x)) / sigma02e; })
  : nWeightsEval;
localLinAtEval = weightedSphereFrechetMean(regData.Y, llWeightsEval, 200, 0.3);

trueAtEval = trueRegFunc(eval_x, reg_type);
globalErr = sphereDist(globalAtEval, trueAtEval);
kernelErr = sphereDist(kernelAtEval, trueAtEval);
localLinErr = sphereDist(localLinAtEval, trueAtEval);

// Overall Fréchet mean
unifWeights = regData.Y.map(function() { return 1; });
overallMean = weightedSphereFrechetMean(regData.Y, unifWeights, 200, 0.3);

// Diagnostics
effN = rawEvalW.filter(function(w) { return w > 1e-6; }).length;
llNegCount = llWeightsEval.filter(function(w) { return w < 0; }).length;
llWeightSum = llWeightsEval.reduce(function(a,b) { return a+b; }, 0);

// ---- 3D Visualization of the sphere octant ----
function build3DView(regData, trueCurve, globalCurve, kernelCurve, localLinCurve, xGridEval, globalAtEval, kernelAtEval, localLinAtEval, 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)
  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, thinner)
  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.2" stroke-linecap="round" stroke-linejoin="round" opacity="0.7"/>');

  // Kernel / local-constant curve (orange)
  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="#f08c00" stroke-width="2.6" stroke-linecap="round" stroke-linejoin="round" opacity="0.8"/>');

  // Local-linear curve (blue, thicker — the featured method)
  var llPtsStr = localLinCurve.map(function(p) { var proj = project(p); return proj[0].toFixed(1) + "," + proj[1].toFixed(1); }).join(" ");
  parts.push('<polyline points="' + llPtsStr + '" fill="none" stroke="#1971c2" stroke-width="3.2" stroke-linecap="round" stroke-linejoin="round" opacity="0.9"/>');

  // 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 lEvalProj = project(localLinAtEval);
  var tEvalProj = project(trueAtEval);

  // Dashed connectors from each fit 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.0" stroke-dasharray="3,2" opacity="0.45"/>');
  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.0" stroke-dasharray="3,2" opacity="0.45"/>');
  parts.push('<line x1="' + lEvalProj[0].toFixed(1) + '" y1="' + lEvalProj[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="6" fill="none" stroke="#e53935" stroke-width="2.5" opacity="0.7"/>');
  parts.push('<circle cx="' + kEvalProj[0].toFixed(1) + '" cy="' + kEvalProj[1].toFixed(1) +
    '" r="6" fill="none" stroke="#f08c00" stroke-width="2.5" opacity="0.8"/>');
  parts.push('<circle cx="' + lEvalProj[0].toFixed(1) + '" cy="' + lEvalProj[1].toFixed(1) +
    '" r="7" fill="none" stroke="#1971c2" stroke-width="3.5"/>');
  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 = 35;
  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+16) + '" x2="' + (lx+14) + '" y2="' + (ly0+16) +
    '" stroke="#2e7d32" stroke-width="2.5" stroke-dasharray="5,3"/>');
  parts.push('<text x="' + (lx+22) + '" y="' + (ly0+20) + '" font-size="9" fill="#37474f">True μ(x)</text>');
  parts.push('<line x1="' + (lx-6) + '" y1="' + (ly0+32) + '" x2="' + (lx+14) + '" y2="' + (ly0+32) +
    '" stroke="#e53935" stroke-width="2.0"/>');
  parts.push('<text x="' + (lx+22) + '" y="' + (ly0+36) + '" font-size="9" fill="#37474f">Global (L6)</text>');
  parts.push('<line x1="' + (lx-6) + '" y1="' + (ly0+48) + '" x2="' + (lx+14) + '" y2="' + (ly0+48) +
    '" stroke="#f08c00" stroke-width="2.2"/>');
  parts.push('<text x="' + (lx+22) + '" y="' + (ly0+52) + '" font-size="9" fill="#37474f">Kernel (L7)</text>');
  parts.push('<line x1="' + (lx-6) + '" y1="' + (ly0+64) + '" x2="' + (lx+14) + '" y2="' + (ly0+64) +
    '" stroke="#1971c2" stroke-width="2.8"/>');
  parts.push('<text x="' + (lx+22) + '" y="' + (ly0+68) + '" font-size="9" fill="#37474f">Local-linear (L8)</text>');
  parts.push('<polygon points="' + lx + ',' + (ly0+78) + ' ' + (lx+4) + ',' + (ly0+83) +
    ' ' + lx + ',' + (ly0+88) + ' ' + (lx-4) + ',' + (ly0+83) + '" fill="#f9a825"/>');
  parts.push('<text x="' + (lx+12) + '" y="' + (ly0+87) + '" 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:700px;">\n' +
    parts.join('\n') + '\n</svg>';
}

sphere_l8_view = html`<div style="width:80%; max-width:700px; margin:0 auto; text-align:center;">${build3DView(regData, trueCurve, globalCurve, kernelCurve, localLinCurve, xGridEval, globalAtEval, kernelAtEval, localLinAtEval, 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); });
localLinAngles = localLinCurve.map(function(p) { return sphericalAngles(p); });
dataAngles = regData.Y.map(function(yi) { return sphericalAngles(yi); });

evalTrueAng = sphericalAngles(trueAtEval);
evalGlobalAng = sphericalAngles(globalAtEval);
evalKernelAng = sphericalAngles(kernelAtEval);
evalLocalLinAng = sphericalAngles(localLinAtEval);

function buildUnwrappedPlot(title, yLabel, trueVals, globalVals, kernelVals, llVals, dataVals, xVals, xGrid, yMin, yMax) {
  var W = 620, 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"/>');

  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.4"/>');
  });

  // True curve (green, dashed)
  var line = "";
  for (var j = 0; j < xGrid.length; j++) {
    line += (j === 0 ? 'M' : 'L') + " " + xToSvg(xGrid[j]).toFixed(1) + " " + yToSvg(trueVals[j]).toFixed(1);
  }
  parts.push('<path d="' + line + '" fill="none" stroke="#2e7d32" stroke-width="2.2" stroke-dasharray="7,4"/>');

  // Global (red)
  line = "";
  for (var j = 0; j < xGrid.length; j++) {
    line += (j === 0 ? 'M' : 'L') + " " + xToSvg(xGrid[j]).toFixed(1) + " " + yToSvg(globalVals[j]).toFixed(1);
  }
  parts.push('<path d="' + line + '" fill="none" stroke="#e53935" stroke-width="1.8" opacity="0.7"/>');

  // Kernel (orange)
  line = "";
  for (var j = 0; j < xGrid.length; j++) {
    line += (j === 0 ? 'M' : 'L') + " " + xToSvg(xGrid[j]).toFixed(1) + " " + yToSvg(kernelVals[j]).toFixed(1);
  }
  parts.push('<path d="' + line + '" fill="none" stroke="#f08c00" stroke-width="2.0"/>');

  // Local-linear (blue, thicker)
  line = "";
  for (var j = 0; j < xGrid.length; j++) {
    line += (j === 0 ? 'M' : 'L') + " " + xToSvg(xGrid[j]).toFixed(1) + " " + yToSvg(llVals[j]).toFixed(1);
  }
  parts.push('<path d="' + line + '" fill="none" stroke="#1971c2" stroke-width="2.4"/>');

  // Evaluation point
  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"/>');

  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>');

  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 - 178, ly = margin.top + 4;
  parts.push('<line x1="' + lx + '" y1="' + ly + '" x2="' + (lx+12) + '" y2="' + ly +
    '" stroke="#2e7d32" stroke-width="1.6" stroke-dasharray="4,3"/>');
  parts.push('<text x="' + (lx+15) + '" y="' + (ly+4) + '" font-size="7.5" fill="#616161">True</text>');
  parts.push('<line x1="' + (lx+44) + '" y1="' + ly + '" x2="' + (lx+56) + '" y2="' + ly +
    '" stroke="#e53935" stroke-width="1.6"/>');
  parts.push('<text x="' + (lx+59) + '" y="' + (ly+4) + '" font-size="7.5" fill="#616161">Global</text>');
  parts.push('<line x1="' + (lx+90) + '" y1="' + ly + '" x2="' + (lx+102) + '" y2="' + ly +
    '" stroke="#f08c00" stroke-width="1.8"/>');
  parts.push('<text x="' + (lx+105) + '" y="' + (ly+4) + '" font-size="7.5" fill="#616161">Kernel</text>');
  parts.push('<line x1="' + (lx+138) + '" y1="' + ly + '" x2="' + (lx+150) + '" y2="' + ly +
    '" stroke="#1971c2" stroke-width="2"/>');
  parts.push('<text x="' + (lx+153) + '" y="' + (ly+4) + '" font-size="7.5" fill="#616161">LocalLin</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:700px;">\n' +
    parts.join('\n') + '\n</svg>';
}

unwrapped_l8_view = html`
<div style="width:80%; max-width:700px; 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; }),
    localLinAngles.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; }),
    localLinAngles.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_l8_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)}, h = ${bandwidth.toFixed(2)})</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 style="border-bottom:1px solid #dee2e6;">
        <td style="padding:5px 8px;" colspan="2"><b>Local-linear (L8)</b></td>
      </tr>
      <tr><td style="padding:4px 8px;">Local-linear μ̂(x)</td>
        <td style="padding:4px 8px; text-align:right; color:#1971c2; font-weight:bold;">θ=${evalLocalLinAng.theta.toFixed(3)}, φ=${evalLocalLinAng.phi.toFixed(3)} rad</td></tr>
      <tr><td style="padding:4px 8px;">Local-linear error</td>
        <td style="padding:4px 8px; text-align:right; color:#1971c2; font-weight:bold;">${localLinErr.toFixed(5)} rad (${(localLinErr*180/Math.PI).toFixed(3)}°)</td></tr>
      <tr><td style="padding:4px 8px;">Negative weights</td>
        <td style="padding:4px 8px; text-align:right;">${llNegCount} of ${regData.X.length} (${(llNegCount/regData.X.length*100).toFixed(1)}%)</td></tr>
      <tr style="border-bottom:1px solid #dee2e6;">
        <td style="padding:5px 8px;" colspan="2"><b>Kernel / local-constant (L7)</b></td>
      </tr>
      <tr><td style="padding:4px 8px;">Kernel μ̂(x)</td>
        <td style="padding:4px 8px; text-align:right; color:#f08c00;">θ=${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:#f08c00;">${kernelErr.toFixed(5)} rad (${(kernelErr*180/Math.PI).toFixed(3)}°)</td></tr>
      <tr style="border-bottom:1px solid #dee2e6;">
        <td style="padding:5px 8px;" colspan="2"><b>Global (L6)</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;">${globalErr.toFixed(5)} rad (${(globalErr*180/Math.PI).toFixed(3)}°)</td></tr>
      <tr style="border-bottom:1px solid #dee2e6;">
        <td style="padding:5px 8px;" colspan="2"><b>Common</b></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><td style="padding:4px 8px;">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_sphere_view}${sphere_l8_view}${unwrapped_l8_view}${statistics_l8_view}</div>`

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

(u)

(v)

(w)

(x)

(y)

(z)

({)

(|)

(})

(~)

()

( )

(¡)

(¢)

(£)

(¤)

(¥)

(¦)

(§)

(¨)

(©)

(ª)

(«)

(¬)

()

Figure 1: Interactive: Global vs. Kernel vs. Local-Linear Fréchet regression on the sphere octant

Try these experiments

Compare all three methods in the 3D view: Global (red) is a straight-line analogue — it misses the sinusoidal oscillation entirely. Kernel (orange) tracks the curve but can be pulled toward the interior at boundaries. Local-linear (blue) tracks the curve everywhere, including near the edges.
Toggle regression types: “Varying θ only” and “Varying φ only” isolate the nonlinearity to one coordinate — observe which methods capture it best in the 2D unwrapped views.
Look near the edges ($x = \pm 1$): Local-linear weights become negative (reported in the stats panel) — this is the mechanism that corrects boundary bias. The kernel weights remain nonnegative and exhibit boundary bias as a result.
Decrease $h$ (0.15–0.25): All local fits become more wiggly. The local-linear advantage at boundaries becomes more visible.
Increase $h$ (0.8–1.2): Kernel and local-linear both approach the unconditional Fréchet mean (gold diamond). At very large $h$, all three methods nearly coincide at $\bar{X}$ (since they all pass through the Fréchet mean there).
Increase $n$ (200+): All fits improve, but the global fit’s bias persists — no amount of data can fix model misspecification.
Study the negative weight count: Local-linear weights are negative near boundaries and in high-curvature regions. This is the mathematical mechanism that makes local-linear unbiased for linear functions — exactly as in classical local polynomial regression.

9 Application: Stress-Testing at the Boundary — Why Local-Linear Matters

Lecture 6 gave us a global linear model of how market covariance changes with stress. Lecture 7 replaced the linearity assumption with kernel smoothing, capturing nonlinear relationships but introducing $O(h)$ bias at the edges of the predictor range. This lecture’s local-linear method corrects that boundary bias — and in the market risk context, the boundaries are precisely where we care most.

9.1 The Boundary Problem in Risk Management

In the market risk application, the predictor is realized SPY volatility $X_t = \hat{\sigma}_t$, and the response is a rolling covariance matrix $\Sigma_t \in \mathcal{S}_{++}^p$. Stress-testing asks: what is the typical covariance structure at an extreme stress level $x^*$ near or beyond the edge of observed data?

At such boundary points:

Global Fréchet regression (Lecture 6) extrapolates linearly — its weights grow with distance from $\bar{X}$, producing estimates that can be far from the local data.
Kernel Fréchet regression (Lecture 7) uses only nearby observations — but at the boundary, all nearby observations lie on one side of $x^*$. The kernel-weighted average is pulled inward toward the interior, creating $O(h)$ boundary bias.
Local-linear Fréchet regression (this lecture) corrects this by allowing weights to become negative for observations on the dense side, effectively extrapolating the local slope back to the boundary. The bias drops to $O(h^2)$ everywhere.

In stress-testing, an $O(h)$ error in the estimated covariance matrix can translate to a material misestimate of portfolio risk at the very moment it matters most.

9.2 Local-Linear Fréchet Regression on SPD Matrices

Under the log-Euclidean metric, the local-linear estimator at stress level $x$ is

\[ \log \hat{\mu}_{\text{LL}}(x) = \frac{1}{n}\sum_{i=1}^n s_{in}(x, h)\, \log \Sigma_i, \qquad s_{in}(x, h) = \frac{K_h(X_i - x)\{\hat{\nu}_2(x, h) - \hat{\nu}_1(x, h)(X_i - x)\}}{\hat{\nu}_0(x, h)\hat{\nu}_2(x, h) - \hat{\nu}_1(x, h)^2}, \]

where $\hat{\nu}_j(x, h) = \frac{1}{n}\sum_{i=1}^n K_h(X_i - x)(X_i - x)^j$ are the empirical local moments. As in the global and kernel cases, the computational core is simply a weighted average of matrix logarithms — the geometry is handled by the log-Euclidean framework, and the regression method only changes the weights.

Three methods, one computational primitive

All three Fréchet regression methods share the same form under the log-Euclidean metric:

\[\log \hat{\mu}(x) = \sum_{i=1}^n w_i(x)\, \log \Sigma_i, \qquad \sum_i w_i(x) = 1.\]

Method	Weight $w_i(x)$	Sign	Boundary bias
Global (L6)	$\frac{1}{n}\bigl[1 + \frac{(X_i - \bar{X})(x - \bar{X})}{\hat{\sigma}_X^2}\bigr]$	Any	$O(1)$ if extrapolating
Kernel (L7)	$\frac{K_h(X_i - x)}{\sum_j K_h(X_j - x)}$	$\ge 0$	$O(h)$
Local-linear (L8)	$\frac{K_h(X_i - x)\{\hat{\nu}_2 - \hat{\nu}_1(X_i - x)\}}{n(\hat{\nu}_0\hat{\nu}_2 - \hat{\nu}_1^2)}$	Any	$O(h^2)$

The only difference between methods is the weight function — a purely predictor-side calculation. Switching between them requires no new geometry code.

9.3 What the Local-Linear Correction Buys You

At an extreme stress level $x^*$ (say, the 95th percentile of observed volatility):

Kernel regression sees mostly moderate-stress observations to the left of $x^*$. Its estimate is a weighted average of their covariance matrices — biased toward calmer-market structures. The resulting risk model underestimates the stress-typical correlations.
Local-linear regression sees the same data but fits a line through them. If covariances are trending upward with stress, the negative weights on moderate-stress observations pull the estimate further in the stress direction, extrapolating the trend. The risk model reflects what the data suggest about the boundary, not just what they show at the interior.
Global regression also extrapolates, but with a single straight line fitted to all data — it cannot adapt to local trends that may differ at the left and right boundaries.

The practical consequence: for a portfolio weight vector $w$, the stress-conditional variance estimate $\sigma_P^2(x^*) = w^\top \hat{\mu}_{\text{LL}}(x^*) w$ from local-linear regression is the most trustworthy of the three at extreme stress levels, because it alone combines local adaptivity with boundary debiasing.

9.4 Interactive Exploration: All Three Methods on Simulated SPD Data

The following demo extends the application from Lecture 7 by adding the local-linear estimator. We simulate $2 \times 2$ SPD matrices with a nonlinear relationship to a scalar stress predictor, and fit all three methods. Pay attention to the edges ($x \approx 0$ and $x \approx 1$) where local-linear should outperform kernel regression.

Code

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.18, 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_spd_view = html`
<style>
  .lfr-spd-slider-grid {
    display: flex;
    flex-wrap: wrap;
    gap: 6px 20px;
    width: 100%;
    margin: 0 0 12px;
    font-family: system-ui, sans-serif;
    font-size: 0.85em;
    container-type: inline-size;
  }
  .lfr-spd-slider-grid > * {
    flex: 1 1 calc((100% - 20px) / 2);
    min-width: 0;
    margin: 0;
  }
  .lfr-spd-slider-grid input[type="number"] {
    width: 6.5rem !important;
  }
  @container (max-width: 480px) {
    .lfr-spd-slider-grid > * {
      flex-basis: 100%;
    }
  }
</style>
<div class="lfr-spd-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 ----
function symTrace(m) { return m[0] + m[2]; }
function symDet(m) { return m[0]*m[2] - m[1]*m[1]; }

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];
}

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];
  var vx = b, vy = evals[0] - a;
  var vn = Math.sqrt(vx*vx + vy*vy);
  if (vn < 1e-14) return [Math.exp(a), 0, Math.exp(c)];
  vx /= vn; vy /= vn;
  return [
    e1*vx*vx + e2*vy*vy,
    e1*vx*(-vy) + e2*vy*vx,
    e1*vy*vy + e2*vx*vx
  ];
}

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];
  }
  if (Math.abs(wSum) < 1e-12) {
    // Degenerate: return unweighted mean
    return [a/n, b/n, c/n];
  }
  return [a/wSum, b/wSum, c/wSum];
}

// ---- True log-matrix as nonlinear function of x ----
// Uses a sigmoid-like transition to highlight boundary effects
function trueLogMatrix(x) {
  var a = 0.5 + 1.5*x + 0.6*Math.sin(2.5*Math.PI*x);
  var b = -0.3 + 0.8*x + 0.4*Math.cos(2.0*Math.PI*x);
  var c = 0.8 + 1.0*x + 0.5*Math.sin(1.8*Math.PI*x + 0.5);
  return [a, b, c];
}

// ---- Data generation ----
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();
    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 };
}

// ---- 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; });

trueLogs_app = xGrid_app.map(function(x) { return trueLogMatrix(x); });
trueTraces = trueLogs_app.map(symTrace);
trueDets = trueLogs_app.map(symDet);

// ---- Global Fréchet (Lecture 6) ----
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);
});

// ---- Kernel Fréchet (Lecture 7) ----
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);
});

// ---- Local-linear Fréchet (Lecture 8) ----
localLinLogs = 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);
  // Local moments
  var nu0 = 0, nu1 = 0, nu2 = 0;
  for (var i = 0; i < n_app; i++) {
    var di = appData.X[i] - x;
    nu0 += rawW[i];
    nu1 += rawW[i] * di;
    nu2 += rawW[i] * di * di;
  }
  nu0 /= n_app; nu1 /= n_app; nu2 /= n_app;
  var sigma02 = nu0 * nu2 - nu1 * nu1;
  var weights;
  if (Math.abs(sigma02) > 1e-10) {
    weights = appData.X.map(function(xi, i) { return rawW[i] * (nu2 - nu1*(xi - x)) / sigma02; });
  } else {
    weights = wSum > 1e-10 ? rawW.map(function(w) { return w/wSum; }) : rawW.map(function() { return 1/n_app; });
  }
  return symWeightedMean(appData.logs, weights);
});

// ---- Scalar summaries ----
globalTraces = globalLogs.map(symTrace); globalDets = globalLogs.map(symDet);
kernelTraces = kernelLogs.map(symTrace); kernelDets = kernelLogs.map(symDet);
llTraces = localLinLogs.map(symTrace);   llDets = localLinLogs.map(symDet);
dataTraces = appData.logs.map(symTrace); dataDets = appData.logs.map(symDet);

// Unconditional Fréchet mean
unifW = appData.logs.map(function() { return 1; });
meanLog = symWeightedMean(appData.logs, unifW);
meanSPD = symExp(meanLog);

// Local-linear weight diagnostics at left and right boundaries
function computeLLWeightsAt(xEval) {
  var rawW = appData.X.map(function(xi) { return gaussKernel((xi - xEval)/h_app); });
  var nu0 = 0, nu1 = 0, nu2 = 0;
  for (var i = 0; i < n_app; i++) {
    var di = appData.X[i] - xEval;
    nu0 += rawW[i]; nu1 += rawW[i]*di; nu2 += rawW[i]*di*di;
  }
  nu0 /= n_app; nu1 /= n_app; nu2 /= n_app;
  var sigma02 = nu0*nu2 - nu1*nu1;
  if (Math.abs(sigma02) > 1e-10) {
    return appData.X.map(function(xi, i) { return rawW[i] * (nu2 - nu1*(xi - xEval)) / sigma02; });
  }
  var wSum = rawW.reduce(function(a,b) { return a+b; }, 0);
  return rawW.map(function(w) { return w/wSum; });
}
llW_left = computeLLWeightsAt(0.05);
llW_right = computeLLWeightsAt(0.95);
negLeft = llW_left.filter(function(w) { return w < 0; }).length;
negRight = llW_right.filter(function(w) { return w < 0; }).length;

// ---- Visualization ----
function buildAppPlot(title, yLabel, trueVals, globalVals, kernelVals, llVals, dataVals, xVals, xGrid, yMin, yMax) {
  var W = 500, H = 210;
  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"/>');

  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.4"/>');
  });

  // True (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 (red, thinner)
  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="1.6" opacity="0.65"/>');

  // Kernel (orange)
  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="#f08c00" stroke-width="2.0"/>');

  // Local-linear (blue, thicker)
  d = "";
  for (var j = 0; j < xGrid.length; j++) d += (j===0?'M':'L') + " " + xS(xGrid[j]).toFixed(1) + " " + yS(llVals[j]).toFixed(1);
  parts.push('<path d="' + d + '" fill="none" stroke="#1971c2" stroke-width="2.6"/>');

  // Boundary shading (left and right 10% of predictor range)
  parts.push('<rect x="' + margin.left + '" y="' + margin.top + '" width="' + (0.1*pw) +
    '" height="' + ph + '" fill="#795548" opacity="0.05"/>');
  parts.push('<rect x="' + (margin.left + 0.9*pw) + '" y="' + margin.top + '" width="' + (0.1*pw) +
    '" height="' + ph + '" fill="#795548" opacity="0.05"/>');

  // 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>');

  // Boundary labels
  parts.push('<text x="' + (margin.left + 0.05*pw) + '" y="' + (margin.top+ph-4) +
    '" text-anchor="middle" font-size="7" fill="#795548" opacity="0.7">boundary</text>');
  parts.push('<text x="' + (margin.left + 0.95*pw) + '" y="' + (margin.top+ph-4) +
    '" text-anchor="middle" font-size="7" fill="#795548" opacity="0.7">boundary</text>');

  // Legend
  var lx = margin.left + pw - 170, ly = margin.top + 3;
  parts.push('<line x1="' + lx + '" y1="' + ly + '" x2="' + (lx+12) + '" y2="' + ly +
    '" stroke="#2e7d32" stroke-width="1.6" stroke-dasharray="4,3"/>');
  parts.push('<text x="' + (lx+15) + '" y="' + (ly+4) + '" font-size="7.5" fill="#616161">True</text>');
  parts.push('<line x1="' + (lx+42) + '" y1="' + ly + '" x2="' + (lx+54) + '" y2="' + ly +
    '" stroke="#e53935" stroke-width="1.6"/>');
  parts.push('<text x="' + (lx+57) + '" y="' + (ly+4) + '" font-size="7.5" fill="#616161">Glob</text>');
  parts.push('<line x1="' + (lx+80) + '" y1="' + ly + '" x2="' + (lx+92) + '" y2="' + ly +
    '" stroke="#f08c00" stroke-width="1.8"/>');
  parts.push('<text x="' + (lx+95) + '" y="' + (ly+4) + '" font-size="7.5" fill="#616161">Kern</text>');
  parts.push('<line x1="' + (lx+120) + '" y1="' + ly + '" x2="' + (lx+132) + '" y2="' + ly +
    '" stroke="#1971c2" stroke-width="2.2"/>');
  parts.push('<text x="' + (lx+135) + '" y="' + (ly+4) + '" font-size="7.5" fill="#616161">LL</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:700px;">\n' +
    parts.join('\n') + '\n</svg>';
}

// Compute y ranges
allTraceVals = dataTraces.concat(trueTraces, globalTraces, kernelTraces, llTraces);
traceMin = Math.min.apply(null, allTraceVals);
traceMax = Math.max.apply(null, allTraceVals);
tracePad = (traceMax - traceMin) * 0.12;
traceMinAdj = traceMin - tracePad; traceMaxAdj = traceMax + tracePad;

allDetVals = dataDets.concat(trueDets, globalDets, kernelDets, llDets);
detMin = Math.min.apply(null, allDetVals);
detMax = Math.max.apply(null, allDetVals);
detPad = (detMax - detMin) * 0.12;
detMinAdj = detMin - detPad; detMaxAdj = detMax + detPad;

plot_spd_view = html`
<div style="width:80%; max-width:700px; margin:0 auto; font-family:system-ui,sans-serif;">
  ${buildAppPlot("Trace of log Σ(x) vs. Stress", "tr(log Σ)", trueTraces, globalTraces, kernelTraces, llTraces, dataTraces, appData.X, xGrid_app, traceMinAdj, traceMaxAdj)}
  <div style="margin-top:10px;">
  ${buildAppPlot("Determinant of log Σ(x) vs. Stress", "det(log Σ)", trueDets, globalDets, kernelDets, llDets, dataDets, appData.X, xGrid_app, detMinAdj, detMaxAdj)}
  </div>
  <p style="margin-top:12px; font-size:0.85em; color:#495057;">
    <b>Observe the boundary regions</b> (shaded at x ≈ 0 and x ≈ 1):
    Local-linear (blue) tracks the true curve at the edges where kernel (orange) is pulled inward.
    Global (red) misses the nonlinearity entirely. Bandwidth h = ${h_app.toFixed(2)}.
  </p>
  <div style="background:#f8f9fa; border-radius:6px; padding:10px 14px; margin-top:8px; font-size:0.85em;">
    <b>Local-linear weight diagnostics at boundaries:</b>
    <table style="width:100%; border-collapse:collapse; margin-top:4px;">
      <tr><td style="padding:2px 8px;">At left boundary (x = 0.05)</td>
        <td style="padding:2px 8px; text-align:right;">${negLeft} of ${n_app} weights negative (${(negLeft/n_app*100).toFixed(1)}%)</td></tr>
      <tr><td style="padding:2px 8px;">At right boundary (x = 0.95)</td>
        <td style="padding:2px 8px; text-align:right;">${negRight} of ${n_app} weights negative (${(negRight/n_app*100).toFixed(1)}%)</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_spd_view}${plot_spd_view}</div>`

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

(u)

(v)

(w)

(x)

(y)

(z)

({)

(|)

(})

(~)

()

Figure 2: Interactive: Global vs. Kernel vs. Local-Linear Fréchet regression on simulated 2×2 SPD matrices

Try these experiments

Focus on the shaded boundary regions ($x$ near 0 and $x$ near 1): The kernel fit (orange) consistently lies inside the true curve — this is the $O(h)$ boundary bias. The local-linear fit (blue) stays closer to the true curve.
Decrease $h$ (0.05–0.10): The boundary bias of kernel regression becomes more localized but more severe. The local-linear weights become more negative (check the diagnostics panel) as the slope correction becomes more aggressive.
Increase $h$ (0.4–0.6): All three methods converge toward the unconditional Fréchet mean. The boundary advantage of local-linear diminishes as the local neighborhood expands to include the whole dataset.
Increase $n$ (150–200): The local-linear advantage at boundaries persists even with large samples — boundary bias is a bias problem, not a variance problem, so it doesn’t vanish with more data.
Compare trace and determinant: In some datasets, boundary bias may be more visible in one scalar summary than the other — this reflects which aspects of the covariance structure are most sensitive to the local slope at the edge.

9.5 Takeaway

Local-linear Fréchet regression inherits the computational simplicity of kernel regression (weighted Fréchet means on the SPD manifold) while correcting its primary weakness: $O(h)$ boundary bias. In the market risk context, where the boundaries correspond to extreme market conditions, this correction is not a theoretical nicety — it directly affects the accuracy of stress-tested portfolio risk estimates.

10 Key Takeaways

Local Fréchet regression extends local linear regression to metric spaces by constructing a predictor-side weight function $s(z, x, h)$ that automatically adjusts for the local slope.
The weights can be negative — this is the mechanism that corrects boundary bias.
The convergence rate $O_p(h^{2/(\beta_1-1)} + (nh)^{-1/[2(\beta_2-1)]})$ generalizes the classical bias–variance decomposition, with the exponents $\beta_1, \beta_2$ determined by the local curvature of the Fréchet objective.
For quadratic curvature ($\beta_1 = \beta_2 = 2$), the optimal bandwidth is $h \asymp n^{-1/5}$, yielding the familiar nonparametric rate $n^{-2/5}$.
Local-linear is almost always preferred over local-constant in practice — the computational cost is nearly identical, but the boundary behavior is substantially better.

11 Exercises

Local-linear weights sum to one. Show that $\sum_{i=1}^n s_{in}(x, h) = 1$. Also show that in the Euclidean case, $\sum_i s_{in}(x, h) X_i = x$ when the model includes an intercept and slope. 📝 Show Solution
Optimal bandwidth derivation. Derive the optimal bandwidth $h \asymp n^{-1/5}$ for the case $\beta_1 = \beta_2 = 2$ by balancing the squared bias and variance terms. 📝 Show Solution
Higher dimensions. Discuss the challenges of extending local Fréchet regression to $p > 1$ predictors. What is the bias rate in terms of $h$ and $p$? 📝 Show Solution
Computation. Discuss how to implement local Fréchet regression in general metric spaces and manifolds. Note that the weights could be negative, and thus the weighted proximal point algorithm may not converge in general.

12 Self-Assessment Quiz

Test your understanding of this lecture with the interactive MCQ quiz:

👉 Lecture 8 Quiz — 10 Multiple-Choice Questions

13 Further Reading

Petersen and Müller (2019) — Sections 4 and S.3: the full theory of local Fréchet regression.

References

Petersen, A., and H.-G. Müller. 2019. “Fréchet Regression for Random Objects with Euclidean Predictors.” The Annals of Statistics 47 (2): 691–719.

--- title: "Lecture 8: Fréchet Regression — Local Fréchet Regression" subtitle: "Local linear smoothing in metric spaces" 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 local-constant kernel regression suffers from boundary bias and how local linear regression corrects it. - Define the local moments $\nu_j(x,h)$ and the local weight function $s(z, x, h)$. - Derive the weight function $s(z, x, h)$ from the normal equations of Euclidean local linear regression. - Explain why $\mathbb{E}[s(X, x, h)] = 1$ and $\mathbb{E}[s(X, x, h)(X - x)] = 0$ are the two key properties enabling the metric-space generalization. - State the local Fréchet regression target and sample estimator. - Interpret the bias, stochastic error, and overall convergence rate of the local Fréchet regression estimator. - Derive the optimal bandwidth $h \asymp n^{-1/5}$ for the quadratic-curvature case $\beta_1 = \beta_2 = 2$. - Contrast local-constant (kernel) and local-linear Fréchet regression. ## Why Local Linear? Local-constant (kernel) Fréchet regression, discussed in Lecture 7, approximates the regression function as **constant** within each kernel window. This leads to: - $O(h^2)$ bias in the interior (acceptable). - $O(h)$ bias near **boundaries** of the predictor support (problematic). Local linear regression fixes this by fitting a **line** within each window, reducing the bias to $O(h^2)$ everywhere — including at boundaries. @PetersenMueller2019 extend this idea to metric spaces. For exposition, we consider the case $p = 1$ (scalar predictor). ## Local Moments and Weight Functions Let $K$ be a kernel density and $K_h(u) = h^{-1}K(u/h)$ for bandwidth $h > 0$. Define the **local moments** $$ \nu_j(x, h) = \mathbb{E}\bigl\{K_h(X - x)(X - x)^j\bigr\}, \qquad j = 0, 1, 2, $$ and $$ \sigma_0^2(x, h) = \nu_0(x, h) \nu_2(x, h) - \nu_1(x, h)^2. $$ First consider the ordinary Euclidean case $Y \in \mathbb{R}$. The local linear estimator at $x$ is the fitted intercept $\hat{\mu}(x) = \hat{\beta}_0$, where $(\hat{\beta}_0, \hat{\beta}_1)$ minimizes $$ \frac{1}{n}\sum_{i=1}^n K_h(X_i - x)\bigl\{Y_i - \beta_0 - \beta_1(X_i - x)\bigr\}^2. $$ At the population level, the corresponding target $(\bar{\beta}_0, \bar{\beta}_1)$ minimizes $$ \mathbb{E}\Bigl[K_h(X - x)\bigl\{Y - \beta_0 - \beta_1(X - x)\bigr\}^2\Bigr]. $$ ## From Euclidean Local Linear to the Weight Function {#sec-euclidean-to-weight} The key algebraic step — and the one that makes the metric-space generalization possible — is to rewrite the local linear intercept as a *weighted expectation* with a data-dependent weight function $s(z, x, h)$. We now derive this transformation in detail. ### Step 1: The Normal Equations Define the cross-moments $$ r_j(x, h) = \mathbb{E}\bigl\{K_h(X - x)(X - x)^j Y\bigr\}, \qquad j = 0, 1. $$ At the population level, $(\bar{\beta}_0, \bar{\beta}_1)$ minimizes $$ Q(\beta_0, \beta_1) = \mathbb{E}\Bigl[K_h(X - x)\bigl\{Y - \beta_0 - \beta_1(X - x)\bigr\}^2\Bigr]. $$ Differentiating with respect to $\beta_0$ and $\beta_1$ and setting to zero yields the normal equations: $$ \begin{aligned} \mathbb{E}\bigl[K_h(X - x)\{Y - \bar{\beta}_0 - \bar{\beta}_1(X - x)\}\bigr] &= 0, \\ \mathbb{E}\bigl[K_h(X - x)(X - x)\{Y - \bar{\beta}_0 - \bar{\beta}_1(X - x)\}\bigr] &= 0. \end{aligned} $$ Expanding and using the definitions of $\nu_j$ and $r_j$: $$ \boxed{ \begin{aligned} \nu_0(x, h)\,\bar{\beta}_0 + \nu_1(x, h)\,\bar{\beta}_1 &= r_0(x, h), \\ \nu_1(x, h)\,\bar{\beta}_0 + \nu_2(x, h)\,\bar{\beta}_1 &= r_1(x, h). \end{aligned}} \tag{1} $$ ### Step 2: Solving for the Intercept The system (1) is a $2 \times 2$ linear system whose coefficient matrix has determinant $\sigma_0^2 = \nu_0\nu_2 - \nu_1^2 > 0$ (strictly positive since $K \ge 0$ and $X$ has a continuous density). By Cramer's rule, the intercept is $$ \tilde{\mu}(x) = \bar{\beta}_0 = \frac{\det\begin{pmatrix} r_0 & \nu_1 \\ r_1 & \nu_2 \end{pmatrix}}{\det\begin{pmatrix} \nu_0 & \nu_1 \\ \nu_1 & \nu_2 \end{pmatrix}} = \frac{\nu_2(x, h)\, r_0(x, h) - \nu_1(x, h)\, r_1(x, h)}{\sigma_0^2(x, h)}. \tag{2} $$ ### Step 3: The Weight Function Emerges Now substitute the definition $r_j(x, h) = \mathbb{E}[K_h(X - x)(X - x)^j Y]$ into (2): $$ \begin{aligned} \tilde{\mu}(x) &= \frac{\nu_2\,\mathbb{E}[K_h(X - x)Y] - \nu_1\,\mathbb{E}[K_h(X - x)(X - x)Y]}{\sigma_0^2} \\[4pt] &= \mathbb{E}\!\left[\frac{K_h(X - x)\bigl\{\nu_2(x, h) - \nu_1(x, h)(X - x)\bigr\}}{\sigma_0^2(x, h)}\, Y\right]. \end{aligned} $$ This is of the form $\tilde{\mu}(x) = \mathbb{E}[s(X, x, h)\,Y]$, where the **weight function** is $$ \boxed{s(z, x, h) = \frac{K_h(z - x)\bigl\{\nu_2(x, h) - \nu_1(x, h)(z - x)\bigr\}}{\sigma_0^2(x, h)}}. \tag{3} $$ ### Step 4: Properties of the Weight Function Two properties of $s(z, x, h)$ are essential for what follows. **Sum-to-one.** Take $Y \equiv 1$ in the derivation above. Then $r_0 = \nu_0$ and $r_1 = \nu_1$, so (2) gives $\tilde{\mu}(x) = 1$. But $\tilde{\mu}(x) = \mathbb{E}[s(X, x, h) \cdot 1]$, hence $$ \boxed{\mathbb{E}[s(X, x, h)] = 1}. \tag{4} $$ **First-moment annihilation.** Take $Y = X - x$. Then $r_0 = \nu_1$ and $r_1 = \nu_2$, so (2) gives $\tilde{\mu}(x) = 0$. But $\tilde{\mu}(x) = \mathbb{E}[s(X, x, h)(X - x)]$, hence $$ \boxed{\mathbb{E}[s(X, x, h)(X - x)] = 0}. \tag{5} $$ Equation (5) is the mathematical expression of the bias-correction mechanism: the weights are orthogonal to linear trends, so the local-linear estimator is **exactly unbiased** for linear conditional mean functions — a property the local-constant (kernel) estimator lacks. ### Step 5: Reformulation as Weighted Least Squares Because $\mathbb{E}[s(X, x, h)] = 1$, we can add a constant term (with respect to $y$) without changing the minimizer: $$ \begin{aligned} \tilde{\mu}(x) &= \mathbb{E}[s(X, x, h) Y] \\ &= \arg\min_{y \in \mathbb{R}} \Bigl\{\mathbb{E}[s(X, x, h) Y^2] - 2y\,\mathbb{E}[s(X, x, h) Y] + y^2\Bigr\} \\ &= \arg\min_{y \in \mathbb{R}} \mathbb{E}\bigl[s(X, x, h)(Y - y)^2\bigr]. \end{aligned} $$ Thus $$ \boxed{\tilde{\mu}(x) = \arg\min_{y \in \mathbb{R}} \mathbb{E}\bigl\{s(X, x, h)(Y - y)^2\bigr\}}. \tag{6} $$ ::: {.callout-tip title="Why this form matters"} Equation (6) expresses the local linear target purely in terms of **a weighted expectation of a squared distance**. This is the form that generalizes: the weight function $s(z, x, h)$ is built entirely from the predictor $X$ and depends on the response $Y$ only through the squared-error criterion. To extend to a metric-space-valued response, we simply: 1. Replace the squared Euclidean distance $(Y - y)^2$ by the squared metric $d^2(Y, \omega)$, 2. Minimize over $\omega$ in the metric space $\mathcal{M}$ rather than over $y \in \mathbb{R}$. The weight function $s(z, x, h)$ remains **unchanged** — it is computed from the scalar predictors exactly as in the Euclidean case. ::: ## The Local Fréchet Regression Target and Estimator Now consider data $(X_i, Y_i) \in \mathbb{R} \times \mathcal{M}$, $i = 1, \dots, n$, where $\mathcal{M}$ is a metric space with distance $d$. ::: {#def-local-frechet .definition title="Local Fréchet regression target"} Let $$ \widetilde{L}_h(\omega, x) = \mathbb{E}\{s(X, x, h) d^2(Y, \omega)\}, $$ The **local Fréchet regression target** is $$ \tilde{\mu}(x) = \arg\min_{\omega \in \mathcal{M}} \widetilde{L}_h(\omega, x). $$ ::: This target is not exactly the conditional Fréchet mean $\mu(x)$; it is the metric-space analogue of the usual local linear population target. The key fact proved by Petersen and Müller is that the induced weighted measure differs from the true conditional law by $O(h^2)$, so minimizing $\widetilde{L}_h(\cdot, x)$ approximates minimizing the conditional Fréchet function. Given data, define empirical local moments $$ \hat{\nu}_j(x, h) = \frac{1}{n}\sum_{i=1}^n K_h(X_i - x)(X_i - x)^j, \qquad j = 0, 1, 2, $$ $$ \hat{\sigma}_0^2(x, h) = \hat{\nu}_0(x, h) \hat{\nu}_2(x, h) - \hat{\nu}_1(x, h)^2, $$ and weights $$ s_{in}(x, h) = \frac{K_h(X_i - x)\{\hat{\nu}_2(x, h) - \hat{\nu}_1(x, h)(X_i - x)\}}{\hat{\sigma}_0^2(x, h)}. $$ The sample local objective is $$ \hat{L}_n(\omega, x) = \frac{1}{n}\sum_{i=1}^n s_{in}(x, h) d^2(Y_i, \omega), $$ and the **local Fréchet regression estimator** is $$ \hat{\mu}(x) = \arg\min_{\omega \in \mathcal{M}} \hat{L}_n(\omega, x). $$ ::: {.callout-note title="Important distinction"} Unlike kernel regression weights, the local linear weights $s_{in}(x, h)$ can be **negative**. This is the price for the bias correction — the weights automatically adjust for the local slope. The empirical weights satisfy $\sum_i s_{in}(x, h) = 1$ and $\sum_i s_{in}(x, h)(X_i - x) = 0$. ::: ## Bias, Stochastic Error, and Convergence Rates The following regularity assumptions are from @PetersenMueller2019, Section 4: ::: {.callout-note title="Assumptions (K0), (L0)–(L3), and (P1)"} **(K0) Kernel regularity.** $K$ is a symmetric probability density with finite moments $|K_{14}| < \infty$ and $|K_{26}| < \infty$, where $K_{kj} = \int_{\mathbb{R}} K(u)^k u^j\, du$. **(L0) Existence and separation.** $\mu(x)$ and $\tilde{\mu}(x)$ exist and are unique. Also, $\hat{\mu}(x)$ exists and is unique, almost surely. For any $\epsilon > 0$, $$ \inf_{d(\omega, \mu(x)) > \epsilon}\{F(\omega, x) - F(\mu(x), x)\} > 0, \qquad \liminf_{n \to \infty} \inf_{d(\omega, \tilde{\mu}(x)) > \epsilon}\{\widetilde{L}_n(\omega, x) - \widetilde{L}_n(\tilde{\mu}(x), x)\} > 0. $$ **(L1) Smooth predictor law.** The marginal density $f$ of $X$ and the conditional densities $g_y$ of $X \mid Y = y$ exist and are twice continuously differentiable, with $\sup_{x, y} |g_y''(x)| < \infty$. **(L2) Bias curvature.** There exist $\eta_1 > 0$, $C_1 > 0$, and $\beta_1 > 1$ such that $$ F(\omega, x) - F(\mu(x), x) \ge C_1 d(\omega, \mu(x))^{\beta_1} \quad \text{whenever } d(\omega, \mu(x)) < \eta_1. $$ **(L3) Stochastic curvature.** There exist $\eta_2 > 0$, $C_2 > 0$, and $\beta_2 > 1$ such that $$ \liminf_{n \to \infty} \bigl[\widetilde{L}_n(\omega, x) - \widetilde{L}_n(\tilde{\mu}(x), x)\bigr] \ge C_2 d(\omega, \tilde{\mu}(x))^{\beta_2} \quad \text{whenever } d(\omega, \tilde{\mu}(x)) < \eta_2. $$ **(P1) Local entropy.** Let $B_\delta(\mu(x)) \subset \mathcal{M}$ be the ball of radius $\delta$ centered at $\mu(x)$, and let $N(\varepsilon, B_\delta(\mu(x)), d)$ denote its covering number by balls of radius $\varepsilon$. Then $$ \int_0^1 \bigl\{1 + \log N(\delta\varepsilon, B_\delta(\mu(x)), d)\bigr\}^{1/2}\, d\varepsilon = O(1) \qquad \text{as } \delta \downarrow 0. $$ This controls the local complexity of the metric space near the target, ensuring that the stochastic fluctuation of the empirical objective on shrinking neighborhoods of $\mu(x)$ remains bounded. For finite-dimensional manifolds (e.g., $\mathbb{S}^d$, $\mathcal{S}_{++}^m$), the covering number grows as $\varepsilon^{-\dim}$, and the integral condition is automatically satisfied. For infinite-dimensional function spaces, it may fail and requires verification. ::: ::: {#thm-local-frechet .theorem title="Local Fréchet regression: bias, stochastic error, and convergence rate (Theorems 3–4, Corollary 1 of Petersen & Müller 2019)"} Suppose $(\mathcal{M}, d)$ is totally bounded. - **Bias:** If (K0), (L0), (L1), (L2), and (P1) hold, then as $h = h_n \to 0$, $$ d(\mu(x), \tilde{\mu}(x)) = O\bigl(h^{2/(\beta_1 - 1)}\bigr). $$ - **Stochastic error:** If (K0), (L0), (L3), and (P1) hold, and $h \to 0$ with $nh \to \infty$, then $$ d(\tilde{\mu}(x), \hat{\mu}(x)) = O_p\bigl((nh)^{-1/[2(\beta_2 - 1)]}\bigr). $$ - **Overall rate:** Under all assumptions, $$ d(\mu(x), \hat{\mu}(x)) = O_p\bigl(h^{2/(\beta_1 - 1)} + (nh)^{-1/[2(\beta_2 - 1)]}\bigr). $$ Choosing $h = n^{-\gamma}$ with $\gamma^* = (\beta_1 - 1)/(4\beta_2 + \beta_1 - 5)$ yields the optimal rate $$ d(\mu(x), \hat{\mu}(x)) = O_p\bigl(n^{-2/(\beta_1 + 4\beta_2 - 5)}\bigr). $$ ::: ::: {.callout-important title="The quadratic-curvature case $\beta_1 = \beta_2 = 2$"} In the important special case where both curvature exponents are 2 (which holds for SPD spaces, BHV space, and Euclidean spaces): $$ d(\mu(x), \hat{\mu}(x)) = O_p\bigl(h^2 + (nh)^{-1/2}\bigr), $$ and the optimal bandwidth $h \asymp n^{-1/5}$ yields the familiar nonparametric rate $n^{-2/5}$. This shows that Fréchet regression **inherits the same bias–variance tradeoff** as classical local linear regression, but with the exponents $\beta_1$ and $\beta_2$ determined by the local geometry of the response space. ::: ## Comparison: Local-Constant vs. Local-Linear | Property | Local-Constant (Kernel) | Local-Linear | |---|---|---| | **Bias (interior)** | $O(h^2)$ | $O(h^2)$ | | **Bias (boundary)** | $O(h)$ | $O(h^2)$ | | **Weights** | Always nonnegative | Can be negative | | **Finite-sample** | Simpler, more stable | Better at boundaries | | **Computational** | Weighted Fréchet mean | Weighted Fréchet mean (same complexity) | ## Interactive Exploration: Global vs. Kernel vs. Local-Linear on the Sphere The following demo compares all three Fréchet regression methods on the positive quadrant (octant) of the unit sphere $\mathbb{S}^2$: - **Global Fréchet regression** (Lecture 6) — the metric analogue of a global linear model. - **Kernel (local-constant) Fréchet regression** (Lecture 7) — the metric analogue of Nadaraya–Watson smoothing. - **Local-linear Fréchet regression** (this lecture) — the metric analogue of local linear regression. The true regression curve is **nonlinear** (sinusoidal), so the global fit is misspecified. The kernel fit captures the nonlinearity but suffers from boundary bias. The local-linear fit corrects the boundary bias by adjusting weights to account for the local slope — some weights become **negative** as a result. **Data generation.** Same as Lectures 6–7: $X_i \sim \text{Uniform}(-1, 1)$, true $\mu(x)$ is a sinusoidal curve on the sphere octant, and $Y_i$ is generated with geodesic noise. **Three weight functions compared:** | Method | Weights $w_i(x)$ | Can be negative? | |---|---|---| | Global (L6) | $1 + \frac{(X_i - \bar{X})(x - \bar{X})}{\hat{\sigma}_X^2}$ | Yes | | Kernel (L7) | $\frac{K_h(X_i - x)}{\sum_j K_h(X_j - x)}$ | No | | Local-linear (L8) | $\frac{K_h(X_i - x)\{\hat{\nu}_2 - \hat{\nu}_1(X_i - x)\}}{\hat{\nu}_0\hat{\nu}_2 - \hat{\nu}_1^2}$ | Yes | All three then compute the weighted Fréchet mean on $\mathbb{S}^2$ via Riemannian gradient descent. **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 (L6) - **Orange curve** = kernel / local-constant fit (L7) - **Blue curve** = local-linear fit (L8) - **Gold diamond** = unconditional Fréchet mean of all $Y_i$ ```{ojs} //| label: fig-local-frechet-sphere //| fig-cap: "Interactive: Global vs. Kernel vs. Local-Linear 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_sphere_view = html` <style> .lfr-control-panel { width: 100%; margin: 0 0 12px; font-family: system-ui, sans-serif; font-size: 0.85em; container-type: inline-size; } .lfr-control-panel > *, .lfr-slider-grid > * { min-width: 0; margin: 0; } .lfr-slider-grid { display: flex; flex-wrap: wrap; gap: 6px 20px; margin-top: 6px; } .lfr-slider-grid > * { flex: 1 1 calc((100% - 40px) / 3); } .lfr-slider-grid input[type="number"] { width: 7.5rem !important; } @container (max-width: 700px) { .lfr-slider-grid > * { flex-basis: calc((100% - 20px) / 2); } } @container (max-width: 480px) { .lfr-slider-grid > * { flex-basis: 100%; } } </style> <div class="lfr-control-panel"> <div>${reg_type_control}</div> <div class="lfr-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 lectures 6–7) ---- 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) 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 }; } // ---- Gaussian kernel ---- function gaussKernel(u) { return Math.exp(-0.5 * u * u); } // ---- 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 (local-constant) Fréchet regression (Lecture 7) ---- 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; } // ---- Local-linear Fréchet regression (Lecture 8) ---- function computeLocalLinearCurve(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 = []; for (var i = 0; i < n; i++) { rawW.push(gaussKernel((X[i] - x) / h)); } // Compute local moments var nu0 = 0, nu1 = 0, nu2 = 0; for (var i = 0; i < n; i++) { var di = X[i] - x; nu0 += rawW[i]; nu1 += rawW[i] * di; nu2 += rawW[i] * di * di; } nu0 /= n; nu1 /= n; nu2 /= n; var sigma02 = nu0 * nu2 - nu1 * nu1; var weights; if (Math.abs(sigma02) > 1e-10) { // Local-linear weights (can be negative — the bias-correction mechanism) weights = X.map(function(xi, i) { return rawW[i] * (nu2 - nu1 * (xi - x)) / sigma02; }); } else { // Fallback to kernel (local-constant) var wSum = rawW.reduce(function(a,b) { return a+b; }, 0); 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); localLinCurve = computeLocalLinearCurve(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); // Local-linear at eval nu0e = rawEvalW.reduce(function(a,b) { return a+b; }, 0) / regData.X.length; nu1e = rawEvalW.reduce(function(s, w, i) { return s + w * (regData.X[i] - eval_x); }, 0) / regData.X.length; nu2e = rawEvalW.reduce(function(s, w, i) { return s + w * Math.pow(regData.X[i] - eval_x, 2); }, 0) / regData.X.length; sigma02e = nu0e * nu2e - nu1e * nu1e; llWeightsEval = (Math.abs(sigma02e) > 1e-10) ? regData.X.map(function(xi, i) { return rawEvalW[i] * (nu2e - nu1e * (xi - eval_x)) / sigma02e; }) : nWeightsEval; localLinAtEval = weightedSphereFrechetMean(regData.Y, llWeightsEval, 200, 0.3); trueAtEval = trueRegFunc(eval_x, reg_type); globalErr = sphereDist(globalAtEval, trueAtEval); kernelErr = sphereDist(kernelAtEval, trueAtEval); localLinErr = sphereDist(localLinAtEval, trueAtEval); // Overall Fréchet mean unifWeights = regData.Y.map(function() { return 1; }); overallMean = weightedSphereFrechetMean(regData.Y, unifWeights, 200, 0.3); // Diagnostics effN = rawEvalW.filter(function(w) { return w > 1e-6; }).length; llNegCount = llWeightsEval.filter(function(w) { return w < 0; }).length; llWeightSum = llWeightsEval.reduce(function(a,b) { return a+b; }, 0); // ---- 3D Visualization of the sphere octant ---- function build3DView(regData, trueCurve, globalCurve, kernelCurve, localLinCurve, xGridEval, globalAtEval, kernelAtEval, localLinAtEval, 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) 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, thinner) 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.2" stroke-linecap="round" stroke-linejoin="round" opacity="0.7"/>'); // Kernel / local-constant curve (orange) 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="#f08c00" stroke-width="2.6" stroke-linecap="round" stroke-linejoin="round" opacity="0.8"/>'); // Local-linear curve (blue, thicker — the featured method) var llPtsStr = localLinCurve.map(function(p) { var proj = project(p); return proj[0].toFixed(1) + "," + proj[1].toFixed(1); }).join(" "); parts.push('<polyline points="' + llPtsStr + '" fill="none" stroke="#1971c2" stroke-width="3.2" stroke-linecap="round" stroke-linejoin="round" opacity="0.9"/>'); // 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 lEvalProj = project(localLinAtEval); var tEvalProj = project(trueAtEval); // Dashed connectors from each fit 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.0" stroke-dasharray="3,2" opacity="0.45"/>'); 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.0" stroke-dasharray="3,2" opacity="0.45"/>'); parts.push('<line x1="' + lEvalProj[0].toFixed(1) + '" y1="' + lEvalProj[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="6" fill="none" stroke="#e53935" stroke-width="2.5" opacity="0.7"/>'); parts.push('<circle cx="' + kEvalProj[0].toFixed(1) + '" cy="' + kEvalProj[1].toFixed(1) + '" r="6" fill="none" stroke="#f08c00" stroke-width="2.5" opacity="0.8"/>'); parts.push('<circle cx="' + lEvalProj[0].toFixed(1) + '" cy="' + lEvalProj[1].toFixed(1) + '" r="7" fill="none" stroke="#1971c2" stroke-width="3.5"/>'); 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 = 35; 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+16) + '" x2="' + (lx+14) + '" y2="' + (ly0+16) + '" stroke="#2e7d32" stroke-width="2.5" stroke-dasharray="5,3"/>'); parts.push('<text x="' + (lx+22) + '" y="' + (ly0+20) + '" font-size="9" fill="#37474f">True μ(x)</text>'); parts.push('<line x1="' + (lx-6) + '" y1="' + (ly0+32) + '" x2="' + (lx+14) + '" y2="' + (ly0+32) + '" stroke="#e53935" stroke-width="2.0"/>'); parts.push('<text x="' + (lx+22) + '" y="' + (ly0+36) + '" font-size="9" fill="#37474f">Global (L6)</text>'); parts.push('<line x1="' + (lx-6) + '" y1="' + (ly0+48) + '" x2="' + (lx+14) + '" y2="' + (ly0+48) + '" stroke="#f08c00" stroke-width="2.2"/>'); parts.push('<text x="' + (lx+22) + '" y="' + (ly0+52) + '" font-size="9" fill="#37474f">Kernel (L7)</text>'); parts.push('<line x1="' + (lx-6) + '" y1="' + (ly0+64) + '" x2="' + (lx+14) + '" y2="' + (ly0+64) + '" stroke="#1971c2" stroke-width="2.8"/>'); parts.push('<text x="' + (lx+22) + '" y="' + (ly0+68) + '" font-size="9" fill="#37474f">Local-linear (L8)</text>'); parts.push('<polygon points="' + lx + ',' + (ly0+78) + ' ' + (lx+4) + ',' + (ly0+83) + ' ' + lx + ',' + (ly0+88) + ' ' + (lx-4) + ',' + (ly0+83) + '" fill="#f9a825"/>'); parts.push('<text x="' + (lx+12) + '" y="' + (ly0+87) + '" 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:700px;">\n' + parts.join('\n') + '\n</svg>'; } sphere_l8_view = html`<div style="width:80%; max-width:700px; margin:0 auto; text-align:center;">${build3DView(regData, trueCurve, globalCurve, kernelCurve, localLinCurve, xGridEval, globalAtEval, kernelAtEval, localLinAtEval, 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); }); localLinAngles = localLinCurve.map(function(p) { return sphericalAngles(p); }); dataAngles = regData.Y.map(function(yi) { return sphericalAngles(yi); }); evalTrueAng = sphericalAngles(trueAtEval); evalGlobalAng = sphericalAngles(globalAtEval); evalKernelAng = sphericalAngles(kernelAtEval); evalLocalLinAng = sphericalAngles(localLinAtEval); function buildUnwrappedPlot(title, yLabel, trueVals, globalVals, kernelVals, llVals, dataVals, xVals, xGrid, yMin, yMax) { var W = 620, 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"/>'); 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.4"/>'); }); // True curve (green, dashed) var line = ""; for (var j = 0; j < xGrid.length; j++) { line += (j === 0 ? 'M' : 'L') + " " + xToSvg(xGrid[j]).toFixed(1) + " " + yToSvg(trueVals[j]).toFixed(1); } parts.push('<path d="' + line + '" fill="none" stroke="#2e7d32" stroke-width="2.2" stroke-dasharray="7,4"/>'); // Global (red) line = ""; for (var j = 0; j < xGrid.length; j++) { line += (j === 0 ? 'M' : 'L') + " " + xToSvg(xGrid[j]).toFixed(1) + " " + yToSvg(globalVals[j]).toFixed(1); } parts.push('<path d="' + line + '" fill="none" stroke="#e53935" stroke-width="1.8" opacity="0.7"/>'); // Kernel (orange) line = ""; for (var j = 0; j < xGrid.length; j++) { line += (j === 0 ? 'M' : 'L') + " " + xToSvg(xGrid[j]).toFixed(1) + " " + yToSvg(kernelVals[j]).toFixed(1); } parts.push('<path d="' + line + '" fill="none" stroke="#f08c00" stroke-width="2.0"/>'); // Local-linear (blue, thicker) line = ""; for (var j = 0; j < xGrid.length; j++) { line += (j === 0 ? 'M' : 'L') + " " + xToSvg(xGrid[j]).toFixed(1) + " " + yToSvg(llVals[j]).toFixed(1); } parts.push('<path d="' + line + '" fill="none" stroke="#1971c2" stroke-width="2.4"/>'); // Evaluation point 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"/>'); 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>'); 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 - 178, ly = margin.top + 4; parts.push('<line x1="' + lx + '" y1="' + ly + '" x2="' + (lx+12) + '" y2="' + ly + '" stroke="#2e7d32" stroke-width="1.6" stroke-dasharray="4,3"/>'); parts.push('<text x="' + (lx+15) + '" y="' + (ly+4) + '" font-size="7.5" fill="#616161">True</text>'); parts.push('<line x1="' + (lx+44) + '" y1="' + ly + '" x2="' + (lx+56) + '" y2="' + ly + '" stroke="#e53935" stroke-width="1.6"/>'); parts.push('<text x="' + (lx+59) + '" y="' + (ly+4) + '" font-size="7.5" fill="#616161">Global</text>'); parts.push('<line x1="' + (lx+90) + '" y1="' + ly + '" x2="' + (lx+102) + '" y2="' + ly + '" stroke="#f08c00" stroke-width="1.8"/>'); parts.push('<text x="' + (lx+105) + '" y="' + (ly+4) + '" font-size="7.5" fill="#616161">Kernel</text>'); parts.push('<line x1="' + (lx+138) + '" y1="' + ly + '" x2="' + (lx+150) + '" y2="' + ly + '" stroke="#1971c2" stroke-width="2"/>'); parts.push('<text x="' + (lx+153) + '" y="' + (ly+4) + '" font-size="7.5" fill="#616161">LocalLin</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:700px;">\n' + parts.join('\n') + '\n</svg>'; } unwrapped_l8_view = html` <div style="width:80%; max-width:700px; 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; }), localLinAngles.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; }), localLinAngles.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_l8_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)}, h = ${bandwidth.toFixed(2)})</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 style="border-bottom:1px solid #dee2e6;"> <td style="padding:5px 8px;" colspan="2"><b>Local-linear (L8)</b></td> </tr> <tr><td style="padding:4px 8px;">Local-linear μ̂(x)</td> <td style="padding:4px 8px; text-align:right; color:#1971c2; font-weight:bold;">θ=${evalLocalLinAng.theta.toFixed(3)}, φ=${evalLocalLinAng.phi.toFixed(3)} rad</td></tr> <tr><td style="padding:4px 8px;">Local-linear error</td> <td style="padding:4px 8px; text-align:right; color:#1971c2; font-weight:bold;">${localLinErr.toFixed(5)} rad (${(localLinErr*180/Math.PI).toFixed(3)}°)</td></tr> <tr><td style="padding:4px 8px;">Negative weights</td> <td style="padding:4px 8px; text-align:right;">${llNegCount} of ${regData.X.length} (${(llNegCount/regData.X.length*100).toFixed(1)}%)</td></tr> <tr style="border-bottom:1px solid #dee2e6;"> <td style="padding:5px 8px;" colspan="2"><b>Kernel / local-constant (L7)</b></td> </tr> <tr><td style="padding:4px 8px;">Kernel μ̂(x)</td> <td style="padding:4px 8px; text-align:right; color:#f08c00;">θ=${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:#f08c00;">${kernelErr.toFixed(5)} rad (${(kernelErr*180/Math.PI).toFixed(3)}°)</td></tr> <tr style="border-bottom:1px solid #dee2e6;"> <td style="padding:5px 8px;" colspan="2"><b>Global (L6)</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;">${globalErr.toFixed(5)} rad (${(globalErr*180/Math.PI).toFixed(3)}°)</td></tr> <tr style="border-bottom:1px solid #dee2e6;"> <td style="padding:5px 8px;" colspan="2"><b>Common</b></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><td style="padding:4px 8px;">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_sphere_view}${sphere_l8_view}${unwrapped_l8_view}${statistics_l8_view}</div>` ``` ::: {.callout-tip title="Try these experiments"} - **Compare all three methods in the 3D view:** Global (red) is a straight-line analogue — it misses the sinusoidal oscillation entirely. Kernel (orange) tracks the curve but can be pulled toward the interior at boundaries. Local-linear (blue) tracks the curve everywhere, including near the edges. - **Toggle regression types:** "Varying θ only" and "Varying φ only" isolate the nonlinearity to one coordinate — observe which methods capture it best in the 2D unwrapped views. - **Look near the edges ($x = \pm 1$):** Local-linear weights become negative (reported in the stats panel) — this is the mechanism that corrects boundary bias. The kernel weights remain nonnegative and exhibit boundary bias as a result. - **Decrease $h$ (0.15–0.25):** All local fits become more wiggly. The local-linear advantage at boundaries becomes more visible. - **Increase $h$ (0.8–1.2):** Kernel and local-linear both approach the unconditional Fréchet mean (gold diamond). At very large $h$, all three methods nearly coincide at $\bar{X}$ (since they all pass through the Fréchet mean there). - **Increase $n$ (200+):** All fits improve, but the global fit's bias persists — no amount of data can fix model misspecification. - **Study the negative weight count:** Local-linear weights are negative near boundaries and in high-curvature regions. This is the mathematical mechanism that makes local-linear unbiased for linear functions — exactly as in classical local polynomial regression. ::: ## Application: Stress-Testing at the Boundary — Why Local-Linear Matters Lecture 6 gave us a **global linear** model of how market covariance changes with stress. Lecture 7 replaced the linearity assumption with **kernel smoothing**, capturing nonlinear relationships but introducing $O(h)$ bias at the edges of the predictor range. This lecture's local-linear method corrects that boundary bias — and in the market risk context, the boundaries are precisely where we care most. ### The Boundary Problem in Risk Management In the market risk application, the predictor is realized SPY volatility $X_t = \hat{\sigma}_t$, and the response is a rolling covariance matrix $\Sigma_t \in \mathcal{S}_{++}^p$. Stress-testing asks: *what is the typical covariance structure at an extreme stress level $x^*$ near or beyond the edge of observed data?* At such boundary points: - **Global Fréchet regression** (Lecture 6) extrapolates linearly — its weights grow with distance from $\bar{X}$, producing estimates that can be far from the local data. - **Kernel Fréchet regression** (Lecture 7) uses only nearby observations — but at the boundary, all nearby observations lie on one side of $x^*$. The kernel-weighted average is pulled *inward* toward the interior, creating $O(h)$ boundary bias. - **Local-linear Fréchet regression** (this lecture) corrects this by allowing weights to become **negative** for observations on the dense side, effectively extrapolating the local slope back to the boundary. The bias drops to $O(h^2)$ everywhere. In stress-testing, an $O(h)$ error in the estimated covariance matrix can translate to a material misestimate of portfolio risk at the very moment it matters most. ### Local-Linear Fréchet Regression on SPD Matrices Under the log-Euclidean metric, the local-linear estimator at stress level $x$ is $$ \log \hat{\mu}_{\text{LL}}(x) = \frac{1}{n}\sum_{i=1}^n s_{in}(x, h)\, \log \Sigma_i, \qquad s_{in}(x, h) = \frac{K_h(X_i - x)\{\hat{\nu}_2(x, h) - \hat{\nu}_1(x, h)(X_i - x)\}}{\hat{\nu}_0(x, h)\hat{\nu}_2(x, h) - \hat{\nu}_1(x, h)^2}, $$ where $\hat{\nu}_j(x, h) = \frac{1}{n}\sum_{i=1}^n K_h(X_i - x)(X_i - x)^j$ are the empirical local moments. As in the global and kernel cases, the computational core is simply a **weighted average of matrix logarithms** — the geometry is handled by the log-Euclidean framework, and the regression method only changes the weights. ::: {.callout-important title="Three methods, one computational primitive"} All three Fréchet regression methods share the same form under the log-Euclidean metric: $$\log \hat{\mu}(x) = \sum_{i=1}^n w_i(x)\, \log \Sigma_i, \qquad \sum_i w_i(x) = 1.$$ | Method | Weight $w_i(x)$ | Sign | Boundary bias | |---|---|---|---| | Global (L6) | $\frac{1}{n}\bigl[1 + \frac{(X_i - \bar{X})(x - \bar{X})}{\hat{\sigma}_X^2}\bigr]$ | Any | $O(1)$ if extrapolating | | Kernel (L7) | $\frac{K_h(X_i - x)}{\sum_j K_h(X_j - x)}$ | $\ge 0$ | $O(h)$ | | Local-linear (L8) | $\frac{K_h(X_i - x)\{\hat{\nu}_2 - \hat{\nu}_1(X_i - x)\}}{n(\hat{\nu}_0\hat{\nu}_2 - \hat{\nu}_1^2)}$ | Any | $O(h^2)$ | The only difference between methods is the weight function — a purely predictor-side calculation. Switching between them requires no new geometry code. ::: ### What the Local-Linear Correction Buys You At an extreme stress level $x^*$ (say, the 95th percentile of observed volatility): 1. **Kernel regression** sees mostly moderate-stress observations to the left of $x^*$. Its estimate is a weighted average of their covariance matrices — biased toward calmer-market structures. The resulting risk model *underestimates* the stress-typical correlations. 2. **Local-linear regression** sees the same data but fits a line through them. If covariances are trending upward with stress, the negative weights on moderate-stress observations pull the estimate *further* in the stress direction, extrapolating the trend. The risk model reflects what the data suggest about the boundary, not just what they show at the interior. 3. **Global regression** also extrapolates, but with a single straight line fitted to all data — it cannot adapt to local trends that may differ at the left and right boundaries. The practical consequence: for a portfolio weight vector $w$, the stress-conditional variance estimate $\sigma_P^2(x^*) = w^\top \hat{\mu}_{\text{LL}}(x^*) w$ from local-linear regression is the most trustworthy of the three at extreme stress levels, because it alone combines **local adaptivity** with **boundary debiasing**. ### Interactive Exploration: All Three Methods on Simulated SPD Data The following demo extends the application from Lecture 7 by adding the local-linear estimator. We simulate $2 \times 2$ SPD matrices with a nonlinear relationship to a scalar stress predictor, and fit all three methods. Pay attention to the **edges** ($x \approx 0$ and $x \approx 1$) where local-linear should outperform kernel regression. ```{ojs} //| label: fig-local-app-spd //| fig-cap: "Interactive: Global vs. Kernel vs. Local-Linear 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.18, 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_spd_view = html` <style> .lfr-spd-slider-grid { display: flex; flex-wrap: wrap; gap: 6px 20px; width: 100%; margin: 0 0 12px; font-family: system-ui, sans-serif; font-size: 0.85em; container-type: inline-size; } .lfr-spd-slider-grid > * { flex: 1 1 calc((100% - 20px) / 2); min-width: 0; margin: 0; } .lfr-spd-slider-grid input[type="number"] { width: 6.5rem !important; } @container (max-width: 480px) { .lfr-spd-slider-grid > * { flex-basis: 100%; } } </style> <div class="lfr-spd-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 ---- function symTrace(m) { return m[0] + m[2]; } function symDet(m) { return m[0]*m[2] - m[1]*m[1]; } 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]; } 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]; var vx = b, vy = evals[0] - a; var vn = Math.sqrt(vx*vx + vy*vy); if (vn < 1e-14) return [Math.exp(a), 0, Math.exp(c)]; vx /= vn; vy /= vn; return [ e1*vx*vx + e2*vy*vy, e1*vx*(-vy) + e2*vy*vx, e1*vy*vy + e2*vx*vx ]; } 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]; } if (Math.abs(wSum) < 1e-12) { // Degenerate: return unweighted mean return [a/n, b/n, c/n]; } return [a/wSum, b/wSum, c/wSum]; } // ---- True log-matrix as nonlinear function of x ---- // Uses a sigmoid-like transition to highlight boundary effects function trueLogMatrix(x) { var a = 0.5 + 1.5*x + 0.6*Math.sin(2.5*Math.PI*x); var b = -0.3 + 0.8*x + 0.4*Math.cos(2.0*Math.PI*x); var c = 0.8 + 1.0*x + 0.5*Math.sin(1.8*Math.PI*x + 0.5); return [a, b, c]; } // ---- Data generation ---- 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(); 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 }; } // ---- 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; }); trueLogs_app = xGrid_app.map(function(x) { return trueLogMatrix(x); }); trueTraces = trueLogs_app.map(symTrace); trueDets = trueLogs_app.map(symDet); // ---- Global Fréchet (Lecture 6) ---- 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); }); // ---- Kernel Fréchet (Lecture 7) ---- 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); }); // ---- Local-linear Fréchet (Lecture 8) ---- localLinLogs = 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); // Local moments var nu0 = 0, nu1 = 0, nu2 = 0; for (var i = 0; i < n_app; i++) { var di = appData.X[i] - x; nu0 += rawW[i]; nu1 += rawW[i] * di; nu2 += rawW[i] * di * di; } nu0 /= n_app; nu1 /= n_app; nu2 /= n_app; var sigma02 = nu0 * nu2 - nu1 * nu1; var weights; if (Math.abs(sigma02) > 1e-10) { weights = appData.X.map(function(xi, i) { return rawW[i] * (nu2 - nu1*(xi - x)) / sigma02; }); } else { weights = wSum > 1e-10 ? rawW.map(function(w) { return w/wSum; }) : rawW.map(function() { return 1/n_app; }); } return symWeightedMean(appData.logs, weights); }); // ---- Scalar summaries ---- globalTraces = globalLogs.map(symTrace); globalDets = globalLogs.map(symDet); kernelTraces = kernelLogs.map(symTrace); kernelDets = kernelLogs.map(symDet); llTraces = localLinLogs.map(symTrace); llDets = localLinLogs.map(symDet); dataTraces = appData.logs.map(symTrace); dataDets = appData.logs.map(symDet); // Unconditional Fréchet mean unifW = appData.logs.map(function() { return 1; }); meanLog = symWeightedMean(appData.logs, unifW); meanSPD = symExp(meanLog); // Local-linear weight diagnostics at left and right boundaries function computeLLWeightsAt(xEval) { var rawW = appData.X.map(function(xi) { return gaussKernel((xi - xEval)/h_app); }); var nu0 = 0, nu1 = 0, nu2 = 0; for (var i = 0; i < n_app; i++) { var di = appData.X[i] - xEval; nu0 += rawW[i]; nu1 += rawW[i]*di; nu2 += rawW[i]*di*di; } nu0 /= n_app; nu1 /= n_app; nu2 /= n_app; var sigma02 = nu0*nu2 - nu1*nu1; if (Math.abs(sigma02) > 1e-10) { return appData.X.map(function(xi, i) { return rawW[i] * (nu2 - nu1*(xi - xEval)) / sigma02; }); } var wSum = rawW.reduce(function(a,b) { return a+b; }, 0); return rawW.map(function(w) { return w/wSum; }); } llW_left = computeLLWeightsAt(0.05); llW_right = computeLLWeightsAt(0.95); negLeft = llW_left.filter(function(w) { return w < 0; }).length; negRight = llW_right.filter(function(w) { return w < 0; }).length; // ---- Visualization ---- function buildAppPlot(title, yLabel, trueVals, globalVals, kernelVals, llVals, dataVals, xVals, xGrid, yMin, yMax) { var W = 500, H = 210; 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"/>'); 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.4"/>'); }); // True (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 (red, thinner) 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="1.6" opacity="0.65"/>'); // Kernel (orange) 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="#f08c00" stroke-width="2.0"/>'); // Local-linear (blue, thicker) d = ""; for (var j = 0; j < xGrid.length; j++) d += (j===0?'M':'L') + " " + xS(xGrid[j]).toFixed(1) + " " + yS(llVals[j]).toFixed(1); parts.push('<path d="' + d + '" fill="none" stroke="#1971c2" stroke-width="2.6"/>'); // Boundary shading (left and right 10% of predictor range) parts.push('<rect x="' + margin.left + '" y="' + margin.top + '" width="' + (0.1*pw) + '" height="' + ph + '" fill="#795548" opacity="0.05"/>'); parts.push('<rect x="' + (margin.left + 0.9*pw) + '" y="' + margin.top + '" width="' + (0.1*pw) + '" height="' + ph + '" fill="#795548" opacity="0.05"/>'); // 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>'); // Boundary labels parts.push('<text x="' + (margin.left + 0.05*pw) + '" y="' + (margin.top+ph-4) + '" text-anchor="middle" font-size="7" fill="#795548" opacity="0.7">boundary</text>'); parts.push('<text x="' + (margin.left + 0.95*pw) + '" y="' + (margin.top+ph-4) + '" text-anchor="middle" font-size="7" fill="#795548" opacity="0.7">boundary</text>'); // Legend var lx = margin.left + pw - 170, ly = margin.top + 3; parts.push('<line x1="' + lx + '" y1="' + ly + '" x2="' + (lx+12) + '" y2="' + ly + '" stroke="#2e7d32" stroke-width="1.6" stroke-dasharray="4,3"/>'); parts.push('<text x="' + (lx+15) + '" y="' + (ly+4) + '" font-size="7.5" fill="#616161">True</text>'); parts.push('<line x1="' + (lx+42) + '" y1="' + ly + '" x2="' + (lx+54) + '" y2="' + ly + '" stroke="#e53935" stroke-width="1.6"/>'); parts.push('<text x="' + (lx+57) + '" y="' + (ly+4) + '" font-size="7.5" fill="#616161">Glob</text>'); parts.push('<line x1="' + (lx+80) + '" y1="' + ly + '" x2="' + (lx+92) + '" y2="' + ly + '" stroke="#f08c00" stroke-width="1.8"/>'); parts.push('<text x="' + (lx+95) + '" y="' + (ly+4) + '" font-size="7.5" fill="#616161">Kern</text>'); parts.push('<line x1="' + (lx+120) + '" y1="' + ly + '" x2="' + (lx+132) + '" y2="' + ly + '" stroke="#1971c2" stroke-width="2.2"/>'); parts.push('<text x="' + (lx+135) + '" y="' + (ly+4) + '" font-size="7.5" fill="#616161">LL</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:700px;">\n' + parts.join('\n') + '\n</svg>'; } // Compute y ranges allTraceVals = dataTraces.concat(trueTraces, globalTraces, kernelTraces, llTraces); traceMin = Math.min.apply(null, allTraceVals); traceMax = Math.max.apply(null, allTraceVals); tracePad = (traceMax - traceMin) * 0.12; traceMinAdj = traceMin - tracePad; traceMaxAdj = traceMax + tracePad; allDetVals = dataDets.concat(trueDets, globalDets, kernelDets, llDets); detMin = Math.min.apply(null, allDetVals); detMax = Math.max.apply(null, allDetVals); detPad = (detMax - detMin) * 0.12; detMinAdj = detMin - detPad; detMaxAdj = detMax + detPad; plot_spd_view = html` <div style="width:80%; max-width:700px; margin:0 auto; font-family:system-ui,sans-serif;"> ${buildAppPlot("Trace of log Σ(x) vs. Stress", "tr(log Σ)", trueTraces, globalTraces, kernelTraces, llTraces, dataTraces, appData.X, xGrid_app, traceMinAdj, traceMaxAdj)} <div style="margin-top:10px;"> ${buildAppPlot("Determinant of log Σ(x) vs. Stress", "det(log Σ)", trueDets, globalDets, kernelDets, llDets, dataDets, appData.X, xGrid_app, detMinAdj, detMaxAdj)} </div> <p style="margin-top:12px; font-size:0.85em; color:#495057;"> <b>Observe the boundary regions</b> (shaded at x ≈ 0 and x ≈ 1): Local-linear (blue) tracks the true curve at the edges where kernel (orange) is pulled inward. Global (red) misses the nonlinearity entirely. Bandwidth h = ${h_app.toFixed(2)}. </p> <div style="background:#f8f9fa; border-radius:6px; padding:10px 14px; margin-top:8px; font-size:0.85em;"> <b>Local-linear weight diagnostics at boundaries:</b> <table style="width:100%; border-collapse:collapse; margin-top:4px;"> <tr><td style="padding:2px 8px;">At left boundary (x = 0.05)</td> <td style="padding:2px 8px; text-align:right;">${negLeft} of ${n_app} weights negative (${(negLeft/n_app*100).toFixed(1)}%)</td></tr> <tr><td style="padding:2px 8px;">At right boundary (x = 0.95)</td> <td style="padding:2px 8px; text-align:right;">${negRight} of ${n_app} weights negative (${(negRight/n_app*100).toFixed(1)}%)</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_spd_view}${plot_spd_view}</div>` ``` ::: {.callout-tip title="Try these experiments"} - **Focus on the shaded boundary regions** ($x$ near 0 and $x$ near 1): The kernel fit (orange) consistently lies *inside* the true curve — this is the $O(h)$ boundary bias. The local-linear fit (blue) stays closer to the true curve. - **Decrease $h$ (0.05–0.10):** The boundary bias of kernel regression becomes more localized but more severe. The local-linear weights become more negative (check the diagnostics panel) as the slope correction becomes more aggressive. - **Increase $h$ (0.4–0.6):** All three methods converge toward the unconditional Fréchet mean. The boundary advantage of local-linear diminishes as the local neighborhood expands to include the whole dataset. - **Increase $n$ (150–200):** The local-linear advantage at boundaries persists even with large samples — boundary bias is a *bias* problem, not a variance problem, so it doesn't vanish with more data. - **Compare trace and determinant:** In some datasets, boundary bias may be more visible in one scalar summary than the other — this reflects which aspects of the covariance structure are most sensitive to the local slope at the edge. ::: ### Takeaway Local-linear Fréchet regression inherits the computational simplicity of kernel regression (weighted Fréchet means on the SPD manifold) while correcting its primary weakness: $O(h)$ boundary bias. In the market risk context, where the boundaries correspond to extreme market conditions, this correction is not a theoretical nicety — it directly affects the accuracy of stress-tested portfolio risk estimates. ## Key Takeaways - **Local Fréchet regression** extends local linear regression to metric spaces by constructing a predictor-side weight function $s(z, x, h)$ that automatically adjusts for the local slope. - The **weights can be negative** — this is the mechanism that corrects boundary bias. - The convergence rate $O_p(h^{2/(\beta_1-1)} + (nh)^{-1/[2(\beta_2-1)]})$ generalizes the classical bias–variance decomposition, with the exponents $\beta_1, \beta_2$ determined by the local curvature of the Fréchet objective. - For **quadratic curvature** ($\beta_1 = \beta_2 = 2$), the optimal bandwidth is $h \asymp n^{-1/5}$, yielding the familiar nonparametric rate $n^{-2/5}$. - Local-linear is almost always preferred over local-constant in practice — the computational cost is nearly identical, but the boundary behavior is substantially better. ## Exercises 1. **Local-linear weights sum to one.** Show that $\sum_{i=1}^n s_{in}(x, h) = 1$. Also show that in the Euclidean case, $\sum_i s_{in}(x, h) X_i = x$ when the model includes an intercept and slope. <a href="javascript:void(0)" onclick="showSolution('l8-sol-1')" class="solution-link">📝 Show Solution</a> 2. **Optimal bandwidth derivation.** Derive the optimal bandwidth $h \asymp n^{-1/5}$ for the case $\beta_1 = \beta_2 = 2$ by balancing the squared bias and variance terms. <a href="javascript:void(0)" onclick="showSolution('l8-sol-2')" class="solution-link">📝 Show Solution</a> 3. **Higher dimensions.** Discuss the challenges of extending local Fréchet regression to $p > 1$ predictors. What is the bias rate in terms of $h$ and $p$? <a href="javascript:void(0)" onclick="showSolution('l8-sol-3')" class="solution-link">📝 Show Solution</a> 4. **Computation.** Discuss how to implement local Fréchet regression in general metric spaces and manifolds. Note that the weights could be negative, and thus the weighted proximal point algorithm may not converge in general. <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="l8-sol-1"> <div class="solution-dialog"> <div class="solution-dialog-header"> <h4>Exercise 1: Weight Sum</h4> <button onclick="closeSolution('l8-sol-1')" class="solution-dialog-close">✕ Close</button> </div> <div class="solution-original"><strong>Exercise:</strong> Show that $\sum_i s_{in}(x, h) = 1$ and the first-moment property.</div> <div class="solution-answer"> <strong>Solution:</strong> For the sum: $$ \sum_i s_{in}(x, h) = \sum_i \frac{K_h(X_i - x)\{\hat{\nu}_2 - \hat{\nu}_1(X_i - x)\}}{\hat{\sigma}_0^2} = \frac{\hat{\nu}_2\sum_i K_h(X_i - x) - \hat{\nu}_1\sum_i K_h(X_i - x)(X_i - x)}{\hat{\sigma}_0^2}. $$ But $\sum_i K_h(X_i - x) = n\hat{\nu}_0$ and $\sum_i K_h(X_i - x)(X_i - x) = n\hat{\nu}_1$. So the numerator is $n(\hat{\nu}_2\hat{\nu}_0 - \hat{\nu}_1\hat{\nu}_1) = n\hat{\sigma}_0^2$, and the sum is 1. For the first-moment property: $$ \begin{aligned} \sum_i s_{in}(x, h) X_i &= \sum_i \frac{K_h(X_i - x)\{\hat{\nu}_2 - \hat{\nu}_1(X_i - x)\}}{\hat{\sigma}_0^2} X_i \\[4pt] &= \frac{1}{\hat{\sigma}_0^2}\left[ \hat{\nu}_2\sum_i K_h(X_i - x) X_i - \hat{\nu}_1\sum_i K_h(X_i - x)(X_i - x) X_i \right]. \end{aligned} $$ Write $X_i = (X_i - x) + x$ in both inner sums: $$ \begin{aligned} \sum_i K_h(X_i - x) X_i &= \sum_i K_h(X_i - x)(X_i - x) + x\sum_i K_h(X_i - x) = n\hat{\nu}_1 + x n\hat{\nu}_0, \\[4pt] \sum_i K_h(X_i - x)(X_i - x) X_i &= \sum_i K_h(X_i - x)(X_i - x)^2 + x\sum_i K_h(X_i - x)(X_i - x) = n\hat{\nu}_2 + x n\hat{\nu}_1. \end{aligned} $$ Substituting back: $$ \begin{aligned} \sum_i s_{in}(x, h) X_i &= \frac{1}{\hat{\sigma}_0^2}\left[ \hat{\nu}_2(n\hat{\nu}_1 + x n\hat{\nu}_0) - \hat{\nu}_1(n\hat{\nu}_2 + x n\hat{\nu}_1) \right] \\[4pt] &= \frac{n}{\hat{\sigma}_0^2}\left[ \hat{\nu}_1\hat{\nu}_2 + x\hat{\nu}_0\hat{\nu}_2 - \hat{\nu}_1\hat{\nu}_2 - x\hat{\nu}_1^2 \right] \\[4pt] &= \frac{n}{\hat{\sigma}_0^2}\left[ x(\hat{\nu}_0\hat{\nu}_2 - \hat{\nu}_1^2) \right] = \frac{x n \hat{\sigma}_0^2}{\hat{\sigma}_0^2} = x. \end{aligned} $$ This means the local-linear weights exactly reproduce a linear function — a property local-constant weights do not possess. </div> </div> </dialog> <dialog id="l8-sol-2"> <div class="solution-dialog"> <div class="solution-dialog-header"> <h4>Exercise 2: Optimal Bandwidth</h4> <button onclick="closeSolution('l8-sol-2')" class="solution-dialog-close">✕ Close</button> </div> <div class="solution-original"><strong>Exercise:</strong> Derive the optimal bandwidth $h \asymp n^{-1/5}$.</div> <div class="solution-answer"> <strong>Solution:</strong> For $\beta_1 = \beta_2 = 2$, the bias rate is $h^2$ and the variance rate is $(nh)^{-1/2}$. The MSE is: $$ \text{MSE} \asymp (h^2)^2 + ((nh)^{-1/2})^2 = h^4 + (nh)^{-1}. $$ To find the optimal $h$, differentiate with respect to $h$ and set to zero: $$ \frac{d}{dh}\bigl(h^4 + n^{-1}h^{-1}\bigr) = 4h^3 - n^{-1}h^{-2} = 0. $$ Solving: $4h^5 = n^{-1}$, so $h = (4n)^{-1/5} \asymp n^{-1/5}$. Substituting back: $\text{MSE} \asymp n^{-4/5} + n^{-1} \cdot n^{1/5} = n^{-4/5}$, giving the rate $n^{-2/5}$ for $d(\hat{\mu}(x), \mu(x))$. </div> </div> </dialog> <dialog id="l8-sol-3"> <div class="solution-dialog"> <div class="solution-dialog-header"> <h4>Exercise 3: Higher Predictor Dimensions</h4> <button onclick="closeSolution('l8-sol-3')" class="solution-dialog-close">✕ Close</button> </div> <div class="solution-original"><strong>Exercise:</strong> Discuss challenges of extending local Fréchet regression to $p > 1$ predictors.</div> <div class="solution-answer"> <strong>Solution:</strong> For $p > 1$, local linear regression fits a hyperplane: $\beta_0 + \beta_1^\top(X - x)$. The weight function generalizes to $$ s(z, x, h) = \frac{1}{n} K_h(z - x) \bigl[1 - \hat{\nu}_1(x, h)^\top \hat{\Sigma}(x, h)^{-1}(z - x)\bigr], $$ where $\hat{\nu}_1(x, h)$ is a $p$-vector and $\hat{\Sigma}(x, h)$ is a $p \times p$ local covariance matrix. **Challenges:** 1. **Curse of dimensionality:** The bias rate is $O(h^2)$ (unchanged), but the variance rate becomes $(nh^p)^{-1/2}$. The optimal bandwidth is $h \asymp n^{-1/(p+4)}$, and the optimal MSE rate degrades to $n^{-4/(p+4)}$ — much slower for large $p$. 2. **Effective sample size:** For reasonable bandwidths, the number of observations in a $p$-dimensional neighborhood shrinks exponentially with $p$. **Practical implications:** Local Fréchet regression may work well for $p \le 5$. For larger $p$, dimension reduction (e.g., sufficient dimension reduction, additive models) or global Fréchet regression with regularization is recommended. The @PetersenMueller2019 paper handles general $p$; the weight formula involves the full $p$-dimensional local linear expansion. </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 8 Quiz — 10 Multiple-Choice Questions](../quizzes/lecture-08-quiz.qmd)** ## Further Reading - @PetersenMueller2019 — Sections 4 and S.3: the full theory of local Fréchet regression.