Lecture 5: Fréchet Means — Computation

Proximal algorithms for Fréchet means in Hadamard spaces

1 Learning Goals

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

Define the proximal map (resolvent) of a convex function on a Hadamard space.
Describe the proximal-point algorithm (PPA) and explain why each step is well defined.
Derive the closed-form component resolvent for the Fréchet mean objective.
Formulate both cyclic and randomized split proximal-point algorithms for computing sample Fréchet means.
State the convergence theorem for the PPA in locally compact Hadamard spaces.
Implement a basic PPA for computing Fréchet means with geodesic interpolation.

2 Why Computation Is Geometry-Dependent

The computational side of Fréchet means is equally geometry-dependent as the theoretical side. In a general metric space, computing a Fréchet mean involves minimizing

\[ F_n(x) = \frac{1}{n}\sum_{i=1}^n d^2(x, X_i) \]

over $x \in \mathcal{M}$, which may be non-convex and may have local minima.

However, in locally compact Hadamard spaces, the empirical Fréchet function is convex, and proximal algorithms become available. These algorithms exploit the fact that squared distances in Hadamard spaces are strongly convex along geodesics.

3 Proximal Maps in Hadamard Spaces

The basic variational object is the proximal subproblem.

Definition 1 Let $f : \mathcal{M} \to (-\infty, \infty]$ be a proper lower-semicontinuous convex function on a Hadamard space $(\mathcal{M}, d)$. For $\lambda > 0$ and $x \in \mathcal{M}$, define its proximal map (or resolvent) by

\[ J_{\lambda}^{f}(x) = \arg\min_{y \in \mathcal{M}}\left\{f(y) + \frac{1}{2\lambda}d^2(y, x)\right\}. \]

Because the squared-distance term is strongly convex in a Hadamard space, the minimizer is unique.

The proximal map generalizes the Euclidean proximal operator $J_{\lambda}^f(x) = \arg\min_y\{f(y) + \frac{1}{2\lambda}\|y-x\|_2^2\}$ to Hadamard spaces. The term $\frac{1}{2\lambda}d^2(y, x)$ penalizes movement away from $x$, so each step balances descent of $f$ against excessive movement.

4 The Proximal-Point Algorithm (PPA)

The proximal-point algorithm (PPA) is the iteration

\[ x_{k+1} = J_{\lambda_k}^{f}(x_k), \qquad k = 0, 1, 2, \ldots, \]

started from an arbitrary $x_0 \in \mathcal{M}$. Intuitively, the iterate $x_{k+1}$ minimizes the original objective $f$ penalized by the squared distance to the previous iterate.

For the weighted sample Fréchet mean problem,

\[ f(x) = \sum_{i=1}^n w_i\,d^2(x, X_i), \qquad w_i > 0, \qquad \sum_{i=1}^n w_i = 1, \]

it is natural to split the objective into summands:

\[ f = \sum_{i=1}^n f_i, \qquad f_i(x) = w_i\,d^2(x, X_i). \]

Note that taking $w_i = 1/n$ for all $i$ corresponds to computing the ordinary sample Fréchet mean.

The full resolvent $J_{\lambda}^{f}$ is globally well defined, but evaluating it may be almost as hard as minimizing $f$ directly. The practical idea is therefore to apply the proximal step to one summand at a time.

5 Split Proximal-Point Algorithms

There are three closely related formulations.

Full Proximal-Point Iteration

\[ x_{k+1} = J_{\lambda_k}^{f}(x_k). \]

Conceptually simplest, but each step requires solving the whole regularized Fréchet minimization problem — often impractical.

Cyclic Split Proximal-Point Iteration

For $k = 0, 1, \ldots$ (indexing the cycle), set

\[ \begin{aligned} x_{kn+1} &= J_{\lambda_k}^{f_1}(x_{kn}), \\ x_{kn+2} &= J_{\lambda_k}^{f_2}(x_{kn+1}), \\ &\;\;\vdots \\ x_{kn+n} &= J_{\lambda_k}^{f_n}(x_{kn+n-1}). \end{aligned} \]

The step size parameter $\lambda_k$ remains fixed within the $k$th cycle.

Randomized Split Proximal-Point Iteration

Choose an index $r_k \in \{1, \ldots, n\}$ at iteration $k$ according to a prescribed distribution (typically uniform) and set

\[ x_{k+1} = J_{\lambda_k}^{f_{r_k}}(x_k). \]

6 Closed Form of the Component Resolvent

For the Fréchet mean objective, every component resolvent has a closed form. Consider the objective

\[ w\,d^2(\cdot, a) + \frac{1}{2\lambda}d^2(\cdot, y), \]

where $a$ is a data point, $y$ is the current iterate, and $w$ is the weight. The minimizer $z$ lies on the geodesic $[y, a]$ connecting $y$ to $a$, i.e.,

\[ z = (1-t)y + t a, \]

where the “$+$” operation represents geodesic interpolation. A simple calculation shows that the optimal $t$ is

\[ t = \frac{2\lambda w}{1 + 2\lambda w} \in (0, 1). \]

Derivation

On the geodesic, the objective becomes a quadratic function of $t$. Since $d(z, a) = (1-t)d(y, a)$ and $d(z, y) = t\,d(y, a)$ for a constant-speed geodesic, we minimize

\[ w(1-t)^2 d^2(y, a) + \frac{1}{2\lambda}t^2 d^2(y, a). \]

Differentiating with respect to $t$ yields $t = \frac{2\lambda w}{1 + 2\lambda w}$.

Therefore, the update in the cyclic PPA becomes

\[ x_{kn+i} = (1 - t_{k,i})\,x_{kn+i-1} + t_{k,i}\,X_i, \]

where $t_{k,i} = \frac{2\lambda_k w_i}{1 + 2\lambda_k w_i}$. A split proximal step is just a geodesic interpolation toward one data point.

The step-size factor $t_{k,i}$ shows the basic behavior clearly:

For small $\lambda_k$, $t_{k,i} \approx 2\lambda_k w_i$ is small, so the move is conservative.
For large $\lambda_k$, $t_{k,i} \to 1$, so the update moves farther toward the data point $X_i$.

7 Algorithm Summary

A randomized proximal-point scheme for weighted sample Fréchet means:

Step 1: Choose an initial point $x_0 \in \mathcal{M}$ and a step-size sequence $(\lambda_k)$ satisfying

\[ \sum_{k=0}^{\infty} \lambda_k = \infty, \qquad \sum_{k=0}^{\infty} \lambda_k^2 < \infty. \]

Step 2: At iteration $k$, sample an index $r_k \in \{1, \ldots, N\}$ (e.g., uniformly).

Step 3: Compute the geodesic from $x_k$ to $X_{r_k}$ and move along it by

\[ t_k = \frac{2\lambda_k w_{r_k}}{1 + 2\lambda_k w_{r_k}}. \]

Equivalently, set $x_{k+1} = J_{\lambda_k}^{f_{r_k}}(x_k)$.

Step 4: Stop when $d(x_{k+1}, x_k)$ is below a chosen tolerance, or when the objective decrease becomes negligible.

The cyclic counterpart is analogous, cycling through the data points in a fixed order.

Step-Size Conditions

The condition $\sum \lambda_k = \infty$ ensures the algorithm can move arbitrarily far (sufficient exploration), while $\sum \lambda_k^2 < \infty$ ensures the steps eventually become small enough to settle (sufficient damping). A typical choice is $\lambda_k = c/k$ for some $c > 0$.

8 Convergence

The following result specializes Corollary 3.9 of Bačák (2014) to the sample Fréchet objective.

Proposition 1 Let $(\mathcal{M}, d)$ be a locally compact Hadamard space, let $X_1, \ldots, X_n \in \mathcal{M}$, and let $w_i = 1/n$. Consider the weighted sample Fréchet function

\[ f(x) = \sum_{i=1}^n w_i\,d^2(x, X_i). \]

Let $(\lambda_k)$ be a sequence of positive numbers such that

\[ \sum_{k=0}^{\infty} \lambda_k = \infty, \qquad \sum_{k=0}^{\infty} \lambda_k^2 < \infty. \]

Starting from $x_0 \in \mathcal{M}$, the sequences defined by the cyclic and randomized split PPA converge to the sample Fréchet mean as $k \to \infty$ (for the randomized scheme, almost surely).

The proximal-point algorithm is particularly attractive in Hadamard spaces: convexity guarantees a unique barycenter, while geodesic interpolation makes every split proximal step explicit.

9 Interactive Exploration: Proximal Algorithm on the Sphere

The demo below runs the randomized split PPA for computing the Fréchet mean of points on a quadrant of the unit sphere $\mathbb{S}^2$.

Data generation. $n$ points $\{X_1, \ldots, X_n\}$ are sampled on the sphere octant $\{x, y, z \ge 0\}$ clustered around a randomly chosen center $c \in \mathbb{S}^2$: \[X_i = \gamma_{c \leftarrow U_i}(W_i), \qquad U_i \sim \text{Uniform}(\text{octant}), \quad W_i \sim \text{Uniform}(0.4, 1.0),\] where $\gamma_{c \leftarrow u}(w)$ denotes the point at fraction $w$ along the great-circle geodesic from $u$ to $c$. This produces a random cluster in the octant whose Fréchet mean is near $c$.

Algorithm. At each iteration $k$, the PPA picks a random data point $X_{r_k}$ and moves the current iterate $x_k$ along the geodesic toward it: \[x_{k+1} = \gamma_{x_k \to X_{r_k}}(t_k), \qquad t_k = \frac{2\lambda_k / n}{1 + 2\lambda_k / n}, \qquad \lambda_k = \frac{c}{k+1}.\] On the sphere, the geodesic from $p$ to $q$ with angular distance $\theta = \arccos\langle p, q\rangle$ is \[\gamma_{p \to q}(t) = \frac{\sin((1-t)\theta)}{\sin\theta}\,p \;+\; \frac{\sin(t\theta)}{\sin\theta}\,q,\] which stays on the unit sphere for all $t \in [0,1]$. Drag the iteration slider to step through each update and watch how the formulas translate into geometry.

