---
title: "Lecture 6: Fréchet Regression — Global Fréchet Regression"
subtitle: "Extending linear regression to metric-space responses"
format:
html:
code-fold: true
code-tools: true
code-copy: true
pdf:
documentclass: scrartcl
pdf-engine: xelatex
toc: true
number-sections: true
geometry:
- margin=1in
colorlinks: true
bibliography: source/ref.bib
---
## Learning Goals
By the end of this lecture, learners should be able to:
- Define the conditional Fréchet mean and Fréchet regression function in a general metric space.
- Explain how ordinary linear regression can be rewritten as a weighted least-squares problem to facilitate metric-space generalization.
- State the definition of the global Fréchet regression function using the weight function $s(z,x)$.
- Interpret the three key properties of global Fréchet regression (recovery of the Fréchet mean at $\mathbb{E}X$, nonlocal weights, and recovery of Euclidean linear regression).
- Formulate the sample global Fréchet regression estimator as a weighted Fréchet mean problem.
- Describe the consistency and convergence-rate theory, including the roles of the entropy condition (P1) and curvature condition (P2).
## Motivation: Regression with Non-Euclidean Responses
The basic regression problem is to relate a Euclidean predictor to a response that lives in a nonlinear metric space. The framework of @PetersenMueller2019 extends conditional mean regression to this setting by replacing ordinary squared-error loss with squared metric distance.
Throughout this lecture, let $X \in \mathbb{R}^p$ be a predictor and let $Y$ be a random element of a metric space $(\mathcal{M}, d)$. We assume that the conditional laws of $Y$ given $X = x$ exist.
::: {#def-conditional-frechet-mean .definition title="Conditional Fréchet mean and Fréchet regression function"}
For each $x \in \mathbb{R}^p$, define the **conditional Fréchet function**
$$
F(\omega, x) = \mathbb{E}\{d^2(Y, \omega) \mid X = x\}, \qquad \omega \in \mathcal{M}.
$$
Any minimizer of $F(\cdot, x)$ is called a **conditional Fréchet mean** of $Y$ given $X = x$. If the minimizer is unique, we denote it by
$$
\mu(x) = \arg\min_{\omega \in \mathcal{M}} F(\omega, x)
$$
and call $\mu$ the **Fréchet regression function** @PetersenMueller2019.
:::
This is the direct analogue of the conditional mean. When $\mathcal{M} = \mathbb{R}^p$ with the Euclidean metric, $\mu(x)$ is just the usual conditional expectation $\mathbb{E}(Y \mid X = x)$.
## From Linear Regression to Weighted Least Squares
Petersen and Müller propose a global model that mimics ordinary multiple linear regression while remaining intrinsic to the metric. To see their idea clearly, suppose first that $Y \in \mathbb{R}$. Write the usual linear regression function in centered form as
$$
\mu(x) = \bar{\beta}_0 + \bar{\beta}_1^\top(x - \mathbb{E}X), \qquad
\bar{\beta}_0 = \mathbb{E}(Y), \qquad
\bar{\beta}_1 = \Sigma^{-1}\mathbb{E}\{(X - \mathbb{E}X)Y\},
$$
where $\bar{\beta}_0$ denotes the true intercept, $\bar{\beta}_1$ denotes the slope vector, and $\Sigma$ denotes the covariance matrix of $X$.
Observe that
$$
(\bar{\beta}_0, \bar{\beta}_1) = \arg\min_{\beta_0 \in \mathbb{R}, \beta_1 \in \mathbb{R}^p} \mathbb{E}\left[\mathbb{E}(Y|X) - (\beta_0 + \beta_1^\top(X - \mathbb{E}X))\right]^2.
$$
::: {.callout-note collapse="true" title="Derivation: Why the minimizer recovers the OLS coefficients"}
Define the objective function
$$
g(\beta_0, \beta_1) = \mathbb{E}\Bigl[\mathbb{E}(Y|X) - \beta_0 - \beta_1^\top(X - \mathbb{E}X)\Bigr]^2.
$$
This is a convex quadratic in $(\beta_0, \beta_1)$, so the unique global minimizer is found by setting the gradient to zero.
**Step 1 — Minimize with respect to $\beta_0$:**
$$
\frac{\partial g}{\partial \beta_0} = -2\,\mathbb{E}\Bigl[\mathbb{E}(Y|X) - \beta_0 - \beta_1^\top(X - \mathbb{E}X)\Bigr] = 0.
$$
Distribute the expectation and use $\mathbb{E}[X - \mathbb{E}X] = 0$ together with $\mathbb{E}[\mathbb{E}(Y|X)] = \mathbb{E}Y$ (law of iterated expectations):
$$
\mathbb{E}Y - \beta_0 = 0 \quad\Longrightarrow\quad \bar{\beta}_0 = \mathbb{E}Y.
$$
**Step 2 — Minimize with respect to $\beta_1$:**
$$
\nabla_{\beta_1} g = -2\,\mathbb{E}\Bigl[(X - \mathbb{E}X)\bigl(\mathbb{E}(Y|X) - \beta_0 - \beta_1^\top(X - \mathbb{E}X)\bigr)\Bigr] = 0.
$$
Substitute $\beta_0 = \mathbb{E}Y$ and split:
$$
\mathbb{E}\Bigl[(X - \mathbb{E}X)\bigl(\mathbb{E}(Y|X) - \mathbb{E}Y\bigr)\Bigr]
- \mathbb{E}\Bigl[(X - \mathbb{E}X)\,\beta_1^\top(X - \mathbb{E}X)\Bigr] = 0.
$$
**First term:** Since $(X - \mathbb{E}X)$ is $X$-measurable, pull it inside the conditional expectation:
$$
\begin{aligned}
\mathbb{E}\Bigl[(X - \mathbb{E}X)\bigl(\mathbb{E}(Y|X) - \mathbb{E}Y\bigr)\Bigr]
&= \mathbb{E}\Bigl[\mathbb{E}\bigl((X - \mathbb{E}X)(Y - \mathbb{E}Y) \mid X\bigr)\Bigr] \\
&= \mathbb{E}\bigl[(X - \mathbb{E}X)(Y - \mathbb{E}Y)\bigr] = \operatorname{Cov}(X, Y) \equiv \sigma_{YX}.
\end{aligned}
$$
**Second term:** Since $\beta_1^\top(X - \mathbb{E}X)$ is a scalar, rewrite as an outer product:
$$
(X - \mathbb{E}X)\,\beta_1^\top(X - \mathbb{E}X)
= (X - \mathbb{E}X)(X - \mathbb{E}X)^\top \beta_1,
$$
so taking expectations gives $\mathbb{E}[(X - \mathbb{E}X)(X - \mathbb{E}X)^\top]\,\beta_1 = \operatorname{Var}(X)\,\beta_1 \equiv \Sigma\beta_1$.
**Step 3 — Solve the normal equation:**
$$
\sigma_{YX} - \Sigma\beta_1 = 0 \quad\Longrightarrow\quad \bar{\beta}_1 = \Sigma^{-1}\sigma_{YX}.
$$
Thus the minimizer of $g(\beta_0, \beta_1)$ is exactly $(\mathbb{E}Y,\; \Sigma^{-1}\sigma_{YX})$, the standard OLS coefficients in centered form.
:::
Denoting $\sigma_{YX} = \mathbb{E}[Y(X - \mathbb{E}X)]$, we have the normal equations:
$$
\mathbb{E}Y - \beta_0 = 0, \qquad \sigma_{YX} - \Sigma\beta_1 = 0,
$$
and the solutions are $\bar{\beta}_1 = \Sigma^{-1}\sigma_{YX}$ and $\bar{\beta}_0 = \mathbb{E}Y$. Therefore,
$$
\begin{aligned}
\mu(x) &= \bar{\beta}_0 + \bar{\beta}_1^\top(x - \mathbb{E}X) \\
&= \mathbb{E}Y + \sigma_{YX}^\top \Sigma^{-1}(x - \mathbb{E}X) \\
&= \mathbb{E}\{Y[1 + (X - \mathbb{E}X)^\top\Sigma^{-1}(x - \mathbb{E}X)]\} \\
&= \mathbb{E}[Y s(X, x)],
\end{aligned}
$$
where $s(z, x) = 1 + (z - \mathbb{E}X)^\top \Sigma^{-1}(x - \mathbb{E}X)$. Since $\mathbb{E}\, s(X, x) = 1$, for fixed $x$, the minimizer of $\mathbb{E}[s(X, x)(Y - y)^2]$ over $y$ is $\mathbb{E}[Y s(X, x)]$, i.e.,
$$
\mu(x) = \arg\min_{y \in \mathbb{R}} \mathbb{E}[s(X, x) d_E^2(Y, y)],
$$
where $d_E$ denotes the Euclidean distance.
::: {.callout-tip title="The key insight"}
This reformulation is the key observation of Section 2 in @PetersenMueller2019. Once linear regression is written this way, all linear structure is absorbed into the predictor-side weight $s(X, x)$, while the response enters only through squared distance. The extension to a general metric space is then immediate: replace $d_E(Y, y)$ by $d(Y, \omega)$ and minimize over $\omega \in \mathcal{M}$.
:::
## Global Fréchet Regression
::: {#def-global-frechet-regression .definition title="Global Fréchet regression"}
Assume that $\mathbb{E}X$ and $\Sigma = \operatorname{var}(X)$ exist and that $\Sigma$ is positive definite. Define the weight function
$$
s(z, x) = 1 + (z - \mathbb{E}X)^\top\Sigma^{-1}(x - \mathbb{E}X), \qquad z, x \in \mathbb{R}^p.
$$
The **global Fréchet regression function** is then defined by
$$
\mu(x) = \arg\min_{\omega \in \mathcal{M}} M(\omega, x), \qquad
M(\omega, x) = \mathbb{E}\{s(X, x) d^2(Y, \omega)\}.
$$
:::
Several features are immediate:
1. **Recovery of the Fréchet mean at $\mathbb{E}X$.** Because $s(\cdot, \mathbb{E}X) \equiv 1$, one has
$$
\mu(\mathbb{E}X) = \arg\min_{\omega \in \mathcal{M}} \mathbb{E}\{d^2(Y, \omega)\},
$$
and thus the global regression curve passes through the unconditional Fréchet mean.
2. **Nonlocal weights.** The weights $s(z, x)$ need not be nonnegative and do not vanish far from $x$; this is the metric analogue of a **global** linear model rather than a local smoother.
3. **Recovery of Euclidean linear regression.** If $\mathcal{M} = \mathbb{R}^q$ with the Euclidean metric, the minimizer is exactly the usual multivariate linear regression fit.
## Sample Estimator
::: {#def-sample-global-frechet .definition title="Sample global Fréchet regression estimator"}
Given independent observations $(X_i, Y_i)$, $i = 1, \ldots, n$, define
$$
\bar{X} = \frac{1}{n}\sum_{i=1}^n X_i, \qquad
\hat{\Sigma} = \frac{1}{n}\sum_{i=1}^n (X_i - \bar{X})(X_i - \bar{X})^\top,
$$
and empirical weights
$$
s_{in}(x) = 1 + (X_i - \bar{X})^\top \hat{\Sigma}^{-1}(x - \bar{X}).
$$
The empirical objective is
$$
M_n(\omega, x) = \frac{1}{n}\sum_{i=1}^n s_{in}(x) d^2(Y_i, \omega),
$$
and the **global sample Fréchet regression estimator** is
$$
\hat{\mu}(x) = \arg\min_{\omega \in \mathcal{M}} M_n(\omega, x).
$$
:::
Thus estimation reduces to a **weighted Fréchet mean** problem for each predictor value $x$. Computationally, once one can minimize weighted Fréchet functions in the response space, one can fit the regression estimator.
## Consistency and Convergence Rates
Section 3 of @PetersenMueller2019 studies consistency and convergence rates of the estimator on a totally bounded metric space $(\mathcal{M}, d)$. For a fixed $x \in \mathbb{R}^p$, the imposed assumptions are:
::: {.callout-note title="Assumptions (P0)–(P2)"}
**(P0) Existence and separation.** The objects $\mu(x)$ and $\hat{\mu}(x)$ exist and are unique, the latter almost surely, and for any $\varepsilon > 0$,
$$
\inf_{d(\omega, \mu(x)) > \varepsilon} M(\omega, x) > M(\mu(x), x).
$$
This is the standard separation assumption used in the consistency theory of M-estimators [@vanDerVaartWellner1996, Chapter 3.2]. Conceptually, it is the hypothesis behind the usual argmin theorem: once the empirical criterion $M_n(\cdot, x)$ converges suitably to the population criterion $M(\cdot, x)$, the corresponding minimizer $\hat{\mu}(x)$ is forced to lie close to $\mu(x)$.
**(P1) Local entropy near the target.** 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.
$$
**(P2) Curvature (margin) condition.** There exist $\eta > 0$, $C > 0$, and $\beta > 1$, possibly depending on $x$, such that whenever $d(\mu(x), \omega) < \eta$,
$$
M(\omega, x) - M(\mu(x), x) \ge C\, d(\omega, \mu(x))^{\beta}.
$$
:::
The entropy condition (P1) controls the local size of the parameter space near $\mu(x)$, and therefore the stochastic oscillation of $M_n - M$ on shrinking neighborhoods of the minimum. The curvature condition (P2) quantifies how sharply the population criterion rises away from its minimum. Together, these two ingredients are exactly what empirical process arguments need.
::: {#thm-global-consistency .theorem title="Consistency and convergence rate (Theorem 1 of Petersen & Müller 2019)"}
If (P0) holds and $\mathcal{M}$ is bounded, then for any fixed $x \in \mathbb{R}^p$,
$$
d(\hat{\mu}(x), \mu(x)) = o_p(1).
$$
If (P0)–(P2) hold, then for any fixed $x \in \mathbb{R}^p$,
$$
d(\hat{\mu}(x), \mu(x)) = O_p(n^{-1/(2\beta - 2)}).
$$
:::
::: {.callout-important title="Examples"}
For the examples introduced in the introductory lecture:
- For the **sphere**, (P0) and (P2) may not hold.
- For **SPD spaces** endowed with either the affine-invariant, log-Euclidean, or log-Cholesky distance, and for the **BHV space** and the **Wasserstein space $\mathcal{P}(\mathbb{R})$**, (P0)–(P2) hold with $\beta = 2$, yielding the pointwise rate $O_p(n^{-1/2})$.
:::
## Interactive Exploration: Global Fréchet Regression on the Sphere
The following demo implements global Fréchet regression where the response space is the positive quadrant (octant) of the unit sphere $\mathbb{S}^2$. The predictor $X$ is one-dimensional, and the response $Y$ lies on the sphere octant $\{x, y, z \ge 0\}$. We visualize the true regression curve, the data points, and the fitted global Fréchet regression curve together on the sphere.
**Data generation.** $n$ pairs $(X_i, Y_i)$ are generated with $X_i \sim \text{Uniform}(-1, 1)$. The true regression function $\mu(x)$ is a curve on the octant parameterized by spherical angles:
$$\mu(x) = (\sin\theta(x)\cos\phi(x),\; \sin\theta(x)\sin\phi(x),\; \cos\theta(x)),$$
where $\theta(x)$ and $\phi(x)$ are linear (or constant) functions of $x$, chosen by the regression-type selector. The response $Y_i$ is generated by moving a random fraction along the great-circle geodesic from a random octant point $U_i$ toward $\mu(X_i)$:
$$Y_i = \gamma_{U_i \to \mu(X_i)}(1 - \sigma\varepsilon_i), \qquad U_i \sim \text{Uniform(octant)},\; \varepsilon_i \sim \text{Uniform}(0,1).$$
**Global Fréchet regression estimator.** For any predictor value $x$, the sample estimator is the weighted Fréchet mean
$$\hat{\mu}(x) = \arg\min_{\omega \in \mathbb{S}^2} \frac{1}{n}\sum_{i=1}^n s_{in}(x)\, d^2(Y_i, \omega), \qquad s_{in}(x) = 1 + \frac{(X_i - \bar{X})(x - \bar{X})}{\hat{\sigma}_X^2}.$$
We compute each $\hat{\mu}(x)$ using **Riemannian gradient descent** (RGD) on $\mathbb{S}^2$, initialized at the normalized weighted Euclidean mean. The RGD iteratively updates the estimate by moving along geodesics in the direction of the negative Riemannian gradient $\nabla_\omega M_n = -\frac{2}{n}\sum_i s_{in}(x)\,\mathrm{Log}_\omega(Y_i)$.
**Visual guide:**
- **Red dots** = observed responses $Y_i$ (color intensity shows predictor value $X_i$)
- **Green curve** = true regression function $\mu(x)$
- **Red curve** = fitted global Fréchet regression $\hat{\mu}(x)$
- **Dashed line** = connector from fitted to true at the evaluation point
- **Gold diamond** = Fréchet mean of all $Y_i$ (where the regression curve passes through at $\bar{X}$)
```{ojs}
//| label: fig-global-frechet-sphere
//| fig-cap: "Interactive: Global Fréchet regression on the sphere octant"
// ---- Controls ----
n_pts_control = Inputs.range([20, 250], {step: 10, value: 80, label: "Sample size n"})
reg_type_control = Inputs.select(["Linear θ & φ", "Varying θ only", "Varying φ only", "Diagonal (θ=φ)"], {value: "Linear θ & φ", 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)
reg_type = Generators.input(reg_type_control)
noise_level = Generators.input(noise_level_control)
data_seed = Generators.input(data_seed_control)
eval_x = Generators.input(eval_x_control)
controls_view = html`
<style>
.gfr-control-panel {
max-width: 620px;
margin: 0 auto 12px;
font-family: system-ui, sans-serif;
font-size: 0.85em;
}
.gfr-control-panel > *,
.gfr-slider-grid > * {
min-width: 0;
margin: 0;
}
.gfr-slider-grid {
display: grid;
grid-template-columns: repeat(2, minmax(0, 1fr));
gap: 6px 20px;
margin-top: 6px;
}
@media (max-width: 500px) {
.gfr-slider-grid {
grid-template-columns: 1fr;
}
}
</style>
<div class="gfr-control-panel">
<div>${reg_type_control}</div>
<div class="gfr-slider-grid">
<div>${n_pts_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];
var ry = (1-t)*p[1] + t*q[1];
var 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;
var 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)));
}
// Riemannian Log map on S^2: Log_p(q) in T_p S^2
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])];
}
// Riemannian Exp map on S^2: Exp_p(v) for v in T_p S^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;
// Init: normalized weighted Euclidean mean (use absolute weights for safety)
var w = [0, 0, 0];
var 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];
// RGD iterations
var lr = initLr || 0.3;
for (var iter = 0; iter < maxIter; iter++) {
// Compute Riemannian gradient = -2/n * sum_i w_i * Log_w(Y_i)
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;
// Gradient descent step: w_new = Exp_w(-lr * grad)
var step = -lr;
var v = [step * grad[0], step * grad[1], step * grad[2]];
w = sphereExp(w, v);
// Decay learning rate
lr = initLr / (1 + iter / 80);
}
return w;
}
// ---- Orthographic projection (matching lecture 5) ----
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) {
// Rejection-free: sample z ~ Uniform(0,1), phi ~ Uniform(0, π/2)
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: maps x in [-1,1] to a point on the octant
function trueRegFunc(x, regType) {
var xNorm = (x + 1) / 2; // map to [0, 1]
var theta, phi;
if (regType === "Linear θ & φ") {
theta = 0.2 + 0.8 * xNorm * (Math.PI/2 - 0.1);
phi = 0.3 + 0.7 * xNorm * (Math.PI/2 - 0.1);
} else if (regType === "Varying θ only") {
theta = 0.15 + 0.9 * xNorm * (Math.PI/2 - 0.1);
phi = Math.PI/4;
} else if (regType === "Varying φ only") {
theta = Math.PI/4;
phi = 0.15 + 0.9 * xNorm * (Math.PI/2 - 0.1);
} else { // Diagonal
var ang = 0.2 + 0.8 * xNorm * (Math.PI/2 - 0.1);
theta = ang;
phi = ang;
}
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 = [];
var Y = [];
for (var i = 0; i < n; i++) {
var xi = (rng() - 0.5) * 2; // Uniform(-1, 1)
X.push(xi);
var mu_i = trueRegFunc(xi, regType);
var randPt = randomOctantPoint(rng);
var eps = rng();
var w = 1 - noise * eps;
var yi = sphereGeodesic(randPt, mu_i, w);
Y.push(yi);
}
return { X: X, Y: Y };
}
// ---- Global Fréchet regression computation ----
function computeGlobalFrechetCurve(X, Y, xGrid, rngSeed) {
var n = X.length;
// Sample mean and variance of X
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 = [];
var prevFit = null;
for (var j = 0; j < xGrid.length; j++) {
var x = xGrid[j];
// Compute weights s_in(x)
var weights = [];
for (var i = 0; i < n; i++) {
weights.push(1 + (X[i] - meanX) * (x - meanX) / varX);
}
// Use warm start: initialize from previous grid point's fit
var dataPts = Y;
// If we have a previous fit, use it as init (do fewer iters)
var maxIter = 150;
var initLr = 0.25;
// Warm start: override the init inside weightedSphereFrechetMean by using a custom init
var fit;
if (prevFit) {
// Warm start from previous fit with fewer iterations
fit = weightedSphereFrechetMeanWarm(dataPts, weights, prevFit, 80, 0.2);
} else {
fit = weightedSphereFrechetMean(dataPts, weights, maxIter, initLr);
}
fitted.push(fit);
prevFit = fit;
}
return fitted;
}
// Weighted Fréchet mean with explicit warm start
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;
var v = [step * grad[0], step * grad[1], step * grad[2]];
w = sphereExp(w, v);
lr = initLr / (1 + iter / 60);
}
return w;
}
// ---- Run everything ----
regData = generateSphereRegressionData(n_pts, reg_type, noise_level, data_seed);
// Evaluation grid
nGrid = 50;
xGridEval = Array.from({length: nGrid + 1}, function(_, j) { return -1 + j * 2 / nGrid; });
fittedCurve = computeGlobalFrechetCurve(regData.X, regData.Y, xGridEval, data_seed + 1000);
// True curve on the same grid
trueCurve = xGridEval.map(function(x) { return trueRegFunc(x, reg_type); });
// Fit at evaluation point
evalWeights = regData.X.map(function(xi) {
var meanX = regData.X.reduce(function(a,b) { return a+b; }, 0) / regData.X.length;
var varX = regData.X.reduce(function(s, xi) { var d = xi - meanX; return s + d*d; }, 0) / regData.X.length;
return 1 + (xi - meanX) * (eval_x - meanX) / varX;
});
fitAtEval = weightedSphereFrechetMean(regData.Y, evalWeights, 200, 0.3);
trueAtEval = trueRegFunc(eval_x, reg_type);
errorAtEval = sphereDist(fitAtEval, trueAtEval);
// Overall Fréchet mean (unweighted) — where regression passes through at meanX
unifWeights = regData.Y.map(function() { return 1; });
overallMean = weightedSphereFrechetMean(regData.Y, unifWeights, 200, 0.3);
// Sample statistics
meanX_val = regData.X.reduce(function(a,b) { return a+b; }, 0) / regData.X.length;
weightSum = evalWeights.reduce(function(a,b) { return a+b; }, 0);
negCount = evalWeights.filter(function(w) { return w < 0; }).length;
// ---- 3D Visualization of the sphere octant ----
function build3DView(regData, trueCurve, fittedCurve, xGridEval, fitAtEval, 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;
}
// Three boundary arcs of the octant
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);
// Full boundary polygon
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]);
// Grid arcs for constant z
function constantZArc(zVal, nPts) {
var r = Math.sqrt(1 - zVal*zVal);
var pts = [];
for (var i = 0; i <= nPts; i++) {
var alpha = (i / nPts) * Math.PI / 2;
pts.push([r * Math.cos(alpha), r * Math.sin(alpha), zVal]);
}
return pts;
}
var gridLines = [];
[0.2, 0.4, 0.6, 0.8].forEach(function(zv) {
gridLines.push(constantZArc(zv, 30));
});
var parts = [];
// Octant surface
parts.push('<path d="M ' + projectPolyline(boundaryPts) + ' Z" fill="#e3f2fd" stroke="none" opacity="0.45"/>');
// Grid lines
gridLines.forEach(function(arc) {
parts.push('<polyline points="' + projectPolyline(arc) + '" fill="none" stroke="#bbdefb" stroke-width="0.5"/>');
});
// Boundary arcs
[arc_xy, arc_yz, arc_zx].forEach(function(arc) {
parts.push('<polyline points="' + projectPolyline(arc) + '" fill="none" stroke="#64b5f6" stroke-width="2.0" stroke-linecap="round"/>');
});
// True regression curve (green, dashed, thicker)
var trueProj = trueCurve.map(function(p) { return project(p); });
var truePtsStr = trueProj.map(function(p) { return p[0].toFixed(1) + "," + p[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"/>');
// Fitted regression curve (red, solid, thicker)
var fittedProj = fittedCurve.map(function(p) { return project(p); });
var fittedPtsStr = fittedProj.map(function(p) { return p[0].toFixed(1) + "," + p[1].toFixed(1); }).join(" ");
parts.push('<polyline points="' + fittedPtsStr + '" fill="none" stroke="#e53935" stroke-width="3.0" stroke-linecap="round" stroke-linejoin="round" opacity="0.85"/>');
// Data points (colored by X value: low X = cool blue, high X = warm red)
var xMin = -1, xMax = 1;
regData.Y.forEach(function(yi, idx) {
var xi = regData.X[idx];
var t = (xi - xMin) / (xMax - xMin); // 0 (blue) to 1 (red)
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: fitted (red ring) and true (green ring)
var evalFitProj = project(fitAtEval);
var evalTrueProj = project(trueAtEval);
parts.push('<circle cx="' + evalFitProj[0].toFixed(1) + '" cy="' + evalFitProj[1].toFixed(1) +
'" r="8" fill="none" stroke="#e53935" stroke-width="3.5"/>');
parts.push('<circle cx="' + evalTrueProj[0].toFixed(1) + '" cy="' + evalTrueProj[1].toFixed(1) +
'" r="8" fill="none" stroke="#2e7d32" stroke-width="3.5" stroke-dasharray="5,3"/>');
// Dashed connector
parts.push('<line x1="' + evalFitProj[0].toFixed(1) + '" y1="' + evalFitProj[1].toFixed(1) +
'" x2="' + evalTrueProj[0].toFixed(1) + '" y2="' + evalTrueProj[1].toFixed(1) +
'" stroke="#795548" stroke-width="1.5" stroke-dasharray="4,3" opacity="0.7"/>');
// 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]);
var lvy = project([0, 1.15, 0]);
var 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 = 45;
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="10" fill="#37474f">Data point Yᵢ (color = Xᵢ)</text>');
parts.push('<line x1="' + (lx-6) + '" y1="' + (ly0+20) + '" x2="' + (lx+14) + '" y2="' + (ly0+20) +
'" stroke="#2e7d32" stroke-width="2.5" stroke-dasharray="5,3"/>');
parts.push('<text x="' + (lx+22) + '" y="' + (ly0+24) + '" font-size="10" fill="#37474f">True μ(x)</text>');
parts.push('<line x1="' + (lx-6) + '" y1="' + (ly0+38) + '" x2="' + (lx+14) + '" y2="' + (ly0+38) +
'" stroke="#e53935" stroke-width="2.5"/>');
parts.push('<text x="' + (lx+22) + '" y="' + (ly0+42) + '" font-size="10" fill="#37474f">Fitted μ̂(x)</text>');
parts.push('<polygon points="' + lx + ',' + (ly0+53) + ' ' + (lx+4) + ',' + (ly0+58) +
' ' + lx + ',' + (ly0+63) + ' ' + (lx-4) + ',' + (ly0+58) + '" fill="#f9a825"/>');
parts.push('<text x="' + (lx+12) + '" y="' + (ly0+62) + '" font-size="10" fill="#37474f">Fréchet mean μ̂(X̄)</text>');
parts.push('<circle cx="' + lx + '" cy="' + (ly0+78) + '" r="5" fill="none" stroke="#e53935" stroke-width="2.5"/>');
parts.push('<text x="' + (lx+12) + '" y="' + (ly0+82) + '" font-size="10" fill="#37474f">Fitted at eval x</text>');
parts.push('<circle cx="' + lx + '" cy="' + (ly0+96) + '" r="5" fill="none" stroke="#2e7d32" stroke-width="2.5" stroke-dasharray="4,2"/>');
parts.push('<text x="' + (lx+12) + '" y="' + (ly0+100) + '" font-size="10" fill="#37474f">True at eval x</text>');
return '<svg viewBox="0 0 520 420" xmlns="http://www.w3.org/2000/svg" style="width:100%;max-width:540px;">\n' +
parts.join('\n') + '\n</svg>';
}
sphere_view = html`<div style="text-align:center;">${build3DView(regData, trueCurve, fittedCurve, xGridEval, fitAtEval, 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 };
}
// Compute angles for all points
dataAngles = regData.Y.map(function(yi) { return sphericalAngles(yi); });
trueAngles = trueCurve.map(function(p) { return sphericalAngles(p); });
fittedAngles = fittedCurve.map(function(p) { return sphericalAngles(p); });
evalTrueAng = sphericalAngles(trueAtEval);
evalFitAng = sphericalAngles(fitAtEval);
function buildUnwrappedPlot(title, yLabel, trueVals, fittedVals, dataVals, xVals, xGrid, yMin, yMax) {
var W = 520, H = 190;
var margin = {top: 15, 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 = [];
// Background
parts.push('<rect x="' + margin.left + '" y="' + margin.top + '" width="' + pw +
'" height="' + ph + '" fill="#fafafa" rx="2"/>');
// Grid lines
for (var gy = 0; gy <= 4; gy++) {
var yv = yMin + (gy/4) * (yMax - yMin);
var yp = yToSvg(yv);
parts.push('<line x1="' + margin.left + '" y1="' + yp + '" x2="' + (margin.left+pw) +
'" y2="' + yp + '" stroke="#e0e0e0" stroke-width="0.5"/>');
parts.push('<text x="' + (margin.left-4) + '" y="' + (yp+4) +
'" text-anchor="end" font-size="8" fill="#9e9e9e">' + yv.toFixed(2) + '</text>');
}
// Data points
dataVals.forEach(function(dv, idx) {
var sx = xToSvg(xVals[idx]);
var sy = yToSvg(dv);
parts.push('<circle cx="' + sx + '" cy="' + sy + '" r="2.2" fill="#90caf9" opacity="0.55"/>');
});
// True curve
var trueLine = "";
for (var j = 0; j < xGrid.length; j++) {
trueLine += (j === 0 ? 'M' : 'L') + " " + xToSvg(xGrid[j]).toFixed(1) + " " + yToSvg(trueVals[j]).toFixed(1);
}
parts.push('<path d="' + trueLine + '" fill="none" stroke="#2e7d32" stroke-width="2.2" stroke-dasharray="7,4"/>');
// Fitted curve
var fitLine = "";
for (var j = 0; j < xGrid.length; j++) {
fitLine += (j === 0 ? 'M' : 'L') + " " + xToSvg(xGrid[j]).toFixed(1) + " " + yToSvg(fittedVals[j]).toFixed(1);
}
parts.push('<path d="' + fitLine + '" fill="none" stroke="#e53935" stroke-width="2.2"/>');
// Evaluation point markers
var exSvg = xToSvg(eval_x);
parts.push('<line x1="' + exSvg + '" y1="' + margin.top + '" x2="' + exSvg +
'" y2="' + (margin.top+ph) + '" stroke="#795548" stroke-width="1.2" stroke-dasharray="3,3" opacity="0.5"/>');
// Axes
parts.push('<line x1="' + margin.left + '" y1="' + (margin.top+ph) + '" x2="' +
(margin.left+pw) + '" y2="' + (margin.top+ph) + '" stroke="#424242" stroke-width="1"/>');
parts.push('<line x1="' + margin.left + '" y1="' + margin.top + '" x2="' +
margin.left + '" y2="' + (margin.top+ph) + '" stroke="#424242" stroke-width="1"/>');
// Axis labels
parts.push('<text x="' + (margin.left+pw/2) + '" y="' + (H-4) +
'" text-anchor="middle" font-size="10" fill="#616161">Predictor X</text>');
parts.push('<text x="' + (margin.left-34) + '" y="' + (margin.top+ph/2) +
'" text-anchor="middle" font-size="10" fill="#616161" ' +
'transform="rotate(-90,' + (margin.left-34) + ',' + (margin.top+ph/2) + ')">' + yLabel + '</text>');
// Title
parts.push('<text x="' + (margin.left+pw/2) + '" y="' + (margin.top-2) +
'" text-anchor="middle" font-size="11" font-weight="bold" fill="#37474f">' + title + '</text>');
// Legend
var lx = margin.left + pw - 140, ly = margin.top + 4;
parts.push('<line x1="' + lx + '" y1="' + ly + '" x2="' + (lx+18) + '" y2="' + ly +
'" stroke="#2e7d32" stroke-width="2" stroke-dasharray="5,3"/>');
parts.push('<text x="' + (lx+22) + '" y="' + (ly+4) + '" font-size="9" fill="#616161">True</text>');
parts.push('<line x1="' + (lx+60) + '" y1="' + ly + '" x2="' + (lx+78) + '" y2="' + ly +
'" stroke="#e53935" stroke-width="2"/>');
parts.push('<text x="' + (lx+82) + '" y="' + (ly+4) + '" font-size="9" fill="#616161">Fitted</text>');
parts.push('<rect x="' + margin.left + '" y="' + margin.top + '" width="' + pw +
'" height="' + ph + '" fill="none" stroke="#dee2e6"/>');
return '<svg viewBox="0 0 ' + W + ' ' + H + '" xmlns="http://www.w3.org/2000/svg" style="width:100%;max-width:540px;">\n' +
parts.join('\n') + '\n</svg>';
}
unwrapped_view = html`
<div style="max-width:540px; margin:8px auto; font-family:system-ui,sans-serif;">
${buildUnwrappedPlot("Polar angle θ vs X", "θ (rad)",
trueAngles.map(function(a) { return a.theta; }),
fittedAngles.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; }),
fittedAngles.map(function(a) { return a.phi; }),
dataAngles.map(function(a) { return a.phi; }),
regData.X, xGridEval, 0, Math.PI/2)}
</div>
</div>
`
// ---- Statistics panel ----
meanX_val2 = regData.X.reduce(function(a,b) { return a+b; }, 0) / regData.X.length;
meanX_dist = Math.abs(eval_x - meanX_val2);
statistics_view = html`
<div style="font-family:system-ui,sans-serif; max-width:540px; margin:10px auto; font-size:0.88em;">
<div style="background:#f8f9fa; border-radius:6px; padding:12px 16px;">
<table style="width:100%; border-collapse:collapse;">
<tr style="border-bottom:1px solid #dee2e6;">
<td style="padding:5px 8px;" colspan="2"><b>Evaluation at x = ${eval_x.toFixed(2)}</b> (X̄ = ${meanX_val2.toFixed(3)}, distance from X̄ = ${meanX_dist.toFixed(3)})</td>
</tr>
<tr>
<td style="padding:4px 8px;">Fitted μ̂(x)</td>
<td style="padding:4px 8px; text-align:right;">θ=${evalFitAng.theta.toFixed(3)}, φ=${evalFitAng.phi.toFixed(3)} rad</td>
</tr>
<tr>
<td style="padding:4px 8px;">True μ(x)</td>
<td style="padding:4px 8px; text-align:right;">θ=${evalTrueAng.theta.toFixed(3)}, φ=${evalTrueAng.phi.toFixed(3)} rad</td>
</tr>
<tr>
<td style="padding:4px 8px;">Geodesic error</td>
<td style="padding:4px 8px; text-align:right; font-weight:bold; color:#c62828;">${errorAtEval.toFixed(5)} rad (${(errorAtEval*180/Math.PI).toFixed(3)}°)</td>
</tr>
<tr>
<td style="padding:4px 8px;">Weight sum at eval x</td>
<td style="padding:4px 8px; text-align:right;">${weightSum.toFixed(3)} / ${regData.X.length} (mean = ${(weightSum/regData.X.length).toFixed(3)} ≈ 1)</td>
</tr>
<tr>
<td style="padding:4px 8px;">Negative weights at x</td>
<td style="padding:4px 8px; text-align:right;">${negCount} of ${regData.X.length} (${(negCount/regData.X.length*100).toFixed(1)}%)</td>
</tr>
<tr>
<td style="padding:4px 8px;">Global Fréchet mean μ̂(X̄)</td>
<td style="padding:4px 8px; text-align:right;">(${overallMean[0].toFixed(3)}, ${overallMean[1].toFixed(3)}, ${overallMean[2].toFixed(3)})</td>
</tr>
</table>
</div>
</div>
`
html`<div>${controls_view}${sphere_view}${unwrapped_view}${statistics_view}</div>`
```
::: {.callout-tip title="Try these experiments"}
- **Switch regression types** using the dropdown: "Varying θ only" and "Varying φ only" show how the regression captures trends along a single spherical coordinate, while "Diagonal" constrains both angles to move together.
- **Increase noise $\sigma$** (0.3–0.4): The fitted curve shrinks toward the overall Fréchet mean (gold diamond) — this is the metric analogue of regression-to-the-mean.
- **Evaluate near the edges** ($x = \pm 1$): Watch the "Negative weights" row — as $x$ moves away from $\bar{X}$, more weights become negative, reflecting the extrapolatory nature of the global linear model.
- **Increase sample size $n$** (200+): The fitted (red) curve converges to the true (green dashed) curve. The error at the evaluation point decreases.
- **Change the random seed**: See different data configurations — some produce tight clusters where regression is easy, others are more diffuse.
- **Look at the 2D unwrapped views**: Compare how well the fitted $\theta(x)$ and $\phi(x)$ track the true curves. Mismatches in one coordinate explain the geodesic error seen in the 3D view.
- **Observe the effect of global weights**: Unlike kernel regression, global Fréchet regression uses weights that can go negative and do not decay with distance — all data points influence the fit everywhere.
:::
## Application: Market Risk Structure as a Function of Stress
Lecture 4 asked a **categorical** question: do calm, normal, and stress market regimes have different Fréchet mean covariance matrices? Global Fréchet regression lets us ask a richer, **continuous** question: how does the typical covariance structure evolve as market stress increases?
### From Regimes to Regression
In the Fréchet ANOVA application, we discretized market conditions into three bins (calm / normal / stress) by thresholding realized SPY volatility. This raises a natural follow-up: rather than splitting a continuous stress variable into arbitrary categories, why not model the relationship directly?
Global Fréchet regression provides exactly this tool. We set:
- **Predictor** $X_t = \hat{\sigma}_t$, the realized SPY volatility over the trailing window (a continuous scalar measuring market stress).
- **Response** $Y_t = \Sigma_t \in \mathcal{S}_{++}^p$, the rolling covariance matrix of the six-asset universe (SPY, TLT, GLD, XLF, XLE, QQQ), computed over the same window.
The same data that were assigned to discrete regimes in Lecture 4 can now be used to fit a continuous regression curve $\mu(x)$ on the SPD manifold.
### Global Fréchet Regression on SPD Matrices
Under the **log-Euclidean metric** $d_{\mathrm{LE}}(A, B) = \|\log A - \log B\|_F$, the population regression function at stress level $x$ solves
$$
\mu(x) = \arg\min_{\omega \in \mathcal{S}_{++}^p} \mathbb{E}\bigl\{s(X, x)\, \|\log Y - \log \omega\|_F^2\bigr\},
$$
where $s(X, x) = 1 + (X - \mathbb{E}X)(x - \mathbb{E}X)/\operatorname{var}(X)$ is the scalar weight function for a one-dimensional predictor. Because the log-Euclidean metric makes $\mathcal{S}_{++}^p$ isometric to a Euclidean space of symmetric matrices, the minimizer has the **closed form**
$$
\log \hat{\mu}(x) = \frac{1}{n}\sum_{i=1}^n s_{in}(x)\, \log \Sigma_i,
\qquad
\hat{\mu}(x) = \exp\!\bigl(\log\hat{\mu}(x)\bigr),
$$
where $s_{in}(x) = 1 + (X_i - \bar{X})(x - \bar{X}) / \hat{\sigma}_X^2$. Each $\log \Sigma_i$ is a symmetric matrix; the regression estimator simply computes **weighted matrix averages in the log domain** and exponentiates the result back to the SPD manifold.
### What the Regression Reveals
The global Fréchet regression curve $\hat{\mu}(x)$ traces a path through the SPD manifold parameterized by market stress $x$:
- **At low stress** ($x$ small), $\hat{\mu}(x)$ should resemble the calm-regime Fréchet mean from the ANOVA — low volatilities, moderate correlations.
- **At high stress** ($x$ large), $\hat{\mu}(x)$ moves toward the stress-regime mean — elevated volatilities, heightened equity correlations.
- **At $\bar{X}$** (average stress), the curve passes through the **unconditional Fréchet mean** of all covariance matrices, exactly as predicted by the theory.
Crucially, the regression quantifies the **rate of change** of the typical covariance structure with respect to market stress. If $\hat{\mu}(x)$ changes rapidly in a particular region of $x$, this indicates a nonlinear transition in market risk structure — information that is lost when stress is discretized into coarse bins.
### Portfolio-Risk Interpretation
For a portfolio weight vector $w \in \mathbb{R}^p$, the implied portfolio variance at stress level $x$ is
$$
\sigma_P^2(x) = w^\top \hat{\mu}(x) w.
$$
The global Fréchet regression therefore provides a **stress-dependent risk model**: rather than using a single covariance matrix for all conditions, one can evaluate portfolio risk conditional on the current stress level. The difference $\sigma_P^2(x_{\text{high}}) - \sigma_P^2(x_{\text{low}})$ measures how much additional portfolio variance is expected when markets transition from calm to stressed conditions — a quantity of direct interest for risk managers.
::: {.callout-important title="From ANOVA to regression"}
Fréchet ANOVA tests whether group means differ. Global Fréchet regression models **how** the mean varies with a continuous predictor. Together, they form a natural progression: first test for any difference (ANOVA), then model the functional form (regression). Both operate on the same SPD-valued data with the same log-Euclidean geometry — only the question changes.
:::
### Takeaway
Global Fréchet regression turns the discrete regime comparison of Lecture 4 into a continuous model of how market risk structure evolves with stress. The same six-asset covariance matrices, the same log-Euclidean geometry, and the same computational primitives (weighted Fréchet means) serve both analyses. The regression adds the ability to quantify rates of change and to interpolate risk estimates at any stress level, not just at three discrete bins.
## Key Takeaways
- **Global Fréchet regression** extends ordinary multiple linear regression to metric-space responses by rewriting linear regression as a weighted least-squares problem and replacing Euclidean distance with metric distance.
- The **weight function** $s(z, x) = 1 + (z - \mathbb{E}X)^\top\Sigma^{-1}(x - \mathbb{E}X)$ absorbs all linear structure into the predictor side; the response enters only through squared distances.
- At $x = \mathbb{E}X$, the regression curve passes through the **unconditional Fréchet mean** — recovering the classical property that the regression line goes through $(\bar{X}, \bar{Y})$.
- The **sample estimator** is a weighted Fréchet mean, reducing regression to a computational problem already solved for Fréchet means.
- **Consistency** requires a well-separated minimum (P0), and the **convergence rate** $O_p(n^{-1/(2\beta-2)})$ depends on the curvature exponent $\beta$ of the population objective near its minimum.
## Exercises
1. **Weight properties.** Show that for the global Fréchet regression weight function, $\mathbb{E}[s(X, x)] = 1$ for all $x$. What is $\operatorname{var}(s(X, x))$, and how does it behave as $x$ moves away from $\mathbb{E}X$? <a href="javascript:void(0)" onclick="showSolution('l6-sol-1')" class="solution-link">📝 Show Solution</a>
2. **Recovery of Euclidean linear regression.** Verify that when $\mathcal{M} = \mathbb{R}^q$ with the Euclidean metric, the global Fréchet regression minimizer is exactly the usual multivariate linear regression fit componentwise. <a href="javascript:void(0)" onclick="showSolution('l6-sol-2')" class="solution-link">📝 Show Solution</a>
3. **Curvature exponent $\beta$ for SPD spaces.** Explain why $\beta = 2$ for SPD spaces with the log-Euclidean metric. (Hint: consider the second-order Taylor expansion of the Fréchet function near its minimum.) <a href="javascript:void(0)" onclick="showSolution('l6-sol-3')" class="solution-link">📝 Show Solution</a>
4. **Impact of predictor dimension.** The weights involve $\hat{\Sigma}^{-1}$. Discuss the practical challenges when $p$ (the predictor dimension) is large relative to $n$. How might one regularize the estimator? <a href="javascript:void(0)" onclick="showSolution('l6-sol-4')" class="solution-link">📝 Show Solution</a>
<style>
.solution-link {
font-size: 0.9em;
text-decoration: none;
white-space: nowrap;
margin-left: 0.3em;
}
.solution-link:hover {
text-decoration: underline;
}
.solution-dialog {
padding: 0;
max-width: 720px;
}
.solution-dialog-header {
display: flex;
justify-content: space-between;
align-items: flex-start;
border-bottom: 1px solid #dee2e6;
padding: 1.25rem 1.5rem 1rem;
background: #f8f9fa;
border-radius: 8px 8px 0 0;
}
.solution-dialog-header h4 {
margin: 0;
font-size: 1.15rem;
}
.solution-dialog-close {
background: none;
border: 1px solid #adb5bd;
border-radius: 4px;
padding: 0.2rem 0.75rem;
cursor: pointer;
font-size: 0.9rem;
color: #495057;
white-space: nowrap;
flex-shrink: 0;
}
.solution-dialog-close:hover {
background: #e9ecef;
}
.solution-original {
padding: 1rem 1.5rem;
background: #f1f3f5;
border-left: 4px solid #868e96;
margin: 1rem 1.5rem;
border-radius: 4px;
font-size: 0.95rem;
}
.solution-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="l6-sol-1">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 1: Weight Properties</h4>
<button onclick="closeSolution('l6-sol-1')" class="solution-dialog-close">✕ Close</button>
</div>
<div class="solution-original">
<strong>Exercise:</strong> Show that $\mathbb{E}[s(X, x)] = 1$ for all $x$. What is $\operatorname{var}(s(X, x))$, and how does it behave as $x$ moves away from $\mathbb{E}X$?
</div>
<div class="solution-answer">
<strong>Solution:</strong>
Since $\mathbb{E}[(X - \mathbb{E}X)] = 0$, we have
$$
\mathbb{E}[s(X, x)] = \mathbb{E}[1 + (X - \mathbb{E}X)^\top\Sigma^{-1}(x - \mathbb{E}X)] = 1 + \mathbb{E}[(X - \mathbb{E}X)]^\top\Sigma^{-1}(x - \mathbb{E}X) = 1.
$$
For the variance, note that $s(X, x) = 1 + a^\top(X - \mathbb{E}X)$ where $a = \Sigma^{-1}(x - \mathbb{E}X)$. Then
$$
\operatorname{var}(s(X, x)) = \operatorname{var}(a^\top(X - \mathbb{E}X)) = a^\top\Sigma a = (x - \mathbb{E}X)^\top\Sigma^{-1}(x - \mathbb{E}X).
$$
This is the **Mahalanobis distance** between $x$ and $\mathbb{E}X$. It grows quadratically as $x$ moves away from $\mathbb{E}X$ — so the weights become more variable, and some may become negative, reflecting the extrapolatory nature of the global linear model.
</div>
</div>
</dialog>
<dialog id="l6-sol-2">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 2: Recovery of Euclidean Linear Regression</h4>
<button onclick="closeSolution('l6-sol-2')" class="solution-dialog-close">✕ Close</button>
</div>
<div class="solution-original">
<strong>Exercise:</strong> Verify that when $\mathcal{M} = \mathbb{R}^q$ with the Euclidean metric, the global Fréchet regression minimizer is exactly the usual multivariate linear regression fit componentwise.
</div>
<div class="solution-answer">
<strong>Solution:</strong>
For $\mathcal{M} = \mathbb{R}^q$ with $d(y, \omega) = \|y - \omega\|_2$, the objective is
$$
M(\omega, x) = \mathbb{E}\{s(X, x)\|Y - \omega\|_2^2\} = \sum_{j=1}^q \mathbb{E}\{s(X, x)(Y_j - \omega_j)^2\}.
$$
This separates componentwise. For component $j$, differentiate with respect to $\omega_j$:
$$
\frac{\partial}{\partial \omega_j} M(\omega, x) = -2\,\mathbb{E}\{s(X, x)(Y_j - \omega_j)\} = 0.
$$
Since $\mathbb{E}[s(X, x)] = 1$,
$$
\omega_j = \mathbb{E}[s(X, x)Y_j] = \mathbb{E}Y_j + \operatorname{Cov}(X, Y_j)^\top\Sigma^{-1}(x - \mathbb{E}X),
$$
which is exactly the $j$-th component of the multivariate linear regression of $Y$ on $X$. Thus the global Fréchet regression recovers the usual linear model componentwise.
</div>
</div>
</dialog>
<dialog id="l6-sol-3">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 3: Curvature Exponent $\beta$ for SPD Spaces</h4>
<button onclick="closeSolution('l6-sol-3')" class="solution-dialog-close">✕ Close</button>
</div>
<div class="solution-original">
<strong>Exercise:</strong> Explain why $\beta = 2$ for SPD spaces with the log-Euclidean metric.
</div>
<div class="solution-answer">
<strong>Solution:</strong>
Under the log-Euclidean metric, the map $\Sigma \mapsto \log\Sigma$ is an isometric bijection from $(\mathcal{S}_{++}^m, d_{\mathrm{LE}})$ to the Euclidean space of symmetric matrices. After applying the matrix logarithm, the Fréchet function becomes
$$
F(\omega) = \mathbb{E}\|\log Y - \log \omega\|_F^2,
$$
which is a **quadratic** function in $\log\omega$. The minimizer is $\log\mu = \mathbb{E}[\log Y]$, and a Taylor expansion around $\mu$ gives
$$
F(\omega) - F(\mu) = \|\log\omega - \log\mu\|_F^2 = d_{\mathrm{LE}}^2(\omega, \mu).
$$
Thus $F(\omega) - F(\mu) \ge C\, d^2(\omega, \mu)$ with $C = 1$ and $\beta = 2$. The same argument works for any metric space that is isometric to a convex subset of a Hilbert space — the squared distance is exactly quadratic.
</div>
</div>
</dialog>
<dialog id="l6-sol-4">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 4: High-Dimensional Predictors</h4>
<button onclick="closeSolution('l6-sol-4')" class="solution-dialog-close">✕ Close</button>
</div>
<div class="solution-original">
<strong>Exercise:</strong> Discuss the practical challenges when $p$ is large relative to $n$. How might one regularize the estimator?
</div>
<div class="solution-answer">
<strong>Solution:</strong>
When $p$ is large relative to $n$:
1. **$\hat{\Sigma}$ is singular or ill-conditioned.** The sample covariance matrix is not invertible when $p > n$, making the weights $s_{in}(x)$ undefined.
2. **Variance inflation.** Even when $p < n$, the Mahalanobis distance $(x - \bar{X})^\top\hat{\Sigma}^{-1}(x - \bar{X})$ can be very large, leading to extreme weights.
**Regularization strategies:**
- **Ridge-type regularization:** Replace $\hat{\Sigma}^{-1}$ with $(\hat{\Sigma} + \lambda I)^{-1}$. This shrinks the weights toward 1 (the equal-weight case that recovers the ordinary Fréchet mean).
- **Dimension reduction:** Project $X$ onto a lower-dimensional subspace (e.g., via PCA) before applying global Fréchet regression.
- **Sparse regularization:** If only a few predictors are relevant, use a regularized covariance estimator (e.g., graphical lasso) for $\hat{\Sigma}^{-1}$.
- **Switching to local methods:** Kernel or local Fréchet regression avoids the global covariance inversion entirely; their weights are always nonnegative and decay away from $x$.
These strategies trade bias for variance, just as in classical high-dimensional linear regression.
</div>
</div>
</dialog>
<dialog id="l6-sol-5">
<div class="solution-dialog">
<div class="solution-dialog-header">
<h4>Exercise 5: Global vs. Local Methods</h4>
<button onclick="closeSolution('l6-sol-5')" class="solution-dialog-close">✕ Close</button>
</div>
<div class="solution-original">
<strong>Exercise:</strong> Under what circumstances would you prefer global over local (kernel) Fréchet regression, and vice versa?
</div>
<div class="solution-answer">
<strong>Solution:</strong>
**Prefer global Fréchet regression when:**
- The true regression function is approximately linear (in the sense of the weight representation) — i.e., the conditional Fréchet mean varies linearly with $x$ in the intrinsic geometry.
- The predictor dimension $p$ is small to moderate relative to $n$.
- You need interpretability: the global model gives a single "slope" measure of association between $X$ and $Y$, analogous to a regression coefficient.
- You need to extrapolate beyond the range of the observed $X$ values.
**Prefer local (kernel) Fréchet regression when:**
- The relationship is nonlinear — kernel regression makes no global parametric assumption.
- $p$ is large — local methods avoid inverting a large covariance matrix.
- You have dense data and can afford a small bandwidth.
- Interpretability of a global linear trend is not required; prediction accuracy is the goal.
**Hybrid approach:** Local Fréchet regression (Lecture 8) combines the flexibility of kernel smoothing with the bias-reduction of local linear fitting — it is often the best default choice in practice.
</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 6 Quiz — 10 Multiple-Choice Questions](../quizzes/lecture-06-quiz.qmd)**
## Further Reading
- @PetersenMueller2019 — The foundational paper establishing global, local, and kernel Fréchet regression.
- @vanDerVaartWellner1996 — Empirical process theory underlying the argmin theorem and the convergence-rate analysis (Chapters 3.2, 2.7, 2.14).
- @LinMuller2021TV — Total-variation regularized Fréchet regression for data with jumps and changepoints.
- @DavisFletcherBullittJoshi2007 — Early work on manifold-valued kernel regression for shape data.