Code

n_pts_control = Inputs.range([3, 25], {step: 1, value: 10, label: "Number of data points"})
c_val_control = Inputs.range([0.1, 5], {step: 0.1, value: 1.0, label: "Step-size constant c"})
data_seed_control = Inputs.range([1, 100], {step: 1, value: 42, label: "Random seed"})
step_control = Inputs.range([0, MAX_ITER], {step: 1, value: 0, label: "Iteration step"})
n_pts = Generators.input(n_pts_control)
c_val = Generators.input(c_val_control)
data_seed = Generators.input(data_seed_control)
step = Generators.input(step_control)

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

function frechetObj(x, dataPts) {
  var n = dataPts.length;
  var sum = 0;
  for (var i = 0; i < n; i++) {
    var d = sphereDist(x, dataPts[i]);
    sum += d * d;
  }
  return sum / n;
}

// ---- Orthographic projection ----
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];
}

function generateSphereData(n, seed) {
  var rng = mulberry32(seed);
  var center = randomOctantPoint(rng);
  var pts = [];
  for (var i = 0; i < n; i++) {
    var randPt = randomOctantPoint(rng);
    var w = 0.4 + rng() * 0.6;
    pts.push(sphereGeodesic(randPt, center, w));
  }
  return pts;
}

// ---- PPA computation ----
MAX_ITER = 200;

function runSpherePPA(dataPts, cVal, maxIter, baseSeed) {
  var n = dataPts.length;
  var weight = 1/n;
  var rng = mulberry32(baseSeed + 9999);
  var x = [1/Math.sqrt(3), 1/Math.sqrt(3), 1/Math.sqrt(3)];
  var trajectory = [{pos: x.slice(), obj: frechetObj(x, dataPts), stepTo: -1}];
  for (var k = 0; k < maxIter; k++) {
    var lambda_k = cVal / (k + 1);
    var idx = Math.floor(rng() * n);
    var a = dataPts[idx];
    var t_k = 2 * lambda_k * weight / (1 + 2 * lambda_k * weight);
    x = sphereGeodesic(x, a, t_k);
    trajectory.push({pos: x.slice(), obj: frechetObj(x, dataPts), stepTo: idx});
  }
  return trajectory;
}

function referenceMean(dataPts, baseSeed) {
  var n = dataPts.length;
  var weight = 1/n;
  var rng = mulberry32(baseSeed + 77777);
  var x = [1/Math.sqrt(3), 1/Math.sqrt(3), 1/Math.sqrt(3)];
  for (var k = 0; k < 3000; k++) {
    var lambda_k = 0.3 / (1 + k/100);
    var idx = Math.floor(rng() * n);
    var a = dataPts[idx];
    var t_k = 2 * lambda_k * weight / (1 + 2 * lambda_k * weight);
    x = sphereGeodesic(x, a, t_k);
  }
  return x;
}

data_pts = generateSphereData(n_pts, data_seed);
trajectory = runSpherePPA(data_pts, c_val, MAX_ITER, data_seed);
ref_mean = referenceMean(data_pts, data_seed);
cur = trajectory[step];

// ---- Two-column slider panel ----
html`
<style>
  .ppa-slider-grid {
    display: grid;
    grid-template-columns: repeat(2, minmax(0, 1fr));
    gap: 6px 20px;
    max-width: 560px;
    margin: 0 auto 10px auto;
    font-family: system-ui, sans-serif;
    font-size: 0.85em;
  }
  .ppa-slider-grid > * {
    min-width: 0;
    margin: 0;
  }
  @media (max-width: 500px) {
    .ppa-slider-grid { grid-template-columns: 1fr; }
  }
</style>
<div class="ppa-slider-grid">
  <div>${n_pts_control}</div>
  <div>${c_val_control}</div>
  <div>${data_seed_control}</div>
  <div>${step_control}</div>
</div>
`

Code

function buildSVG(dataPts, trajectory, stepIdx, refMean) {
  var curState = trajectory[stepIdx];
  var curPos = curState.pos;

  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.25, 0.5, 0.75].forEach(function(zv) {
    gridLines.push(constantZArc(zv, 30));
  });

  var dataProj = dataPts.map(function(p) { return project(p); });

  var trajProj = [];
  for (var i = 0; i <= stepIdx; i++) {
    trajProj.push(project(trajectory[i].pos));
  }

  var curProj = project(curPos);
  var refProj = project(refMean);
  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.6"/>');
  });

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

  // Trajectory path
  if (trajProj.length > 1) {
    var trajPts = trajProj.map(function(p) { return p[0].toFixed(1) + "," + p[1].toFixed(1); }).join(" ");
    parts.push('<polyline points="' + trajPts + '" fill="none" stroke="#37474f" stroke-width="2.2" stroke-linecap="round" stroke-linejoin="round" opacity="0.75"/>');
  }

  // Geodesic toward next data point
  if (curState.stepTo >= 0 && stepIdx < trajectory.length - 1) {
    var targetPt = dataPts[curState.stepTo];
    var targetProj = project(targetPt);
    parts.push('<line x1="' + curProj[0].toFixed(1) + '" y1="' + curProj[1].toFixed(1) +
      '" x2="' + targetProj[0].toFixed(1) + '" y2="' + targetProj[1].toFixed(1) +
      '" stroke="#ff9800" stroke-width="1.5" stroke-dasharray="5,3" opacity="0.7"/>');
    parts.push('<circle cx="' + targetProj[0].toFixed(1) + '" cy="' + targetProj[1].toFixed(1) +
      '" r="6" fill="none" stroke="#ff9800" stroke-width="2.2"/>');
  }

  // Data points
  dataProj.forEach(function(dp, idx) {
    var isTarget = curState.stepTo === idx;
    var r = isTarget ? 4.5 : 3.3;
    var fill = isTarget ? "#ff9800" : "#e53935";
    var stroke = isTarget ? "#e65100" : "#b71c1c";
    parts.push('<circle cx="' + dp[0].toFixed(1) + '" cy="' + dp[1].toFixed(1) +
      '" r="' + r + '" fill="' + fill + '" stroke="' + stroke + '" stroke-width="0.8"/>');
    parts.push('<title>Data point ' + (idx+1) + '</title>');
  });

  // Current iterate
  parts.push('<circle cx="' + curProj[0].toFixed(1) + '" cy="' + curProj[1].toFixed(1) +
    '" r="8.5" fill="#43a047" stroke="#fff" stroke-width="2.5"/>');
  parts.push('<circle cx="' + curProj[0].toFixed(1) + '" cy="' + curProj[1].toFixed(1) +
    '" r="8.5" fill="none" stroke="#43a047" stroke-width="1" opacity="0.35">' +
    '<animate attributeName="r" values="8.5;15;8.5" dur="2s" repeatCount="indefinite"/></circle>');

  // Reference Fréchet mean
  parts.push('<polygon points="' +
    (refProj[0]).toFixed(1) + ',' + (refProj[1]-7).toFixed(1) + ' ' +
    (refProj[0]+5.5).toFixed(1) + ',' + refProj[1].toFixed(1) + ' ' +
    refProj[0].toFixed(1) + ',' + (refProj[1]+7).toFixed(1) + ' ' +
    (refProj[0]-5.5).toFixed(1) + ',' + refProj[1].toFixed(1) +
    '" fill="#f9a825" stroke="#f57f17" stroke-width="1.3"/>');

  // 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 = 52;
  parts.push('<circle cx="' + lx + '" cy="' + ly0 + '" r="4" fill="#e53935"/>');
  parts.push('<text x="' + (lx+10) + '" y="' + (ly0+4) + '" font-size="10" fill="#37474f">Data point</text>');
  parts.push('<circle cx="' + lx + '" cy="' + (ly0+17) + '" r="6" fill="#43a047" stroke="#fff" stroke-width="1.2"/>');
  parts.push('<text x="' + (lx+10) + '" y="' + (ly0+21) + '" font-size="10" fill="#37474f">Current iterate</text>');
  parts.push('<polygon points="' + lx + ',' + (ly0+31) + ' ' + (lx+4.5) + ',' + (ly0+37) +
    ' ' + lx + ',' + (ly0+43) + ' ' + (lx-4.5) + ',' + (ly0+37) + '" fill="#f9a825"/>');
  parts.push('<text x="' + (lx+10) + '" y="' + (ly0+41) + '" font-size="10" fill="#37474f">Ref. Fréchet mean</text>');
  parts.push('<line x1="' + lx + '" y1="' + (ly0+57) + '" x2="' + (lx+16) + '" y2="' + (ly0+57) +
    '" stroke="#ff9800" stroke-width="1.2" stroke-dasharray="3,2"/>');
  parts.push('<text x="' + (lx+22) + '" y="' + (ly0+61) + '" font-size="10" fill="#37474f">Geodesic step</text>');

  return '<svg viewBox="0 0 500 420" xmlns="http://www.w3.org/2000/svg" style="width:100%;max-width:520px;">\n' +
    parts.join('\n') + '\n</svg>';
}

html`<div style="text-align:center;">${buildSVG(data_pts, trajectory, step, ref_mean)}</div>`

Code

n_weight = 1 / data_pts.length;
cw = cur.stepTo >= 0 ? data_pts[cur.stepTo] : null;
lam_k = step > 0 ? c_val / step : null;
t_k_val = lam_k !== null ? (2 * lam_k * n_weight) / (1 + 2 * lam_k * n_weight) : null;
geo_theta = cw !== null ? sphereDist(cur.pos, cw) : null;
geo_w0 = geo_theta !== null && geo_theta > 1e-12 ? Math.sin((1 - t_k_val) * geo_theta) / Math.sin(geo_theta) : null;
geo_w1 = geo_theta !== null && geo_theta > 1e-12 ? Math.sin(t_k_val * geo_theta) / Math.sin(geo_theta) : null;
next_from_geo = cw !== null && geo_w0 !== null ? [
  geo_w0 * cur.pos[0] + geo_w1 * cw[0],
  geo_w0 * cur.pos[1] + geo_w1 * cw[1],
  geo_w0 * cur.pos[2] + geo_w1 * cw[2]
] : null;

html`
<div style="font-family:system-ui,sans-serif; max-width:520px; margin:0 auto; font-size:0.88em;">

  <div style="display:flex; justify-content:space-between; align-items:center;
    background:#f5f5f5; border-radius:6px; padding:8px 14px; margin-bottom:10px;">
    <span><b>Step</b> ${step} / ${MAX_ITER}</span>
    <span style="color:#e53935;">→ data point <b>${cur.stepTo >= 0 ? "#" + (cur.stepTo + 1) : "—"}</b></span>
    <span><b>F<sub>n</sub></b> = ${cur.obj.toFixed(5)}</span>
  </div>

  ${step === 0 ? html`
  <div style="background:#fff3e0; border-left:3px solid #ff9800; border-radius:4px; padding:10px 14px; margin-bottom:10px;">
    <b>Initialization:</b> x<sub>0</sub> = (1/√3, 1/√3, 1/√3) — the center of the octant.<br>
    Drag the slider to step through the PPA iterations.
  </div>`
  : (lam_k !== null && cw !== null ? html`
  <div style="display:grid; grid-template-columns:1fr 1fr; gap:10px 16px;">

    <div style="background:#fafafa; border-radius:6px; padding:10px 14px; border:1px solid #e0e0e0;">
      <div style="font-weight:600; margin-bottom:6px; color:#1565c0;">⓵ Step size</div>
      <div style="margin-bottom:4px;">
        λ<sub>${step}</sub> = c/k = ${c_val.toFixed(1)}/${step} = <b>${lam_k.toFixed(4)}</b>
      </div>
      <div>
        t<sub>${step}</sub> = <span style="font-size:0.82em;">2λw/(1+2λw)</span>
        = <span style="font-size:0.82em;">2·${lam_k.toFixed(4)}·${n_weight.toFixed(4)}/(1+2·${lam_k.toFixed(4)}·${n_weight.toFixed(4)})</span><br>
        = <b>${t_k_val.toFixed(5)}</b>
      </div>
    </div>

    <div style="background:#fafafa; border-radius:6px; padding:10px 14px; border:1px solid #e0e0e0;">
      <div style="font-weight:600; margin-bottom:6px; color:#1565c0;">⓶ Geodesic on 𝕊²</div>
      <div style="margin-bottom:4px;">
        θ = arccos⟨x<sub>${step-1}</sub>, X<sub>${cur.stepTo+1}</sub>⟩ = <b>${geo_theta.toFixed(4)}</b> rad
      </div>
      <div style="font-size:0.82em;">
        w<sub>0</sub> = sin((1−t)θ)/sinθ = <b>${geo_w0.toFixed(5)}</b><br>
        w<sub>1</sub> = sin(tθ)/sinθ = <b>${geo_w1.toFixed(5)}</b>
      </div>
    </div>

    <div style="grid-column:1/-1; background:#e8f5e9; border-radius:6px; padding:10px 14px; border:1px solid #a5d6a7;">
      <div style="font-weight:600; margin-bottom:4px; color:#2e7d32;">⓷ PPA Update</div>
      <div style="font-size:0.9em;">
        x<sub>${step}</sub> = w<sub>0</sub>·x<sub>${step-1}</sub> + w<sub>1</sub>·X<sub>${cur.stepTo+1}</sub>
      </div>
      <div style="font-size:0.78em; color:#555; margin-top:2px;">
        = ${geo_w0.toFixed(4)}·(${cur.pos[0].toFixed(3)}, ${cur.pos[1].toFixed(3)}, ${cur.pos[2].toFixed(3)})<br>
        + ${geo_w1.toFixed(4)}·(${cw[0].toFixed(3)}, ${cw[1].toFixed(3)}, ${cw[2].toFixed(3)})
      </div>
      ${next_from_geo ? html`
      <div style="font-size:0.82em; margin-top:3px; font-weight:600; color:#2e7d32;">
        = (${next_from_geo[0].toFixed(4)}, ${next_from_geo[1].toFixed(4)}, ${next_from_geo[2].toFixed(4)})  [✓ on 𝕊²]
      </div>` : ""}
    </div>

    <div style="grid-column:1/-1; background:#f5f5f5; border-radius:6px; padding:8px 14px;">
      <span><b>‖x<sub>${step}</sub> − μ<sub>ref</sub>‖</b> = ${sphereDist(cur.pos, ref_mean).toFixed(5)} rad</span>
      <span style="margin-left:16px;"><b>F<sub>n</sub>(x<sub>${step}</sub>)</b> = ${cur.obj.toFixed(5)}</span>
    </div>

  </div>` : "")
  }
</div>
`

Code

function buildConvergencePlot(trajectory, stepIdx, maxIter) {
  var W = 500, H = 180;
  var margin = {top: 15, right: 15, bottom: 28, left: 48};
  var pw = W - margin.left - margin.right;
  var ph = H - margin.top - margin.bottom;

  var objs = trajectory.map(function(t) { return t.obj; });
  var yMin = Math.min.apply(null, objs) * 0.95;
  var yMax = Math.max.apply(null, objs) * 1.02;
  var yRange = yMax - yMin || 0.001;

  function xScale(i) { return margin.left + (i / maxIter) * pw; }
  function yScale(v) { return margin.top + ph - ((v - yMin) / yRange) * 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 <= 3; gy++) {
    var yv = yMin + (gy/3) * yRange;
    var yp = yScale(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="9" fill="#9e9e9e">' + yv.toFixed(4) + '</text>');
  }

  var linePts = "";
  for (var i = 0; i < trajectory.length; i++) {
    linePts += xScale(i).toFixed(1) + "," + yScale(trajectory[i].obj).toFixed(1) + " ";
  }
  parts.push('<polyline points="' + linePts +
    '" fill="none" stroke="#1976d2" stroke-width="1.8" stroke-linejoin="round"/>');

  var cx = xScale(stepIdx);
  var cy = yScale(trajectory[stepIdx].obj);
  parts.push('<line x1="' + cx + '" y1="' + margin.top + '" x2="' + cx +
    '" y2="' + (margin.top+ph) + '" stroke="#c92a2a" stroke-width="1.5" stroke-dasharray="4,3"/>');
  parts.push('<circle cx="' + cx + '" cy="' + cy + '" r="4.5" fill="#c92a2a" stroke="#fff" stroke-width="1.5"/>');

  var sx0 = xScale(0), sy0 = yScale(trajectory[0].obj);
  parts.push('<circle cx="' + sx0 + '" cy="' + sy0 + '" r="2.5" fill="#78909c"/>');

  var lastIdx = trajectory.length - 1;
  var sxf = xScale(lastIdx), syf = yScale(trajectory[lastIdx].obj);
  parts.push('<circle cx="' + sxf + '" cy="' + syf + '" r="2.5" fill="#2e7d32"/>');

  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="11" fill="#616161">Iteration</text>');
  parts.push('<text x="' + (margin.left-38) + '" y="' + (margin.top+ph/2) +
    '" text-anchor="middle" font-size="11" fill="#616161" ' +
    'transform="rotate(-90,' + (margin.left-38) + ',' + (margin.top+ph/2) + ')">F_n(x_k)</text>');

  parts.push('<text x="' + cx + '" y="' + (margin.top-2) +
    '" text-anchor="middle" font-size="10" fill="#c92a2a" font-weight="bold">k=' + stepIdx + '</text>');

  return '<svg viewBox="0 0 ' + W + ' ' + H +
    '" xmlns="http://www.w3.org/2000/svg" style="width:100%;max-width:520px;">\n' +
    parts.join('\n') + '\n</svg>';
}

html`<div style="text-align:center; margin-top:8px;">${buildConvergencePlot(trajectory, step, MAX_ITER)}</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: PPA on a sphere quadrant with step-by-step visualization

Try these experiments

Drag the iteration slider from 0 to 200 and watch the green iterate trace great-circle arcs across the octant toward the gold reference Fréchet mean.
Small $c$ (e.g., 0.3): conservative steps. The iterate creeps along; convergence is slow but monotonic.
Large $c$ (e.g., 3–5): aggressive early moves. The first few steps jump near individual data points, then the algorithm gradually settles.
Change the random seed to see different data configurations — some are tightly clustered, others spread across the octant.
Increase $n$ to make the Fréchet mean more sharply defined; the trajectory tightens around the reference.
Watch the orange dashed geodesic: it shows the great-circle arc to the data point chosen at the current step. The green iterate moves a fraction $t_k$ along it.
Look at the convergence plot: $F_n(x_k)$ decreases over time, with occasional small bumps when a distant data point is randomly selected.
Why the sphere quadrant? The octant $\{x, y, z \ge 0\}$ is geodesically convex — geodesics between points in the octant stay within it, so the PPA never leaves the region.

10 Key Takeaways

In Hadamard spaces, the empirical Fréchet function is convex, enabling proximal algorithms.
The proximal map $J_{\lambda}^f(x)$ generalizes the Euclidean proximal operator and is uniquely defined in Hadamard spaces.
Split proximal-point algorithms (cyclic or randomized) update one data point at a time, avoiding the full minimization.
Each component resolvent has a closed form: geodesic interpolation toward the data point with step size $t = \frac{2\lambda w}{1 + 2\lambda w}$.
The step-size sequence must satisfy $\sum \lambda_k = \infty$ (sufficient exploration) and $\sum \lambda_k^2 < \infty$ (damping).
Convergence is guaranteed in locally compact Hadamard spaces (Bačák 2014).

11 Exercises

Component resolvent derivation: Derive the closed form $t = \frac{2\lambda w}{1 + 2\lambda w}$ for the component resolvent on a geodesic. 📝 Show Solution
Step-size conditions: Show that $\lambda_k = c/k$ satisfies $\sum \lambda_k = \infty$ and $\sum \lambda_k^2 < \infty$. What happens if $\lambda_k = c/k^{1/2}$ or $\lambda_k = c/k^2$? 📝 Show Solution
Cyclic vs. randomized: Implement both the cyclic and randomized split PPA for points in $\mathbb{R}^2$ (with the Euclidean metric). Compare their convergence empirically. Does one consistently outperform the other? 📝 Show Solution
SPD matrices: Describe how you would implement geodesic interpolation on the SPD matrix space under the affine-invariant metric. What is the formula for the geodesic between $\Sigma_0$ and $\Sigma_1$? 📝 Show Solution
Convergence rate: For the simple case of computing the Euclidean mean of $n$ points, what is the convergence rate of the randomized PPA? How does it compare to the standard stochastic gradient descent? 📝 Show Solution

12 Further Reading

Bačák (2014) — The definitive reference for proximal algorithms in Hadamard spaces.
Sturm (2003) — Foundational treatment of convexity and probability on $\mathrm{CAT}(0)$ spaces.
Pennec et al. (2006) — Practical algorithms for Fréchet means on SPD matrices.
Arsigny et al. (2006); Arsigny et al. (2007) — The log-Euclidean framework, enabling vector-space computation for SPD matrices.

13 Self-Assessment Quiz

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

👉 Lecture 5 Quiz — 10 Multiple-Choice Questions

References

Arsigny, V., P. Fillard, X. Pennec, and N. Ayache. 2006. “Log-Euclidean Metrics for Fast and Simple Calculus on Diffusion Tensors.” Magnetic Resonance in Medicine 56 (2): 411–21.

Arsigny, V., P. Fillard, X. Pennec, and N. Ayache. 2007. “Geometric Means in a Novel Vector Space Structure on Symmetric Positive-Definite Matrices.” SIAM Journal on Matrix Analysis and Applications 29 (1): 328–47. https://doi.org/10.1137/050637996.

Bačák, Miroslav. 2014. “Computing Medians and Means in Hadamard Spaces.” SIAM Journal on Optimization 24 (3): 1542–66. https://doi.org/10.1137/140953393.

Pennec, X., P. Fillard, and N. Ayache. 2006. “A Riemannian Framework for Tensor Computing.” International Journal of Computer Vision 66 (1): 41–66. https://doi.org/10.1007/s11263-005-3222-z.

Sturm, Karl-Theodor. 2003. “Probability Measures on Metric Spaces of Nonpositive Curvature.” In Heat Kernels and Analysis on Manifolds, Graphs, and Metric Spaces, vol. 338. Contemporary Mathematics. American Mathematical Society. https://doi.org/10.1090/conm/338/06080.

--- title: "Lecture 5: Fréchet Means — Computation" subtitle: "Proximal algorithms for Fréchet means in Hadamard 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: - Define the proximal map (resolvent) of a convex function on a Hadamard space. - Describe the proximal-point algorithm (PPA) and explain why each step is well defined. - Derive the closed-form component resolvent for the Fréchet mean objective. - Formulate both cyclic and randomized split proximal-point algorithms for computing sample Fréchet means. - State the convergence theorem for the PPA in locally compact Hadamard spaces. - Implement a basic PPA for computing Fréchet means with geodesic interpolation. ## Why Computation Is Geometry-Dependent The computational side of Fréchet means is equally geometry-dependent as the theoretical side. In a general metric space, computing a Fréchet mean involves minimizing $$ F_n(x) = \frac{1}{n}\sum_{i=1}^n d^2(x, X_i) $$ over $x \in \mathcal{M}$, which may be non-convex and may have local minima. However, in **locally compact Hadamard spaces**, the empirical Fréchet function is **convex**, and proximal algorithms become available. These algorithms exploit the fact that squared distances in Hadamard spaces are strongly convex along geodesics. ## Proximal Maps in Hadamard Spaces The basic variational object is the **proximal subproblem**. ::: {#def-proximal-map-resolvent .definition title="Proximal map (resolvent)"} Let $f : \mathcal{M} \to (-\infty, \infty]$ be a proper lower-semicontinuous convex function on a Hadamard space $(\mathcal{M}, d)$. For $\lambda > 0$ and $x \in \mathcal{M}$, define its **proximal map** (or **resolvent**) by $$ J_{\lambda}^{f}(x) = \arg\min_{y \in \mathcal{M}}\left\{f(y) + \frac{1}{2\lambda}d^2(y, x)\right\}. $$ Because the squared-distance term is strongly convex in a Hadamard space, the minimizer is **unique**. ::: The proximal map generalizes the Euclidean proximal operator $J_{\lambda}^f(x) = \arg\min_y\{f(y) + \frac{1}{2\lambda}\|y-x\|_2^2\}$ to Hadamard spaces. The term $\frac{1}{2\lambda}d^2(y, x)$ penalizes movement away from $x$, so each step balances **descent of $f$** against **excessive movement**. ## The Proximal-Point Algorithm (PPA) {#sec-computation-ppa-lecture5} The **proximal-point algorithm** (PPA) is the iteration $$ x_{k+1} = J_{\lambda_k}^{f}(x_k), \qquad k = 0, 1, 2, \ldots, $$ started from an arbitrary $x_0 \in \mathcal{M}$. Intuitively, the iterate $x_{k+1}$ minimizes the original objective $f$ penalized by the squared distance to the previous iterate. For the weighted sample Fréchet mean problem, $$ f(x) = \sum_{i=1}^n w_i\,d^2(x, X_i), \qquad w_i > 0, \qquad \sum_{i=1}^n w_i = 1, $$ it is natural to **split** the objective into summands: $$ f = \sum_{i=1}^n f_i, \qquad f_i(x) = w_i\,d^2(x, X_i). $$ Note that taking $w_i = 1/n$ for all $i$ corresponds to computing the ordinary sample Fréchet mean. The full resolvent $J_{\lambda}^{f}$ is globally well defined, but evaluating it may be almost as hard as minimizing $f$ directly. The practical idea is therefore to apply the proximal step to **one summand at a time**. ## Split Proximal-Point Algorithms There are three closely related formulations. 1. Full Proximal-Point Iteration $$ x_{k+1} = J_{\lambda_k}^{f}(x_k). $$ Conceptually simplest, but each step requires solving the **whole** regularized Fréchet minimization problem — often impractical. 2. Cyclic Split Proximal-Point Iteration For $k = 0, 1, \ldots$ (indexing the cycle), set $$ \begin{aligned} x_{kn+1} &= J_{\lambda_k}^{f_1}(x_{kn}), \\ x_{kn+2} &= J_{\lambda_k}^{f_2}(x_{kn+1}), \\ &\;\;\vdots \\ x_{kn+n} &= J_{\lambda_k}^{f_n}(x_{kn+n-1}). \end{aligned} $$ The step size parameter $\lambda_k$ remains **fixed** within the $k$th cycle. 3. Randomized Split Proximal-Point Iteration Choose an index $r_k \in \{1, \ldots, n\}$ at iteration $k$ according to a prescribed distribution (typically uniform) and set $$ x_{k+1} = J_{\lambda_k}^{f_{r_k}}(x_k). $$ ## Closed Form of the Component Resolvent For the Fréchet mean objective, every component resolvent has a **closed form**. Consider the objective $$ w\,d^2(\cdot, a) + \frac{1}{2\lambda}d^2(\cdot, y), $$ where $a$ is a data point, $y$ is the current iterate, and $w$ is the weight. The minimizer $z$ lies on the geodesic $[y, a]$ connecting $y$ to $a$, i.e., $$ z = (1-t)y + t a, $$ where the "$+$" operation represents geodesic interpolation. A simple calculation shows that the optimal $t$ is $$ t = \frac{2\lambda w}{1 + 2\lambda w} \in (0, 1). $$ ::: {.callout-note title="Derivation"} On the geodesic, the objective becomes a quadratic function of $t$. Since $d(z, a) = (1-t)d(y, a)$ and $d(z, y) = t\,d(y, a)$ for a constant-speed geodesic, we minimize $$ w(1-t)^2 d^2(y, a) + \frac{1}{2\lambda}t^2 d^2(y, a). $$ Differentiating with respect to $t$ yields $t = \frac{2\lambda w}{1 + 2\lambda w}$. ::: Therefore, the update in the cyclic PPA becomes $$ x_{kn+i} = (1 - t_{k,i})\,x_{kn+i-1} + t_{k,i}\,X_i, $$ where $t_{k,i} = \frac{2\lambda_k w_i}{1 + 2\lambda_k w_i}$. A split proximal step is just a **geodesic interpolation** toward one data point. The step-size factor $t_{k,i}$ shows the basic behavior clearly: - For **small $\lambda_k$**, $t_{k,i} \approx 2\lambda_k w_i$ is small, so the move is **conservative**. - For **large $\lambda_k$**, $t_{k,i} \to 1$, so the update moves **farther** toward the data point $X_i$. ## Algorithm Summary A randomized proximal-point scheme for weighted sample Fréchet means: **Step 1:** Choose an initial point $x_0 \in \mathcal{M}$ and a step-size sequence $(\lambda_k)$ satisfying $$ \sum_{k=0}^{\infty} \lambda_k = \infty, \qquad \sum_{k=0}^{\infty} \lambda_k^2 < \infty. $$ **Step 2:** At iteration $k$, sample an index $r_k \in \{1, \ldots, N\}$ (e.g., uniformly). **Step 3:** Compute the geodesic from $x_k$ to $X_{r_k}$ and move along it by $$ t_k = \frac{2\lambda_k w_{r_k}}{1 + 2\lambda_k w_{r_k}}. $$ Equivalently, set $x_{k+1} = J_{\lambda_k}^{f_{r_k}}(x_k)$. **Step 4:** Stop when $d(x_{k+1}, x_k)$ is below a chosen tolerance, or when the objective decrease becomes negligible. The cyclic counterpart is analogous, cycling through the data points in a fixed order. ::: {.callout-tip title="Step-Size Conditions"} The condition $\sum \lambda_k = \infty$ ensures the algorithm can move arbitrarily far (sufficient exploration), while $\sum \lambda_k^2 < \infty$ ensures the steps eventually become small enough to settle (sufficient damping). A typical choice is $\lambda_k = c/k$ for some $c > 0$. ::: ## Convergence The following result specializes Corollary 3.9 of @Bacak2014 to the sample Fréchet objective. ::: {#prp-ppa-frechet-convergence .proposition title="Convergence of PPA for Fréchet means (Bačák 2014)"} Let $(\mathcal{M}, d)$ be a locally compact Hadamard space, let $X_1, \ldots, X_n \in \mathcal{M}$, and let $w_i = 1/n$. Consider the weighted sample Fréchet function $$ f(x) = \sum_{i=1}^n w_i\,d^2(x, X_i). $$ Let $(\lambda_k)$ be a sequence of positive numbers such that $$ \sum_{k=0}^{\infty} \lambda_k = \infty, \qquad \sum_{k=0}^{\infty} \lambda_k^2 < \infty. $$ Starting from $x_0 \in \mathcal{M}$, the sequences defined by the cyclic and randomized split PPA converge to the sample Fréchet mean as $k \to \infty$ (for the randomized scheme, almost surely). ::: The proximal-point algorithm is particularly attractive in Hadamard spaces: convexity guarantees a unique barycenter, while geodesic interpolation makes every split proximal step **explicit**. ## Interactive Exploration: Proximal Algorithm on the Sphere The demo below runs the randomized split PPA for computing the Fréchet mean of points on a **quadrant of the unit sphere** $\mathbb{S}^2$. **Data generation.** $n$ points $\{X_1, \ldots, X_n\}$ are sampled on the sphere octant $\{x, y, z \ge 0\}$ clustered around a randomly chosen center $c \in \mathbb{S}^2$: $$X_i = \gamma_{c \leftarrow U_i}(W_i), \qquad U_i \sim \text{Uniform}(\text{octant}), \quad W_i \sim \text{Uniform}(0.4, 1.0),$$ where $\gamma_{c \leftarrow u}(w)$ denotes the point at fraction $w$ along the great-circle geodesic from $u$ to $c$. This produces a random cluster in the octant whose Fréchet mean is near $c$. **Algorithm.** At each iteration $k$, the PPA picks a random data point $X_{r_k}$ and moves the current iterate $x_k$ along the geodesic toward it: $$x_{k+1} = \gamma_{x_k \to X_{r_k}}(t_k), \qquad t_k = \frac{2\lambda_k / n}{1 + 2\lambda_k / n}, \qquad \lambda_k = \frac{c}{k+1}.$$ On the sphere, the geodesic from $p$ to $q$ with angular distance $\theta = \arccos\langle p, q\rangle$ is $$\gamma_{p \to q}(t) = \frac{\sin((1-t)\theta)}{\sin\theta}\,p \;+\; \frac{\sin(t\theta)}{\sin\theta}\,q,$$ which stays on the unit sphere for all $t \in [0,1]$. **Drag the iteration slider to step through each update** and watch how the formulas translate into geometry. ```{ojs} //| label: fig-ppa-sphere //| fig-cap: "Interactive: PPA on a sphere quadrant with step-by-step visualization" // ---- Controls (2-column grid layout) ---- n_pts_control = Inputs.range([3, 25], {step: 1, value: 10, label: "Number of data points"}) c_val_control = Inputs.range([0.1, 5], {step: 0.1, value: 1.0, label: "Step-size constant c"}) data_seed_control = Inputs.range([1, 100], {step: 1, value: 42, label: "Random seed"}) step_control = Inputs.range([0, MAX_ITER], {step: 1, value: 0, label: "Iteration step"}) n_pts = Generators.input(n_pts_control) c_val = Generators.input(c_val_control) data_seed = Generators.input(data_seed_control) step = Generators.input(step_control) // ---- 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 helpers ---- 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))); } function frechetObj(x, dataPts) { var n = dataPts.length; var sum = 0; for (var i = 0; i < n; i++) { var d = sphereDist(x, dataPts[i]); sum += d * d; } return sum / n; } // ---- Orthographic projection ---- 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]; } function generateSphereData(n, seed) { var rng = mulberry32(seed); var center = randomOctantPoint(rng); var pts = []; for (var i = 0; i < n; i++) { var randPt = randomOctantPoint(rng); var w = 0.4 + rng() * 0.6; pts.push(sphereGeodesic(randPt, center, w)); } return pts; } // ---- PPA computation ---- MAX_ITER = 200; function runSpherePPA(dataPts, cVal, maxIter, baseSeed) { var n = dataPts.length; var weight = 1/n; var rng = mulberry32(baseSeed + 9999); var x = [1/Math.sqrt(3), 1/Math.sqrt(3), 1/Math.sqrt(3)]; var trajectory = [{pos: x.slice(), obj: frechetObj(x, dataPts), stepTo: -1}]; for (var k = 0; k < maxIter; k++) { var lambda_k = cVal / (k + 1); var idx = Math.floor(rng() * n); var a = dataPts[idx]; var t_k = 2 * lambda_k * weight / (1 + 2 * lambda_k * weight); x = sphereGeodesic(x, a, t_k); trajectory.push({pos: x.slice(), obj: frechetObj(x, dataPts), stepTo: idx}); } return trajectory; } function referenceMean(dataPts, baseSeed) { var n = dataPts.length; var weight = 1/n; var rng = mulberry32(baseSeed + 77777); var x = [1/Math.sqrt(3), 1/Math.sqrt(3), 1/Math.sqrt(3)]; for (var k = 0; k < 3000; k++) { var lambda_k = 0.3 / (1 + k/100); var idx = Math.floor(rng() * n); var a = dataPts[idx]; var t_k = 2 * lambda_k * weight / (1 + 2 * lambda_k * weight); x = sphereGeodesic(x, a, t_k); } return x; } data_pts = generateSphereData(n_pts, data_seed); trajectory = runSpherePPA(data_pts, c_val, MAX_ITER, data_seed); ref_mean = referenceMean(data_pts, data_seed); cur = trajectory[step]; // ---- Two-column slider panel ---- html` <style> .ppa-slider-grid { display: grid; grid-template-columns: repeat(2, minmax(0, 1fr)); gap: 6px 20px; max-width: 560px; margin: 0 auto 10px auto; font-family: system-ui, sans-serif; font-size: 0.85em; } .ppa-slider-grid > * { min-width: 0; margin: 0; } @media (max-width: 500px) { .ppa-slider-grid { grid-template-columns: 1fr; } } </style> <div class="ppa-slider-grid"> <div>${n_pts_control}</div> <div>${c_val_control}</div> <div>${data_seed_control}</div> <div>${step_control}</div> </div> ` // ---- SVG visualization ---- function buildSVG(dataPts, trajectory, stepIdx, refMean) { var curState = trajectory[stepIdx]; var curPos = curState.pos; 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.25, 0.5, 0.75].forEach(function(zv) { gridLines.push(constantZArc(zv, 30)); }); var dataProj = dataPts.map(function(p) { return project(p); }); var trajProj = []; for (var i = 0; i <= stepIdx; i++) { trajProj.push(project(trajectory[i].pos)); } var curProj = project(curPos); var refProj = project(refMean); 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.6"/>'); }); // 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"/>'); }); // Trajectory path if (trajProj.length > 1) { var trajPts = trajProj.map(function(p) { return p[0].toFixed(1) + "," + p[1].toFixed(1); }).join(" "); parts.push('<polyline points="' + trajPts + '" fill="none" stroke="#37474f" stroke-width="2.2" stroke-linecap="round" stroke-linejoin="round" opacity="0.75"/>'); } // Geodesic toward next data point if (curState.stepTo >= 0 && stepIdx < trajectory.length - 1) { var targetPt = dataPts[curState.stepTo]; var targetProj = project(targetPt); parts.push('<line x1="' + curProj[0].toFixed(1) + '" y1="' + curProj[1].toFixed(1) + '" x2="' + targetProj[0].toFixed(1) + '" y2="' + targetProj[1].toFixed(1) + '" stroke="#ff9800" stroke-width="1.5" stroke-dasharray="5,3" opacity="0.7"/>'); parts.push('<circle cx="' + targetProj[0].toFixed(1) + '" cy="' + targetProj[1].toFixed(1) + '" r="6" fill="none" stroke="#ff9800" stroke-width="2.2"/>'); } // Data points dataProj.forEach(function(dp, idx) { var isTarget = curState.stepTo === idx; var r = isTarget ? 4.5 : 3.3; var fill = isTarget ? "#ff9800" : "#e53935"; var stroke = isTarget ? "#e65100" : "#b71c1c"; parts.push('<circle cx="' + dp[0].toFixed(1) + '" cy="' + dp[1].toFixed(1) + '" r="' + r + '" fill="' + fill + '" stroke="' + stroke + '" stroke-width="0.8"/>'); parts.push('<title>Data point ' + (idx+1) + '</title>'); }); // Current iterate parts.push('<circle cx="' + curProj[0].toFixed(1) + '" cy="' + curProj[1].toFixed(1) + '" r="8.5" fill="#43a047" stroke="#fff" stroke-width="2.5"/>'); parts.push('<circle cx="' + curProj[0].toFixed(1) + '" cy="' + curProj[1].toFixed(1) + '" r="8.5" fill="none" stroke="#43a047" stroke-width="1" opacity="0.35">' + '<animate attributeName="r" values="8.5;15;8.5" dur="2s" repeatCount="indefinite"/></circle>'); // Reference Fréchet mean parts.push('<polygon points="' + (refProj[0]).toFixed(1) + ',' + (refProj[1]-7).toFixed(1) + ' ' + (refProj[0]+5.5).toFixed(1) + ',' + refProj[1].toFixed(1) + ' ' + refProj[0].toFixed(1) + ',' + (refProj[1]+7).toFixed(1) + ' ' + (refProj[0]-5.5).toFixed(1) + ',' + refProj[1].toFixed(1) + '" fill="#f9a825" stroke="#f57f17" stroke-width="1.3"/>'); // 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 = 52; parts.push('<circle cx="' + lx + '" cy="' + ly0 + '" r="4" fill="#e53935"/>'); parts.push('<text x="' + (lx+10) + '" y="' + (ly0+4) + '" font-size="10" fill="#37474f">Data point</text>'); parts.push('<circle cx="' + lx + '" cy="' + (ly0+17) + '" r="6" fill="#43a047" stroke="#fff" stroke-width="1.2"/>'); parts.push('<text x="' + (lx+10) + '" y="' + (ly0+21) + '" font-size="10" fill="#37474f">Current iterate</text>'); parts.push('<polygon points="' + lx + ',' + (ly0+31) + ' ' + (lx+4.5) + ',' + (ly0+37) + ' ' + lx + ',' + (ly0+43) + ' ' + (lx-4.5) + ',' + (ly0+37) + '" fill="#f9a825"/>'); parts.push('<text x="' + (lx+10) + '" y="' + (ly0+41) + '" font-size="10" fill="#37474f">Ref. Fréchet mean</text>'); parts.push('<line x1="' + lx + '" y1="' + (ly0+57) + '" x2="' + (lx+16) + '" y2="' + (ly0+57) + '" stroke="#ff9800" stroke-width="1.2" stroke-dasharray="3,2"/>'); parts.push('<text x="' + (lx+22) + '" y="' + (ly0+61) + '" font-size="10" fill="#37474f">Geodesic step</text>'); return '<svg viewBox="0 0 500 420" xmlns="http://www.w3.org/2000/svg" style="width:100%;max-width:520px;">\n' + parts.join('\n') + '\n</svg>'; } html`<div style="text-align:center;">${buildSVG(data_pts, trajectory, step, ref_mean)}</div>` // ---- Formula panel ---- n_weight = 1 / data_pts.length; cw = cur.stepTo >= 0 ? data_pts[cur.stepTo] : null; lam_k = step > 0 ? c_val / step : null; t_k_val = lam_k !== null ? (2 * lam_k * n_weight) / (1 + 2 * lam_k * n_weight) : null; geo_theta = cw !== null ? sphereDist(cur.pos, cw) : null; geo_w0 = geo_theta !== null && geo_theta > 1e-12 ? Math.sin((1 - t_k_val) * geo_theta) / Math.sin(geo_theta) : null; geo_w1 = geo_theta !== null && geo_theta > 1e-12 ? Math.sin(t_k_val * geo_theta) / Math.sin(geo_theta) : null; next_from_geo = cw !== null && geo_w0 !== null ? [ geo_w0 * cur.pos[0] + geo_w1 * cw[0], geo_w0 * cur.pos[1] + geo_w1 * cw[1], geo_w0 * cur.pos[2] + geo_w1 * cw[2] ] : null; html` <div style="font-family:system-ui,sans-serif; max-width:520px; margin:0 auto; font-size:0.88em;"> <div style="display:flex; justify-content:space-between; align-items:center; background:#f5f5f5; border-radius:6px; padding:8px 14px; margin-bottom:10px;"> <span><b>Step</b> ${step} / ${MAX_ITER}</span> <span style="color:#e53935;">→ data point <b>${cur.stepTo >= 0 ? "#" + (cur.stepTo + 1) : "—"}</b></span> <span><b>F<sub>n</sub></b> = ${cur.obj.toFixed(5)}</span> </div> ${step === 0 ? html` <div style="background:#fff3e0; border-left:3px solid #ff9800; border-radius:4px; padding:10px 14px; margin-bottom:10px;"> <b>Initialization:</b> x<sub>0</sub> = (1/√3, 1/√3, 1/√3) — the center of the octant.<br> Drag the slider to step through the PPA iterations. </div>` : (lam_k !== null && cw !== null ? html` <div style="display:grid; grid-template-columns:1fr 1fr; gap:10px 16px;"> <div style="background:#fafafa; border-radius:6px; padding:10px 14px; border:1px solid #e0e0e0;"> <div style="font-weight:600; margin-bottom:6px; color:#1565c0;">⓵ Step size</div> <div style="margin-bottom:4px;"> λ<sub>${step}</sub> = c/k = ${c_val.toFixed(1)}/${step} = <b>${lam_k.toFixed(4)}</b> </div> <div> t<sub>${step}</sub> = <span style="font-size:0.82em;">2λw/(1+2λw)</span> = <span style="font-size:0.82em;">2·${lam_k.toFixed(4)}·${n_weight.toFixed(4)}/(1+2·${lam_k.toFixed(4)}·${n_weight.toFixed(4)})</span><br> = <b>${t_k_val.toFixed(5)}</b> </div> </div> <div style="background:#fafafa; border-radius:6px; padding:10px 14px; border:1px solid #e0e0e0;"> <div style="font-weight:600; margin-bottom:6px; color:#1565c0;">⓶ Geodesic on 𝕊²</div> <div style="margin-bottom:4px;"> θ = arccos⟨x<sub>${step-1}</sub>, X<sub>${cur.stepTo+1}</sub>⟩ = <b>${geo_theta.toFixed(4)}</b> rad </div> <div style="font-size:0.82em;"> w<sub>0</sub> = sin((1−t)θ)/sinθ = <b>${geo_w0.toFixed(5)}</b><br> w<sub>1</sub> = sin(tθ)/sinθ = <b>${geo_w1.toFixed(5)}</b> </div> </div> <div style="grid-column:1/-1; background:#e8f5e9; border-radius:6px; padding:10px 14px; border:1px solid #a5d6a7;"> <div style="font-weight:600; margin-bottom:4px; color:#2e7d32;">⓷ PPA Update</div> <div style="font-size:0.9em;"> x<sub>${step}</sub> = w<sub>0</sub>·x<sub>${step-1}</sub> + w<sub>1</sub>·X<sub>${cur.stepTo+1}</sub> </div> <div style="font-size:0.78em; color:#555; margin-top:2px;"> = ${geo_w0.toFixed(4)}·(${cur.pos[0].toFixed(3)}, ${cur.pos[1].toFixed(3)}, ${cur.pos[2].toFixed(3)})<br> + ${geo_w1.toFixed(4)}·(${cw[0].toFixed(3)}, ${cw[1].toFixed(3)}, ${cw[2].toFixed(3)}) </div> ${next_from_geo ? html` <div style="font-size:0.82em; margin-top:3px; font-weight:600; color:#2e7d32;"> = (${next_from_geo[0].toFixed(4)}, ${next_from_geo[1].toFixed(4)}, ${next_from_geo[2].toFixed(4)}) [✓ on 𝕊²] </div>` : ""} </div> <div style="grid-column:1/-1; background:#f5f5f5; border-radius:6px; padding:8px 14px;"> <span><b>‖x<sub>${step}</sub> − μ<sub>ref</sub>‖</b> = ${sphereDist(cur.pos, ref_mean).toFixed(5)} rad</span> <span style="margin-left:16px;"><b>F<sub>n</sub>(x<sub>${step}</sub>)</b> = ${cur.obj.toFixed(5)}</span> </div> </div>` : "") } </div> ` // ---- Convergence plot ---- function buildConvergencePlot(trajectory, stepIdx, maxIter) { var W = 500, H = 180; var margin = {top: 15, right: 15, bottom: 28, left: 48}; var pw = W - margin.left - margin.right; var ph = H - margin.top - margin.bottom; var objs = trajectory.map(function(t) { return t.obj; }); var yMin = Math.min.apply(null, objs) * 0.95; var yMax = Math.max.apply(null, objs) * 1.02; var yRange = yMax - yMin || 0.001; function xScale(i) { return margin.left + (i / maxIter) * pw; } function yScale(v) { return margin.top + ph - ((v - yMin) / yRange) * 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 <= 3; gy++) { var yv = yMin + (gy/3) * yRange; var yp = yScale(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="9" fill="#9e9e9e">' + yv.toFixed(4) + '</text>'); } var linePts = ""; for (var i = 0; i < trajectory.length; i++) { linePts += xScale(i).toFixed(1) + "," + yScale(trajectory[i].obj).toFixed(1) + " "; } parts.push('<polyline points="' + linePts + '" fill="none" stroke="#1976d2" stroke-width="1.8" stroke-linejoin="round"/>'); var cx = xScale(stepIdx); var cy = yScale(trajectory[stepIdx].obj); parts.push('<line x1="' + cx + '" y1="' + margin.top + '" x2="' + cx + '" y2="' + (margin.top+ph) + '" stroke="#c92a2a" stroke-width="1.5" stroke-dasharray="4,3"/>'); parts.push('<circle cx="' + cx + '" cy="' + cy + '" r="4.5" fill="#c92a2a" stroke="#fff" stroke-width="1.5"/>'); var sx0 = xScale(0), sy0 = yScale(trajectory[0].obj); parts.push('<circle cx="' + sx0 + '" cy="' + sy0 + '" r="2.5" fill="#78909c"/>'); var lastIdx = trajectory.length - 1; var sxf = xScale(lastIdx), syf = yScale(trajectory[lastIdx].obj); parts.push('<circle cx="' + sxf + '" cy="' + syf + '" r="2.5" fill="#2e7d32"/>'); 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="11" fill="#616161">Iteration</text>'); parts.push('<text x="' + (margin.left-38) + '" y="' + (margin.top+ph/2) + '" text-anchor="middle" font-size="11" fill="#616161" ' + 'transform="rotate(-90,' + (margin.left-38) + ',' + (margin.top+ph/2) + ')">F_n(x_k)</text>'); parts.push('<text x="' + cx + '" y="' + (margin.top-2) + '" text-anchor="middle" font-size="10" fill="#c92a2a" font-weight="bold">k=' + stepIdx + '</text>'); return '<svg viewBox="0 0 ' + W + ' ' + H + '" xmlns="http://www.w3.org/2000/svg" style="width:100%;max-width:520px;">\n' + parts.join('\n') + '\n</svg>'; } html`<div style="text-align:center; margin-top:8px;">${buildConvergencePlot(trajectory, step, MAX_ITER)}</div>` ``` ::: {.callout-tip title="Try these experiments"} - **Drag the iteration slider from 0 to 200** and watch the green iterate trace great-circle arcs across the octant toward the gold reference Fréchet mean. - **Small $c$ (e.g., 0.3):** conservative steps. The iterate creeps along; convergence is slow but monotonic. - **Large $c$ (e.g., 3--5):** aggressive early moves. The first few steps jump near individual data points, then the algorithm gradually settles. - **Change the random seed** to see different data configurations --- some are tightly clustered, others spread across the octant. - **Increase $n$** to make the Fréchet mean more sharply defined; the trajectory tightens around the reference. - **Watch the orange dashed geodesic**: it shows the great-circle arc to the data point chosen at the current step. The green iterate moves a fraction $t_k$ along it. - **Look at the convergence plot**: $F_n(x_k)$ decreases over time, with occasional small bumps when a distant data point is randomly selected. - **Why the sphere quadrant?** The octant $\{x, y, z \ge 0\}$ is geodesically convex --- geodesics between points in the octant stay within it, so the PPA never leaves the region. ::: ## Key Takeaways - In **Hadamard spaces**, the empirical Fréchet function is convex, enabling proximal algorithms. - The **proximal map** $J_{\lambda}^f(x)$ generalizes the Euclidean proximal operator and is uniquely defined in Hadamard spaces. - **Split proximal-point algorithms** (cyclic or randomized) update one data point at a time, avoiding the full minimization. - Each component resolvent has a **closed form**: geodesic interpolation toward the data point with step size $t = \frac{2\lambda w}{1 + 2\lambda w}$. - The step-size sequence must satisfy $\sum \lambda_k = \infty$ (sufficient exploration) and $\sum \lambda_k^2 < \infty$ (damping). - Convergence is guaranteed in locally compact Hadamard spaces [@Bacak2014]. ## Exercises 1. **Component resolvent derivation**: Derive the closed form $t = \frac{2\lambda w}{1 + 2\lambda w}$ for the component resolvent on a geodesic. <a href="javascript:void(0)" onclick="showSolution('l4-sol-1')" class="solution-link">📝 Show Solution</a> 2. **Step-size conditions**: Show that $\lambda_k = c/k$ satisfies $\sum \lambda_k = \infty$ and $\sum \lambda_k^2 < \infty$. What happens if $\lambda_k = c/k^{1/2}$ or $\lambda_k = c/k^2$? <a href="javascript:void(0)" onclick="showSolution('l4-sol-2')" class="solution-link">📝 Show Solution</a> 3. **Cyclic vs. randomized**: Implement both the cyclic and randomized split PPA for points in $\mathbb{R}^2$ (with the Euclidean metric). Compare their convergence empirically. Does one consistently outperform the other? <a href="javascript:void(0)" onclick="showSolution('l4-sol-3')" class="solution-link">📝 Show Solution</a> 4. **SPD matrices**: Describe how you would implement geodesic interpolation on the SPD matrix space under the affine-invariant metric. What is the formula for the geodesic between $\Sigma_0$ and $\Sigma_1$? <a href="javascript:void(0)" onclick="showSolution('l4-sol-4')" class="solution-link">📝 Show Solution</a> 5. **Convergence rate**: For the simple case of computing the Euclidean mean of $n$ points, what is the convergence rate of the randomized PPA? How does it compare to the standard stochastic gradient descent? <a href="javascript:void(0)" onclick="showSolution('l4-sol-5')" 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="l4-sol-1"> <div class="solution-dialog"> <div class="solution-dialog-header"> <h4>Exercise 1: Component Resolvent Derivation</h4> <button onclick="closeSolution('l4-sol-1')" class="solution-dialog-close">✕ Close</button> </div> <div class="solution-original"> <strong>Exercise:</strong> Derive the closed form $t = \frac{2\lambda w}{1 + 2\lambda w}$ for the component resolvent on a geodesic. </div> <div class="solution-answer"> <strong>Solution:</strong> We wish to minimize $$\phi(z) = w\,d^2(z, a) + \frac{1}{2\lambda}d^2(z, y)$$ over $z \in \mathcal{M}$, where $a$ is a data point and $y$ is the current iterate. In a Hadamard space, the minimizer lies on the unique geodesic connecting $y$ and $a$. Let $\gamma(t) = (1-t)y + t a$ for $t \in [0,1]$, where "$+$" denotes geodesic interpolation. For a constant-speed geodesic, $d(\gamma(t), y) = t\,d(y, a)$ and $d(\gamma(t), a) = (1-t)\,d(y, a)$. Let $D = d(y, a)$. Then the objective restricted to the geodesic becomes: \begin{align*} \phi(\gamma(t)) &= w(1-t)^2 D^2 + \frac{1}{2\lambda}t^2 D^2 \\ &= D^2\left[w(1-t)^2 + \frac{1}{2\lambda}t^2\right]. \end{align*} This is a quadratic function of $t$. To find the minimizer, differentiate with respect to $t$ and set to zero: \begin{align*} \frac{d}{dt}\phi(\gamma(t)) &= D^2\left[2w(1-t)(-1) + \frac{1}{\lambda}t\right] \\ &= D^2\left[-2w(1-t) + \frac{t}{\lambda}\right] = 0. \end{align*} Assuming $D > 0$ (if $D = 0$, then $y = a$ and any $t$ works), we solve: \begin{align*} -2w(1-t) + \frac{t}{\lambda} &= 0 \\ \frac{t}{\lambda} &= 2w(1-t) \\ t &= 2\lambda w - 2\lambda w t \\ t(1 + 2\lambda w) &= 2\lambda w \\ t &= \frac{2\lambda w}{1 + 2\lambda w}. \end{align*} Since $\lambda > 0$ and $w > 0$, we have $t \in (0,1)$, confirming that the minimizer lies strictly between $y$ and $a$. **Behavior:** - As $\lambda \to 0$: $t \approx 2\lambda w$ (small, conservative step). - As $\lambda \to \infty$: $t \to 1$ (move all the way to the data point $a$). </div> </div> </dialog> <dialog id="l4-sol-2"> <div class="solution-dialog"> <div class="solution-dialog-header"> <h4>Exercise 2: Step-Size Conditions</h4> <button onclick="closeSolution('l4-sol-2')" class="solution-dialog-close">✕ Close</button> </div> <div class="solution-original"> <strong>Exercise:</strong> Show that $\lambda_k = c/k$ satisfies $\sum \lambda_k = \infty$ and $\sum \lambda_k^2 < \infty$. What happens if $\lambda_k = c/k^{1/2}$ or $\lambda_k = c/k^2$? </div> <div class="solution-answer"> <strong>Solution:</strong> **Case $\lambda_k = c/k$:** $$\sum_{k=1}^\infty \lambda_k = c\sum_{k=1}^\infty \frac{1}{k} = \infty,$$ since the harmonic series diverges. This ensures the algorithm can travel arbitrarily far — the total available step length is unbounded. $$\sum_{k=1}^\infty \lambda_k^2 = c^2\sum_{k=1}^\infty \frac{1}{k^2} = c^2\frac{\pi^2}{6} < \infty,$$ since the series $\sum 1/k^2$ (a $p$-series with $p=2>1$) converges. This ensures the steps eventually become small enough for the algorithm to settle. ✅ Both conditions satisfied. **Case $\lambda_k = c/k^{1/2}$:** $$\sum_{k=1}^\infty \lambda_k = c\sum_{k=1}^\infty \frac{1}{k^{1/2}} = \infty \quad \text{(diverges, } p=1/2 \le 1\text{)}.$$ $$\sum_{k=1}^\infty \lambda_k^2 = c^2\sum_{k=1}^\infty \frac{1}{k} = \infty \quad \text{(harmonic series)}.$$ ❌ Fails the second condition — the steps do not decay fast enough for the algorithm to stabilize. The algorithm may oscillate indefinitely. **Case $\lambda_k = c/k^2$:** $$\sum_{k=1}^\infty \lambda_k = c\sum_{k=1}^\infty \frac{1}{k^2} = c\frac{\pi^2}{6} < \infty.$$ $$\sum_{k=1}^\infty \lambda_k^2 = c^2\sum_{k=1}^\infty \frac{1}{k^4} < \infty.$$ ❌ Fails the first condition — the total step length is bounded, so if the initial point is far from the true minimizer, the algorithm may never reach it. **Summary:** Only the borderline case $\lambda_k = c/k$ (or more generally $\lambda_k \asymp 1/k$) satisfies both requirements simultaneously. </div> </div> </dialog> <dialog id="l4-sol-3"> <div class="solution-dialog"> <div class="solution-dialog-header"> <h4>Exercise 3: Cyclic vs. Randomized PPA</h4> <button onclick="closeSolution('l4-sol-3')" class="solution-dialog-close">✕ Close</button> </div> <div class="solution-original"> <strong>Exercise:</strong> Implement both the cyclic and randomized split PPA for points in $\mathbb{R}^2$. Compare their convergence empirically. </div> <div class="solution-answer"> <strong>Solution:</strong> **Implementation notes.** In $\mathbb{R}^2$ with the Euclidean metric, the geodesic is the straight line segment, and the component resolvent update for data point $X_i$ with weight $w_i = 1/n$ is: $$x_{k+1} = (1 - t_k)x_k + t_k X_i, \qquad t_k = \frac{2\lambda_k/n}{1 + 2\lambda_k/n}.$$ For the **cyclic** variant, iterate through data points $1, 2, \ldots, n, 1, 2, \ldots$ in order. For the **randomized** variant, sample $i$ uniformly from $\{1, \ldots, n\}$ at each step. **Empirical comparison.** Key observations from numerical experiments: 1. **Randomized typically converges faster** in terms of wall-clock time per effective epoch, because it avoids systematic biases that can arise from ordered data. If the data happen to be sorted (e.g., by some coordinate), the cyclic variant can spend many iterations "drifting" in one direction before cycling back. 2. **Cyclic can be faster** when data are well-shuffled and the objective is well-conditioned, since it ensures every data point is visited exactly once per cycle. This avoids the variance introduced by random sampling. 3. **Pathological cases.** If the data are sorted spatially, cyclic PPA can exhibit slow convergence because it chases points in a systematic pattern. The randomized variant is robust to data ordering. 4. **Theoretical guarantee.** The randomized PPA converges almost surely under the standard step-size conditions. The cyclic PPA converges as well, but theoretical guarantees are slightly more subtle (requiring the cyclic order to be "sufficiently mixing"). **Practical recommendation:** Shuffle the data once and use cyclic PPA for reproducibility, or use randomized PPA for robustness. In practice, the difference is often small compared to the effect of choosing good step sizes. </div> </div> </dialog> <dialog id="l4-sol-4"> <div class="solution-dialog"> <div class="solution-dialog-header"> <h4>Exercise 4: SPD Geodesic Interpolation</h4> <button onclick="closeSolution('l4-sol-4')" class="solution-dialog-close">✕ Close</button> </div> <div class="solution-original"> <strong>Exercise:</strong> Describe how you would implement geodesic interpolation on the SPD matrix space under the affine-invariant metric. What is the formula for the geodesic between $\Sigma_0$ and $\Sigma_1$? </div> <div class="solution-answer"> <strong>Solution:</strong> **The affine-invariant geodesic.** Under the affine-invariant Riemannian metric $$g_{\Sigma}(U, V) = \operatorname{tr}(\Sigma^{-1}U\,\Sigma^{-1}V),$$ the unique geodesic from $\Sigma_0$ to $\Sigma_1$ in $\mathcal{S}_{++}^{m}$ is given by: $$\boxed{\Sigma(t) = \Sigma_0^{1/2}\bigl(\Sigma_0^{-1/2}\,\Sigma_1\,\Sigma_0^{-1/2}\bigr)^t\Sigma_0^{1/2}, \qquad t \in [0,1].}$$ Here $\Sigma_0^{1/2}$ denotes the symmetric matrix square root, and $M^t = \exp(t\log M)$ denotes the matrix power (computed via the spectral decomposition: if $M = U\Lambda U^\top$ with $\Lambda = \operatorname{diag}(\lambda_i)$, then $M^t = U\operatorname{diag}(\lambda_i^t)U^\top$). **Computational steps:** 1. Compute the matrix square root $\Sigma_0^{1/2}$ and its inverse $\Sigma_0^{-1/2}$ via eigendecomposition: $$\Sigma_0 = V D V^\top \;\Longrightarrow\; \Sigma_0^{1/2} = V D^{1/2} V^\top, \quad \Sigma_0^{-1/2} = V D^{-1/2} V^\top.$$ 2. Form the "whitened" matrix $M = \Sigma_0^{-1/2}\,\Sigma_1\,\Sigma_0^{-1/2}$, which is also SPD. 3. Compute $M^t = \exp(t \log M)$ via eigendecomposition $M = U \Lambda U^\top$: $$M^t = U \operatorname{diag}(\lambda_1^t, \ldots, \lambda_m^t) U^\top.$$ 4. Assemble: $\Sigma(t) = \Sigma_0^{1/2}\,M^t\,\Sigma_0^{1/2}$. **PPA update.** For the Fréchet mean computation under the affine-invariant metric, the component PPA step becomes: $$t_k = \frac{2\lambda_k w_i}{1 + 2\lambda_k w_i}, \qquad \Sigma_{k+1} = \Sigma_k^{1/2}\bigl(\Sigma_k^{-1/2}\,\Sigma_i\,\Sigma_k^{-1/2}\bigr)^{t_k}\Sigma_k^{1/2}.$$ Each step requires two matrix square roots and one matrix power — computationally $O(m^3)$ per iteration, dominated by the eigendecompositions. **Alternative (Log-Euclidean).** Under the log-Euclidean metric, the geodesic is simply $\Sigma(t) = \exp((1-t)\log\Sigma_0 + t\log\Sigma_1)$, and the PPA update reduces to vector operations in the log-domain — much cheaper computationally. </div> </div> </dialog> <dialog id="l4-sol-5"> <div class="solution-dialog"> <div class="solution-dialog-header"> <h4>Exercise 5: Convergence Rate of Randomized PPA</h4> <button onclick="closeSolution('l4-sol-5')" class="solution-dialog-close">✕ Close</button> </div> <div class="solution-original"> <strong>Exercise:</strong> For the simple case of computing the Euclidean mean of $n$ points, what is the convergence rate of the randomized PPA? How does it compare to standard SGD? </div> <div class="solution-answer"> <strong>Solution:</strong> **Setup.** For $n$ points $X_1, \ldots, X_n \in \mathbb{R}^p$, the Fréchet objective is $$f(x) = \frac{1}{n}\sum_{i=1}^n \|x - X_i\|^2.$$ The unique minimizer is the sample mean $\bar{X} = \frac{1}{n}\sum_{i=1}^n X_i$. **Randomized PPA.** With weight $w_i = 1/n$ and step size $\lambda_k = c/k$, the update for a randomly chosen data point $X_{r_k}$ is: $$x_{k+1} = (1 - t_k)x_k + t_k X_{r_k}, \qquad t_k = \frac{2c/(nk)}{1 + 2c/(nk)} \approx \frac{2c}{nk}.$$ This is an **implicit** (proximal) stochastic method. Using the step-size analysis for stochastic approximation, one can show: $$\mathbb{E}\|x_k - \bar{X}\|^2 = O\left(\frac{1}{k^{2c/n}}\right)$$ for the expected squared distance, when $2c/n < 1$. Optimizing the constant $c$ yields roughly $O(1/k)$ behavior, matching the rate of SGD. **Comparison with SGD.** Standard SGD for the same problem with learning rate $\eta_k$: $$x_{k+1} = x_k - \eta_k \nabla f_{r_k}(x_k) = x_k - \eta_k \cdot 2(x_k - X_{r_k}) = (1 - 2\eta_k)x_k + 2\eta_k X_{r_k}.$$ This has the same structural form as the PPA update, with the identification $t_k \leftrightarrow 2\eta_k$. For $\eta_k \propto 1/k$, SGD achieves $\mathbb{E}\|x_k - \bar{X}\|^2 = O(1/k)$ as well. **Key difference.** The PPA is **implicit** (uses the gradient at the *next* iterate), which gives it better stability properties, especially for non-smooth or constrained objectives. For smooth, strongly convex problems like the Euclidean mean, both methods achieve the same asymptotic rate. However: - **PPA**: more stable with larger step sizes; step size $t_k \in (0,1)$ is automatically bounded. - **SGD**: can overshoot if $\eta_k$ is too large; requires step-size tuning. - **Both**: achieve $O(1/k)$ rate in squared distance, equivalent to $O(1/\sqrt{k})$ in distance to optimum. This is **not** the parametric $O(1/n)$ rate one gets from the exact sample mean — that requires using all data at each step (the full PPA or batch gradient descent). The split PPA trades per-iteration cost for slower convergence, making it suitable for large $n$. </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> ## Further Reading - @Bacak2014 — The definitive reference for proximal algorithms in Hadamard spaces. - @Sturm2003 — Foundational treatment of convexity and probability on $\mathrm{CAT}(0)$ spaces. - @PennecFillardAyache2006 — Practical algorithms for Fréchet means on SPD matrices. - @ArsignyEtAl2006; @ArsignyEtAl2007 — The log-Euclidean framework, enabling vector-space computation for SPD matrices. ## Self-Assessment Quiz Test your understanding of this lecture with the interactive MCQ quiz: 👉 **[Lecture 5 Quiz — 10 Multiple-Choice Questions](../quizzes/lecture-05-quiz.qmd)**