Lecture 9: Fréchet Regression — Total-Variation Regularized Fréchet Regression

Fused-lasso and change-point detection for metric-space trajectories

1 Learning Goals

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

• Identify scenarios where local smoothing is inappropriate and piecewise-constant structure is expected.
• Define total variation of a metric-space-valued curve using intrinsic metric length, and visualize it on the hemisphere.
• Formulate the TV-regularized Fréchet regression estimator and explain the role of the penalty parameter \(\lambda\).
• State the step-function representation theorem and its implications for computation.
• Describe the cyclic proximal-point algorithm (data-fitting and fusion steps) for computing the estimator.
• Interpret the \(n^{-1/3}\) convergence rate and explain why it is minimax optimal.

2 Motivation: When Smoothness Fails

The local Fréchet regression methods of Petersen and Müller (2019) are designed for smoothly varying regression functions. When the predictor is one-dimensional and ordered, however, the conditional Fréchet mean may exhibit long nearly constant stretches together with abrupt jumps. In such settings, local smoothing tends to blur discontinuities.

To address this, Lin and Müller (2021) proposed a total-variation regularized version of Fréchet regression that replaces local smoothing by a global variational fit with a metric-space total-variation penalty.

TipKey idea

Instead of fitting a smooth curve and hoping it captures jumps, directly fit a curve that is encouraged to be piecewise constant through a penalty on its total variation. Large \(\lambda\) values encourage fusion of adjacent fitted values and fewer jumps; small \(\lambda\) values track the data more closely.

3 Setup

The setup is specialized to a scalar ordered predictor, which we denote by deterministic design points

\[ a \le t_1 \le \cdots \le t_n \le b, \]

and responses \(Y_1, \ldots, Y_n \in \mathcal{M}\). The target mean curve \(\mu : [a, b] \to \mathcal{M}\) is defined at the design points by the conditional Fréchet means

\[ \mu(t_i) = \arg\min_{\omega \in \mathcal{M}} \mathbb{E}\bigl[d^2(\omega, Y_i)\bigr], \qquad i = 1, \ldots, n, \]

assuming these minimizers exist uniquely. In Hadamard spaces, uniqueness follows from the convexity of squared distance whenever a finite second moment exists (Sturm 2003; Afsari 2011).

4 Total Variation in Metric Spaces

For any metric-space-valued curve \(\gamma : [a, b] \to \mathcal{M}\), its total variation is defined intrinsically by the metric length

\[ \operatorname{TV}(\gamma) = \sup_{a = s_0 < \cdots < s_m = b} \sum_{j=1}^m d\bigl(\gamma(s_j), \gamma(s_{j-1})\bigr). \]

This is exactly the natural generalization of Euclidean total variation obtained by replacing absolute increments by metric distances.

To build intuition for total variation in metric spaces, the following interactive demo visualizes the TV of curves on the 2D sphere \(\mathbb{S}^2\), restricted to the upper hemisphere (\(z \ge 0\)).

The definition \(\operatorname{TV}(\gamma) = \sup_{a = s_0 < \cdots < s_m = b} \sum_{j=1}^m d(\gamma(s_j), \gamma(s_{j-1}))\) says: take a curve \(\gamma\), pick a partition of its domain, sum the geodesic distances between consecutive curve points, and take the supremum over all partitions. The geodesic segments between partition points are great-circle shortcuts — they represent the “as-the-crow-flies” distance on the sphere. As the partition is refined, these shortcuts hug the curve more closely, and their total length increases toward the true total variation (the Riemannian arc length of the curve).

Three curves are provided to explore:

• Great circle arc (geodesic): The curve itself is a geodesic — all partition points lie on the same great circle, so every partition gives the same sum, which equals the curve’s total length. No amount of refinement changes the TV.
• Sinusoidal wavy path: The curve meanders as it descends from the north pole to the equator. Coarse partitions produce geodesic segments that cut across the wiggles, underestimating the TV. As the partition is refined, the segments follow the curve more faithfully and the TV approaches its true value — exactly the supremum in action.
• Spiral: The curve spirals outward, creating a long path. A single geodesic segment (m = 1) is a direct shortcut; refining the partition reveals the spiral’s true length.

Visual guide:

• Blue upper hemisphere / gray lower hemisphere = the full 3D sphere, with the curve restricted to \(z \ge 0\)
• Dark blue curve = the curve \(\gamma(t)\)
• Red dots = partition points \(\gamma(s_j)\)
• Orange dashed arcs = geodesic segments \(d(\gamma(s_j), \gamma(s_{j-1}))\) (great-circle shortcuts)
• Gray grid lines = meridians and latitude circles for spatial reference; drag the sphere to change the viewing direction

tvDemoView = {
  var curveControl = Inputs.select(
    ["Great circle arc (geodesic)", "Sinusoidal wavy path", "Spiral"],
    {value: "Sinusoidal wavy path", label: "Curve type"}
  );
  var levelControl = Inputs.range(
    [0, 6],
    {step: 1, value: 2, label: "Refinement level r (segments m = 2ʳ)"}
  );

  var root = document.createElement("div");
  root.className = "tv-sphere-demo";
  root.innerHTML = `
    <style>
      .tv-sphere-demo { width:100%; font-family:system-ui,sans-serif; }
      .tv-sphere-controls {
        display:grid; grid-template-columns:minmax(220px,1fr) minmax(280px,2fr);
        gap:10px 24px; align-items:end; margin-bottom:10px;
      }
      .tv-sphere-controls > * { min-width:0; margin:0; }
      .tv-sphere-stage { width:80%; max-width:700px; margin:0 auto; }
      .tv-sphere-stage svg {
        display:block; width:100%; height:auto; cursor:grab; touch-action:none;
        user-select:none; -webkit-user-select:none;
      }
      .tv-sphere-hint { margin:-2px 0 8px; text-align:center; color:#546e7a; font-size:.82em; }
      .tv-sphere-stats {
        margin-top:8px; padding:12px 16px; background:#f8f9fa;
        border:1px solid #e9ecef; border-radius:7px; font-size:.9em;
      }
      .tv-sphere-stats table { width:100%; border-collapse:collapse; }
      .tv-sphere-stats td { padding:4px 8px; }
      .tv-sphere-stats td:last-child { text-align:right; font-family:monospace; }
      @media (max-width:700px) {
        .tv-sphere-controls { grid-template-columns:1fr; }
        .tv-sphere-stage { width:100%; }
      }
    </style>
    <div class="tv-sphere-controls"></div>
    <div class="tv-sphere-stage">
      <div class="tv-sphere-hint">Drag horizontally and vertically to rotate the sphere</div>
      <svg viewBox="0 0 700 560" role="img" aria-label="Rotatable three-dimensional sphere showing a curve, partition points, and geodesic shortcuts"></svg>
      <div class="tv-sphere-stats"></div>
    </div>`;

  var controls = root.querySelector(".tv-sphere-controls");
  controls.appendChild(curveControl);
  controls.appendChild(levelControl);
  var svg = root.querySelector("svg");
  var stats = root.querySelector(".tv-sphere-stats");
  var state = {rotY: -24, rotX: 56, dragging: false, lastX: 0, lastY: 0};
  var currentData = null;

  function clamp(v, lo, hi) {
    return Math.max(lo, Math.min(hi, v));
  }

  function sphereDist(p, q) {
    return Math.acos(clamp(p[0]*q[0] + p[1]*q[1] + p[2]*q[2], -1, 1));
  }

  function sphereGeodesic(p, q, t) {
    var dot = clamp(p[0]*q[0] + p[1]*q[1] + p[2]*q[2], -1, 1);
    var theta = Math.acos(dot);
    if (theta < 1e-10) return p.slice();
    var sinTheta = Math.sin(theta);
    var a = Math.sin((1-t)*theta) / sinTheta;
    var b = Math.sin(t*theta) / sinTheta;
    return [a*p[0] + b*q[0], a*p[1] + b*q[1], a*p[2] + b*q[2]];
  }

  function curvePoint(type, t) {
    var theta, phi;
    if (type === "Great circle arc (geodesic)") {
      return sphereGeodesic([0, 0, 1], [1, 0, 0], t);
    }
    theta = (Math.PI/2) * t;
    phi = type === "Sinusoidal wavy path"
      ? Math.PI*t + 1.2*Math.sin(2.5*Math.PI*t)
      : 4*Math.PI*t;
    return [
      Math.sin(theta)*Math.cos(phi),
      Math.sin(theta)*Math.sin(phi),
      Math.cos(theta)
    ];
  }

  function sampleCurve(type, n) {
    var points = [];
    for (var i = 0; i <= n; i++) points.push(curvePoint(type, i/n));
    return points;
  }

  function computeData(type, level) {
    var m = Math.pow(2, level);
    var partition = sampleCurve(type, m);
    var curve = sampleCurve(type, 320);
    var reference = sampleCurve(type, 4096);
    var shortcuts = [];
    var tvApprox = 0;
    var tvRef = 0;
    var i, j, segment;

    for (i = 1; i < partition.length; i++) {
      tvApprox += sphereDist(partition[i-1], partition[i]);
      segment = [];
      for (j = 0; j <= 12; j++) {
        segment.push(sphereGeodesic(partition[i-1], partition[i], j/12));
      }
      shortcuts.push(segment);
    }
    for (i = 1; i < reference.length; i++) {
      tvRef += sphereDist(reference[i-1], reference[i]);
    }
    return {
      type: type, level: level, m: m, partition: partition, curve: curve,
      shortcuts: shortcuts, tvApprox: tvApprox, tvRef: tvRef,
      isGeodesic: type === "Great circle arc (geodesic)"
    };
  }

  function rotatePoint(p) {
    var ry = state.rotY * Math.PI/180;
    var rx = state.rotX * Math.PI/180;
    var cy = Math.cos(ry), sy = Math.sin(ry);
    var cx = Math.cos(rx), sx = Math.sin(rx);
    var x1 = cy*p[0] + sy*p[2];
    var y1 = p[1];
    var z1 = -sy*p[0] + cy*p[2];
    return [x1, cx*y1 - sx*z1, sx*y1 + cx*z1];
  }

  function project(p) {
    var q = rotatePoint(p);
    return {x:350 + 238*q[0], y:278 - 238*q[1], z:q[2]};
  }

  function shadeColor(worldZ, depth) {
    var light = clamp(0.72 + 0.20*depth, 0.52, 0.92);
    var base = worldZ >= -1e-8 ? [188, 222, 242] : [220, 228, 234];
    return "rgb(" +
      Math.round(base[0]*light) + "," +
      Math.round(base[1]*light) + "," +
      Math.round(base[2]*light) + ")";
  }

  function addSegment(parts, a, b, style, bias) {
    var pa = project(a), pb = project(b);
    parts.push({
      depth:(pa.z + pb.z)/2 + bias,
      markup:'<line x1="' + pa.x.toFixed(1) + '" y1="' + pa.y.toFixed(1) +
        '" x2="' + pb.x.toFixed(1) + '" y2="' + pb.y.toFixed(1) + '" ' + style + '/>'
    });
  }

  function addPolyline(parts, points, style, bias) {
    for (var i = 1; i < points.length; i++) {
      addSegment(parts, points[i-1], points[i], style, bias);
    }
  }

  function projectedPolyline(points, style) {
    var projected = points.map(function(p) {
      var q = project(p);
      return q.x.toFixed(1) + "," + q.y.toFixed(1);
    }).join(" ");
    return '<polyline points="' + projected + '" fill="none" ' + style + '/>';
  }

  function renderScene() {
    var parts = [];
    var nTheta = 18, nPhi = 36;
    var ti, pi, theta0, theta1, phi0, phi1, p00, p01, p10, p11;

    // Opaque depth-sorted mesh: upper hemisphere is blue; lower hemisphere is gray.
    for (ti = 0; ti < nTheta; ti++) {
      theta0 = Math.PI*ti/nTheta;
      theta1 = Math.PI*(ti+1)/nTheta;
      for (pi = 0; pi < nPhi; pi++) {
        phi0 = 2*Math.PI*pi/nPhi;
        phi1 = 2*Math.PI*(pi+1)/nPhi;
        p00 = [Math.sin(theta0)*Math.cos(phi0), Math.sin(theta0)*Math.sin(phi0), Math.cos(theta0)];
        p01 = [Math.sin(theta0)*Math.cos(phi1), Math.sin(theta0)*Math.sin(phi1), Math.cos(theta0)];
        p10 = [Math.sin(theta1)*Math.cos(phi0), Math.sin(theta1)*Math.sin(phi0), Math.cos(theta1)];
        p11 = [Math.sin(theta1)*Math.cos(phi1), Math.sin(theta1)*Math.sin(phi1), Math.cos(theta1)];
        var q00 = project(p00), q01 = project(p01), q10 = project(p10), q11 = project(p11);
        var depth = (q00.z + q01.z + q10.z + q11.z)/4;
        var worldZ = (p00[2] + p01[2] + p10[2] + p11[2])/4;
        var fill = shadeColor(worldZ, depth);
        parts.push({
          depth:depth,
          markup:'<polygon points="' +
            q00.x.toFixed(1)+','+q00.y.toFixed(1)+' '+
            q01.x.toFixed(1)+','+q01.y.toFixed(1)+' '+
            q11.x.toFixed(1)+','+q11.y.toFixed(1)+' '+
            q10.x.toFixed(1)+','+q10.y.toFixed(1)+
            '" fill="'+fill+'" stroke="'+fill+'" stroke-width="1.2"/>'
        });
      }
    }

    // Full-sphere grid makes the rotation and occlusion visually explicit.
    var gridStyle = 'stroke="#607d8b" stroke-width="0.75" opacity="0.48"';
    for (pi = 0; pi < 12; pi++) {
      var meridian = [];
      phi0 = 2*Math.PI*pi/12;
      for (ti = 0; ti <= 72; ti++) {
        theta0 = Math.PI*ti/72;
        meridian.push([Math.sin(theta0)*Math.cos(phi0), Math.sin(theta0)*Math.sin(phi0), Math.cos(theta0)]);
      }
      addPolyline(parts, meridian, gridStyle, 0.003);
    }
    for (ti = 1; ti < 6; ti++) {
      var latitude = [];
      theta0 = Math.PI*ti/6;
      for (pi = 0; pi <= 96; pi++) {
        phi0 = 2*Math.PI*pi/96;
        latitude.push([Math.sin(theta0)*Math.cos(phi0), Math.sin(theta0)*Math.sin(phi0), Math.cos(theta0)]);
      }
      addPolyline(parts, latitude, gridStyle, 0.003);
    }

    var equator = [];
    for (pi = 0; pi <= 120; pi++) {
      phi0 = 2*Math.PI*pi/120;
      equator.push([Math.cos(phi0), Math.sin(phi0), 0]);
    }
    addPolyline(parts, equator, 'stroke="#455a64" stroke-width="1.8" opacity="0.8"', 0.004);

    parts.sort(function(a, b) { return a.depth - b.depth; });
    // Keep γ(t) as one foreground polyline. Sorting its individual segments with
    // mesh facets can create artificial gaps as the sphere rotates.
    var curveOverlay = projectedPolyline(
      currentData.curve,
      'stroke="#0756a3" stroke-width="4" stroke-linecap="round" stroke-linejoin="round"'
    );
    var shortcutOverlay = currentData.shortcuts.map(function(segment) {
      return projectedPolyline(
        segment,
        'stroke="#f08c00" stroke-width="2.8" stroke-dasharray="7,5" stroke-linecap="round" stroke-linejoin="round"'
      );
    }).join("");
    var pointOverlay = currentData.partition.map(function(p) {
      var q = project(p);
      return '<circle cx="'+q.x.toFixed(1)+'" cy="'+q.y.toFixed(1)+
        '" r="6" fill="#c62828" stroke="#fff" stroke-width="2"/>';
    }).join("");
    svg.innerHTML =
      '<rect width="700" height="560" fill="transparent"/>' +
      parts.map(function(d) { return d.markup; }).join("") +
      curveOverlay +
      shortcutOverlay +
      pointOverlay +
      '<g transform="translate(505,28)" font-family="system-ui,sans-serif" font-size="12" fill="#37474f">' +
        '<rect x="-12" y="-14" width="178" height="78" rx="5" fill="#fff" opacity=".88"/>' +
        '<line x1="0" y1="4" x2="24" y2="4" stroke="#0756a3" stroke-width="4"/>' +
        '<text x="32" y="8">Curve γ(t)</text>' +
        '<line x1="0" y1="28" x2="24" y2="28" stroke="#f08c00" stroke-width="2.5" stroke-dasharray="7,5"/>' +
        '<text x="32" y="32">Geodesic shortcut</text>' +
        '<circle cx="12" cy="52" r="5" fill="#c62828" stroke="#fff" stroke-width="1.5"/>' +
        '<text x="32" y="56">Partition point</text>' +
      '</g>';
  }

  function renderStats() {
    var gap = Math.max(0, currentData.tvRef - currentData.tvApprox);
    var gapPct = currentData.tvRef > 0 ? 100*gap/currentData.tvRef : 0;
    var note = currentData.isGeodesic
      ? '<div style="margin-top:8px;padding:7px 10px;background:#e8f5e9;border-left:3px solid #43a047;">' +
        '<b>Geodesic property:</b> every nested partition has the same length.</div>'
      : '<div style="margin-top:8px;padding:7px 10px;background:#e3f2fd;border-left:3px solid #1971c2;">' +
        '<b>Nested refinement:</b> moving the slider right retains every old partition point, so the approximation cannot decrease.</div>';
    stats.innerHTML =
      '<table><tr style="border-bottom:1px solid #dee2e6"><td colspan="2"><b>' +
      currentData.type + '</b> — level <b>r = ' + currentData.level +
      '</b>, <b>m = ' + currentData.m + '</b> segment' + (currentData.m === 1 ? '' : 's') +
      '</td></tr><tr><td>TV approximation</td><td><b>' +
      currentData.tvApprox.toFixed(4) + ' rad</b></td></tr>' +
      '<tr><td>Reference TV (4,096 segments)</td><td>' +
      currentData.tvRef.toFixed(4) + ' rad</td></tr>' +
      '<tr><td>Gap to reference</td><td>' + gapPct.toFixed(2) +
      '% (' + gap.toFixed(4) + ' rad)</td></tr>' +
      '<tr><td>Average segment length</td><td>' +
      (currentData.tvApprox/currentData.m).toFixed(4) + ' rad</td></tr></table>' + note;
  }

  function updateData() {
    currentData = computeData(curveControl.value, Number(levelControl.value));
    renderScene();
    renderStats();
  }

  curveControl.addEventListener("input", updateData);
  levelControl.addEventListener("input", updateData);
  svg.addEventListener("pointerdown", function(event) {
    state.dragging = true;
    state.lastX = event.clientX;
    state.lastY = event.clientY;
    svg.setPointerCapture(event.pointerId);
    svg.style.cursor = "grabbing";
  });
  svg.addEventListener("pointermove", function(event) {
    if (!state.dragging) return;
    var dx = event.clientX - state.lastX;
    var dy = event.clientY - state.lastY;
    state.lastX = event.clientX;
    state.lastY = event.clientY;
    state.rotY = (state.rotY + 0.55*dx) % 360;
    // Keep the sphere under the pointer: dragging down rotates it down.
    state.rotX = (state.rotX + 0.55*dy) % 360;
    renderScene();
  });
  function stopDrag(event) {
    if (!state.dragging) return;
    state.dragging = false;
    svg.style.cursor = "grab";
    if (svg.hasPointerCapture(event.pointerId)) svg.releasePointerCapture(event.pointerId);
  }
  svg.addEventListener("pointerup", stopDrag);
  svg.addEventListener("pointercancel", stopDrag);
  updateData();
  return root;
}

Figure 1: Interactive: Total variation of curves on the 2D sphere. Drag the sphere to rotate the 3D view.

TipKey observations

• Monotonicity: The demo uses nested dyadic partitions with \(m=2^r\) segments. The TV approximation \(\sum_j d(\gamma(s_j), \gamma(s_{j-1}))\) is therefore non-decreasing as \(r\) increases: every refinement retains the old points and adds new ones, so the triangle inequality can only increase (or leave unchanged) the sum.
• Geodesics: When the curve itself is a geodesic (great circle arc), the triangle inequality is an equality along the curve, so TV is constant regardless of the partition — a distinctive property.
• Convergence from below: The TV approximation always approaches the true total variation from below — the supremum is the least upper bound of all finite-partition sums.
• Metric generality: The same geometry plays out in any metric space: TV is the supremum of polygonal-path lengths, and geodesic segments are replaced by metric distances between consecutive curve points.

5 The TV-Regularized Fréchet Estimator

The total-variation regularized Fréchet estimator is any minimizer of

\[ L_{\lambda}(\gamma) = \frac{1}{n}\sum_{i=1}^n d^2\bigl(\gamma(t_i), Y_i\bigr) + \lambda \operatorname{TV}(\gamma), \]

that is,

\[ \hat{\mu} \in \arg\min_{\gamma : \operatorname{TV}(\gamma) < \infty} L_{\lambda}(\gamma). \]

The first term is the empirical Fréchet least-squares loss, while the penalty couples neighboring fitted values along the ordered predictor and discourages unnecessary oscillation. In contrast to kernel smoothing, this formulation does not impose local continuity/smoothness of \(\mu\) and is well suited to piecewise-constant or jump-like trajectories.

Step-Function Representation. Although the objective is posed over all curves of finite variation, the minimizer can always be chosen to be piecewise constant.

Proposition 1 Any minimizer \(\tilde{\mu}\) of \(L_{\lambda}\) can be replaced by a step function \(\hat{\mu}\) such that \(\hat{\mu}(t_i) = \tilde{\mu}(t_i)\) for all \(i = 1, \ldots, n\) and

\[ \operatorname{TV}(\hat{\mu}) \le \operatorname{TV}(\tilde{\mu}). \]

Consequently, one may take \(\hat{\mu}\) to be constant on each interval \([t_i, t_{i+1})\), with boundary values extended constantly on \([a, t_1)\) and \((t_n, b]\).

This proposition is one of the most important structural properties of the estimator. It shows that TV-regularized Fréchet regression is the metric analogue of one-dimensional total-variation denoising or fused-lasso regression: the fit is intrinsically piecewise constant, and the jump locations are selected adaptively from the data.

According to the step-function proposition, to compute the estimator one only needs the fitted values at the design points. Writing \(p_i = \gamma(t_i)\) reduces the infinite-dimensional problem to

\[ \widetilde{L}_{\lambda}(p_1, \ldots, p_n) = \frac{1}{2}\sum_{i=1}^n d^2(p_i, Y_i) + \frac{n\lambda}{2}\sum_{j=1}^{n-1} d(p_j, p_{j+1}), \]

an optimization problem on the product space \(\mathcal{M}^n\).

6 Computation: Cyclic Proximal-Point Algorithm

When \(\mathcal{M}\) is a Hadamard space, there is a unique geodesic between any two points. Let \(\gamma_{p,q} : [0, 1] \to \mathcal{M}\) denote the geodesic from \(p\) to \(q\). Lin and Müller (2021) adopt the cyclic proximal point scheme of Weinmann et al. (2014), alternating two kinds of updates:

Data-fitting step: Move each \(p_i\) toward \(Y_i\) along the geodesic,

\[ p_i \leftarrow \gamma_{p_i, Y_i}\!\left(\frac{\alpha_r}{1 + \alpha_r}\right). \]

Fusion step: Move adjacent fitted values toward one another,

\[ (p_j, p_{j+1}) \leftarrow \Bigl(\gamma_{p_j, p_{j+1}}(\theta_j), \gamma_{p_{j+1}, p_j}(\theta_j)\Bigr), \qquad \theta_j = \min\!\left\{\frac{\alpha_r n\lambda}{2 d(p_j, p_{j+1})}, \frac{1}{2}\right\}. \]

Here \(\alpha_r > 0\) is a step-size sequence satisfying

\[ \sum_{r=1}^{\infty} \alpha_r = \infty, \qquad \sum_{r=1}^{\infty} \alpha_r^2 < \infty. \]

The first update decreases the Fréchet least-squares loss, while the second update shrinks neighboring fitted values together and therefore reduces total variation. For Hadamard spaces the algorithm converges to a global minimizer of \(\widetilde{L}_{\lambda}\) (Weinmann et al. 2014).

NoteChoosing \(\lambda\) in practice

\(\lambda\) may be chosen by cross-validation, or by selecting the smallest value that yields a desired number of constant pieces. In the fMRI application below, the number of expected jumps is known from the experimental design, so \(\lambda\) is chosen to match that target.

7 Application: Dynamic Functional Connectivity from fMRI

An application in Lin and Müller (2021) concerns dynamic functional connectivity in the human brain. Functional connectivity is naturally summarized by covariance-type objects built from the joint evolution of BOLD (blood-oxygen-level dependent) signals across brain regions. These objects are symmetric positive-definite matrices and therefore live in \(\mathcal{S}_{++}^m\) rather than in an ordinary Euclidean vector space.

7.1 Scientific Question

The key question: do abrupt changes in external visual input induce corresponding abrupt changes in the mean connectivity pattern? If the brain remains in roughly stable connectivity regimes while a given clip is playing, and changes only when the stimulus changes, then the target mean curve should be approximately piecewise constant with a small number of jumps.

7.2 Data and Preprocessing

The data come from the Human Connectome Project social-cognition task. Out of 970 subjects, 850 had usable task-related fMRI data. During the experiment, each participant viewed five short video clips while in the scanner. The fMRI signal was recorded at 274 equally spaced time points. The five clips started and ended at approximately known times, giving 10 known stimulus change points.

For preprocessing:

• The brain was partitioned into 68 regions of interest; 8 regions related to social cognition were retained.
• For subject \(j\) at time \(i\), BOLD signals \(V_{ij} \in \mathbb{R}^8\) were recorded.
• Dynamic connectivity is represented by a moving-window covariance estimate with window \(P = 16\):

\[ \Sigma_{ij} = \frac{1}{P}\sum_{k=i-P}^{i+P-1} (V_{kj} - \bar{V}_{ij})(V_{kj} - \bar{V}_{ij})^\top. \]

• Subject-specific matrices \(\Sigma_{ij} \in \mathcal{S}_{++}^8\) are aggregated at each time point by taking the affine-invariant Fréchet mean:

\[ Y_i = \arg\min_{\Sigma \in \mathcal{S}_{++}^8} \frac{1}{850}\sum_{j=1}^{850} d_{\mathrm{AI}}^2(\Sigma, \Sigma_{ij}), \qquad i = 1, \ldots, 243. \]

7.3 TV-Regularized Fit

A direct specialization of the TV-regularized Fréchet estimator is applied to the SPD cone:

\[ \hat{\mu} \in \arg\min_{\gamma : \operatorname{TV}(\gamma) < \infty} \left[ \frac{1}{243}\sum_{i=1}^{243} d_{\mathrm{AI}}^2(\gamma(t_i), Y_i) + \lambda \operatorname{TV}(\gamma) \right]. \]

Total-variation regularized Fréchet regression for dynamic functional connectivity from task-related fMRI. Source: Lin and Müller (2021).

The fitted estimator with 9 jumps closely matches the known times at which video clips start and end. The solution path across \(\lambda\) values reveals:

• Weaker jumps fuse away first as \(\lambda\) increases.
• The most pronounced changes persist even under heavy regularization.
• Early stimulus changes produce the most persistent shifts in mean connectivity.

This is strong evidence that the method can recover meaningful experimental structure directly from a highly noisy sequence of SPD objects, without supplying the switching times to the estimation procedure.

8 Application: Detecting Market Regime Changes in Portfolio Covariance

The portfolio stress-testing application from Lecture 8 used local-linear Fréchet regression to model how the covariance matrix of asset returns varies smoothly with a market stress indicator. That approach assumes the conditional Fréchet mean \(\mu(x)\) is a smooth function of the stress level \(x\). In reality, financial markets often exhibit regime changes — abrupt transitions between distinct covariance structures triggered by events such as monetary policy shifts, financial crises, or sudden changes in volatility.

This is precisely the setting where TV-regularized Fréchet regression excels: the predictor (time or a market stress index) is ordered, and the response (a covariance matrix in \(\mathcal{S}_{++}^p\)) is expected to be approximately piecewise constant with jumps at regime boundaries.

8.1 From Smooth Variation to Regime Detection

In the Lecture 8 application, we modeled the log-covariance as a smooth function of a stress predictor \(X_t\) (realized SPY volatility). The local-linear estimator produced a curve that bends continuously with the predictor — appropriate when market conditions evolve gradually.

TV-regularized Fréchet regression answers a different question: are there distinct market regimes with sharply different covariance structures, and if so, where are the boundaries? Instead of smoothing through potential discontinuities, the TV penalty explicitly encourages a piecewise-constant fit, and the jump locations are estimated from the data.

8.2 Setup

Let \(t_1 < t_2 < \cdots < t_n\) be ordered observation times (e.g., trading days over a multi-year period). At each time \(t_i\), we observe a \(p \times p\) covariance matrix \(\Sigma_i \in \mathcal{S}_{++}^p\) computed from a rolling window of daily returns. The goal is to estimate a piecewise-constant mean curve \(\mu(t) \in \mathcal{S}_{++}^p\) that captures distinct covariance regimes.

Under the log-Euclidean metric \(d_{\mathrm{LE}}(\Sigma_1, \Sigma_2) = \|\log \Sigma_1 - \log \Sigma_2\|_F\), the TV-regularized estimator on the SPD cone becomes

\[ \hat{\mu} \in \arg\min_{\gamma : \operatorname{TV}(\gamma) < \infty} \left[ \frac{1}{n}\sum_{i=1}^n \|\log \gamma(t_i) - \log \Sigma_i\|_F^2 + \lambda \sum_{j=1}^{n-1} \|\log \gamma(t_{j+1}) - \log \gamma(t_j)\|_F \right]. \]

TipWhy the log-Euclidean metric simplifies TV regularization

Under the log-Euclidean metric, the SPD space is flat — geodesics are straight lines in the log-domain. The problem becomes one-dimensional group fused-lasso denoising of the matrix logarithms. Its fusion operation uses the vector soft-threshold

\[ \operatorname{prox}_{\tau\|\cdot\|_F}(V) = \left(1-\frac{\tau}{\|V\|_F}\right)_+V. \]

Thus the problem is a group fused lasso in the log-domain. The Frobenius norm couples all entries of each adjacent log-matrix difference, so covariance entries change at shared regime boundaries; it is not equivalent to fitting a separate scalar fused lasso to each entry.

8.3 Interactive Exploration: TV-Regularized Regime Detection on Simulated SPD Data

The following demo simulates \(2 \times 2\) SPD covariance matrices evolving through three distinct market regimes (low-volatility, crisis, recovery) with noisy observations. We fit the TV-regularized estimator and compare it to kernel smoothing (Lecture 7) and local-linear regression (Lecture 8).

Visual guide:

• Green dashed = true piecewise-constant mean (3 regimes with 2 jumps)
• Gray dots = noisy observed log-covariance entries
• Blue = TV-regularized fit (piecewise constant, this lecture)
• Orange = kernel-smoothed fit (Lecture 7 — blurs the jumps)
• Red = local-linear fit (Lecture 8 — also smooths through jumps)

tvPortfolioView = {
  var nControl = Inputs.range([40, 200], {step:10, value:100, label:"Number of time points n"});
  var lambdaControl = Inputs.range([0, 1.2], {step:0.02, value:0.46, label:"Regularization λ"});
  var noiseControl = Inputs.range([0.05, 0.5], {step:0.01, value:0.22, label:"Noise level σ"});
  var seedControl = Inputs.range([1, 100], {step:1, value:42, label:"Random seed"});

  var root = document.createElement("div");
  root.className = "tv-spd-demo";
  root.innerHTML = `
    <style>
      .tv-spd-demo { width:100%; font-family:system-ui,sans-serif; container-type:inline-size; }
      .tv-spd-controls {
        display:grid; grid-template-columns:repeat(2,minmax(0,1fr));
        gap:8px 24px; width:100%; margin:0 0 14px;
      }
      .tv-spd-controls > * { min-width:0; margin:0; }
      .tv-spd-controls input[type="number"] { width:7.5rem !important; }
      .tv-spd-output { width:80%; max-width:700px; margin:0 auto; }
      .tv-spd-plot + .tv-spd-plot { margin-top:12px; }
      .tv-spd-plot svg { display:block; width:100%; height:auto; }
      .tv-spd-diagnostics {
        margin-top:12px; padding:12px 16px; background:#f8f9fa;
        border:1px solid #e9ecef; border-radius:7px; font-size:.88em;
      }
      .tv-spd-diagnostics table { width:100%; border-collapse:collapse; }
      .tv-spd-diagnostics td { padding:4px 8px; vertical-align:top; }
      .tv-spd-diagnostics td:last-child { text-align:right; }
      @container (max-width:560px) {
        .tv-spd-controls { grid-template-columns:1fr; }
      }
      @media (max-width:700px) {
        .tv-spd-output { width:100%; }
      }
    </style>
    <div class="tv-spd-controls"></div>
    <div class="tv-spd-output">
      <div class="tv-spd-plot tv-spd-trace"></div>
      <div class="tv-spd-plot tv-spd-offdiag"></div>
      <div class="tv-spd-diagnostics"></div>
    </div>`;

  var controls = root.querySelector(".tv-spd-controls");
  [nControl, lambdaControl, noiseControl, seedControl].forEach(function(control) {
    controls.appendChild(control);
  });
  var traceHost = root.querySelector(".tv-spd-trace");
  var offHost = root.querySelector(".tv-spd-offdiag");
  var diagnostics = root.querySelector(".tv-spd-diagnostics");
  var bandwidth = 0.12;
  var pendingFrame = null;

  function mulberry32(seed) {
    return function() {
      seed |= 0; seed = seed + 0x6D2B79F5 | 0;
      var t = Math.imul(seed ^ seed >>> 15, 1 | seed);
      t = t + Math.imul(t ^ t >>> 7, 61 | t) ^ t;
      return ((t ^ t >>> 14) >>> 0) / 4294967296;
    };
  }

  function generateData(n, noise, seed) {
    var rng = mulberry32(seed);
    function randn() {
      var u1 = Math.max(rng(), 1e-15);
      return Math.sqrt(-2*Math.log(u1))*Math.cos(2*Math.PI*rng());
    }
    var regimes = [[0.5,0.1,0.6], [1.8,0.9,1.5], [0.8,0.3,0.7]];
    var jump1 = Math.floor(0.30*n), jump2 = Math.floor(0.65*n);
    var x = [], logs = [], truth = [];
    for (var i = 0; i < n; i++) {
      var mean = i < jump1 ? regimes[0] : (i < jump2 ? regimes[1] : regimes[2]);
      x.push(i/(n-1));
      truth.push(mean.slice());
      logs.push([
        mean[0] + noise*randn(),
        mean[1] + noise*randn()/Math.sqrt(2),
        mean[2] + noise*randn()
      ]);
    }
    return {x:x, logs:logs, truth:truth, jump1:jump1, jump2:jump2, n:n};
  }

  function toFrobenius(m) {
    return [m[0], Math.sqrt(2)*m[1], m[2]];
  }
  function fromFrobenius(v) {
    return [v[0], v[1]/Math.sqrt(2), v[2]];
  }
  function vectorNorm(v) {
    return Math.sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
  }

  function solveTridiagonal(rhs, rho) {
    var n = rhs.length;
    var cPrime = new Array(n-1), dPrime = new Array(n);
    var diag0 = 1 + rho;
    cPrime[0] = -rho/diag0;
    dPrime[0] = rhs[0]/diag0;
    for (var i = 1; i < n; i++) {
      var diag = 1 + (i === n-1 ? rho : 2*rho);
      var denom = diag + rho*cPrime[i-1];
      if (i < n-1) cPrime[i] = -rho/denom;
      dPrime[i] = (rhs[i] + rho*dPrime[i-1])/denom;
    }
    var solution = new Array(n);
    solution[n-1] = dPrime[n-1];
    for (var j = n-2; j >= 0; j--) {
      solution[j] = dPrime[j] - cPrime[j]*solution[j+1];
    }
    return solution;
  }

  // Exact group fused-lasso formulation in Frobenius coordinates, solved by ADMM.
  function fitGroupTV(logs, lambda) {
    var y = logs.map(toFrobenius);
    var n = y.length, dim = 3;
    if (lambda <= 1e-12) return {fit:logs.map(function(v){return v.slice();}), iterations:0};
    var alpha = 0.5*n*lambda;
    var rho = 1.5;
    var x = y.map(function(v){return v.slice();});
    var z = [], u = [];
    for (var i = 0; i < n-1; i++) {
      z.push([0,0,0]);
      u.push([0,0,0]);
    }
    var iterations = 0;
    for (var iter = 0; iter < 600; iter++) {
      var rhsByDim = [new Array(n), new Array(n), new Array(n)];
      for (var d = 0; d < dim; d++) {
        var w = new Array(n-1);
        for (i = 0; i < n-1; i++) w[i] = z[i][d] - u[i][d];
        rhsByDim[d][0] = y[0][d] - rho*w[0];
        for (i = 1; i < n-1; i++) rhsByDim[d][i] = y[i][d] + rho*(w[i-1] - w[i]);
        rhsByDim[d][n-1] = y[n-1][d] + rho*w[n-2];
      }
      var solved = [
        solveTridiagonal(rhsByDim[0], rho),
        solveTridiagonal(rhsByDim[1], rho),
        solveTridiagonal(rhsByDim[2], rho)
      ];
      for (i = 0; i < n; i++) x[i] = [solved[0][i], solved[1][i], solved[2][i]];

      var primalSq = 0, dualSq = 0;
      for (i = 0; i < n-1; i++) {
        var oldZ = z[i].slice();
        var q = [
          x[i+1][0] - x[i][0] + u[i][0],
          x[i+1][1] - x[i][1] + u[i][1],
          x[i+1][2] - x[i][2] + u[i][2]
        ];
        var normQ = vectorNorm(q);
        var shrink = normQ > 0 ? Math.max(0, 1 - alpha/(rho*normQ)) : 0;
        z[i] = [shrink*q[0], shrink*q[1], shrink*q[2]];
        for (d = 0; d < dim; d++) {
          var residual = x[i+1][d] - x[i][d] - z[i][d];
          u[i][d] += residual;
          primalSq += residual*residual;
          var zChange = z[i][d] - oldZ[d];
          dualSq += rho*rho*zChange*zChange;
        }
      }
      iterations = iter + 1;
      if (iter > 20 && Math.sqrt(primalSq) < 2e-5*Math.sqrt(3*(n-1)) &&
          Math.sqrt(dualSq) < 2e-5*Math.sqrt(3*(n-1))) break;
    }
    return {fit:x.map(fromFrobenius), iterations:iterations};
  }

  function gaussian(u) {
    return Math.exp(-0.5*u*u);
  }
  function weightedMean(logs, weights) {
    var out = [0,0,0];
    for (var i = 0; i < logs.length; i++) {
      out[0] += weights[i]*logs[i][0];
      out[1] += weights[i]*logs[i][1];
      out[2] += weights[i]*logs[i][2];
    }
    return out;
  }
  function fitKernel(logs, x, h) {
    return x.map(function(x0) {
      var weights = x.map(function(xj){return gaussian((xj-x0)/h);});
      var sum = weights.reduce(function(a,b){return a+b;},0);
      return weightedMean(logs, weights.map(function(w){return w/sum;}));
    });
  }
  function fitLocalLinear(logs, x, h) {
    var n = x.length;
    return x.map(function(x0) {
      var kernel = x.map(function(xj){return gaussian((xj-x0)/h);});
      var nu0=0, nu1=0, nu2=0;
      for (var j = 0; j < n; j++) {
        var dx = x[j]-x0;
        nu0 += kernel[j]/n;
        nu1 += kernel[j]*dx/n;
        nu2 += kernel[j]*dx*dx/n;
      }
      var denom = nu0*nu2 - nu1*nu1;
      var weights;
      if (Math.abs(denom) < 1e-12) {
        var sum = kernel.reduce(function(a,b){return a+b;},0);
        weights = kernel.map(function(w){return w/sum;});
      } else {
        weights = x.map(function(xj,j) {
          return kernel[j]*(nu2 - nu1*(xj-x0))/(n*denom);
        });
      }
      return weightedMean(logs, weights);
    });
  }

  function frobeniusDistance(a, b) {
    return vectorNorm([
      a[0]-b[0],
      Math.sqrt(2)*(a[1]-b[1]),
      a[2]-b[2]
    ]);
  }
  function mse(fit, truth) {
    var total = 0;
    for (var i = 0; i < fit.length; i++) {
      var distance = frobeniusDistance(fit[i], truth[i]);
      total += distance*distance;
    }
    return total/fit.length;
  }
  function detectJumps(fit) {
    var jumps = [];
    for (var i = 1; i < fit.length; i++) {
      if (frobeniusDistance(fit[i], fit[i-1]) > 0.06) jumps.push(i);
    }
    return jumps;
  }
  function meanLog(logs) {
    var weights = logs.map(function(){return 1/logs.length;});
    return weightedMean(logs, weights);
  }

  function linePath(values, xScale, yScale) {
    return values.map(function(v,i) {
      return (i === 0 ? "M" : "L") + xScale(i).toFixed(1) + "," + yScale(v).toFixed(1);
    }).join(" ");
  }
  function stepPath(values, xScale, yScale) {
    var path = "M" + xScale(0).toFixed(1) + "," + yScale(values[0]).toFixed(1);
    for (var i = 1; i < values.length; i++) {
      path += "L" + xScale(i).toFixed(1) + "," + yScale(values[i-1]).toFixed(1);
      path += "L" + xScale(i).toFixed(1) + "," + yScale(values[i]).toFixed(1);
    }
    return path;
  }

  function buildPlot(title, yLabel, values, data, jumps) {
    var W=700, H=285, margin={top:42,right:18,bottom:36,left:58};
    var pw=W-margin.left-margin.right, ph=H-margin.top-margin.bottom, n=data.n;
    var all = values.observed.concat(values.truth, values.tv, values.kernel, values.local);
    var yMin=Math.min.apply(null,all), yMax=Math.max.apply(null,all);
    var pad=Math.max(0.08,0.12*(yMax-yMin));
    yMin-=pad; yMax+=pad;
    function xS(i){return margin.left+i*pw/(n-1);}
    function yS(v){return margin.top+ph-(v-yMin)*ph/(yMax-yMin);}
    var p=[];
    p.push('<rect x="'+margin.left+'" y="'+margin.top+'" width="'+pw+'" height="'+ph+'" fill="#fafafa"/>');
    p.push('<rect x="'+margin.left+'" y="'+margin.top+'" width="'+(xS(data.jump1)-margin.left)+'" height="'+ph+'" fill="#e8f5e9" opacity=".38"/>');
    p.push('<rect x="'+xS(data.jump1)+'" y="'+margin.top+'" width="'+(xS(data.jump2)-xS(data.jump1))+'" height="'+ph+'" fill="#fce4ec" opacity=".32"/>');
    p.push('<rect x="'+xS(data.jump2)+'" y="'+margin.top+'" width="'+(margin.left+pw-xS(data.jump2))+'" height="'+ph+'" fill="#e3f2fd" opacity=".32"/>');
    for (var g=0; g<=4; g++) {
      var yv=yMin+g*(yMax-yMin)/4, yp=yS(yv);
      p.push('<line x1="'+margin.left+'" y1="'+yp+'" x2="'+(margin.left+pw)+'" y2="'+yp+'" stroke="#dfe3e6" stroke-width=".7"/>');
      p.push('<text x="'+(margin.left-6)+'" y="'+(yp+4)+'" text-anchor="end" font-size="10" fill="#6c757d">'+yv.toFixed(2)+'</text>');
    }
    values.observed.forEach(function(v,i){
      p.push('<circle cx="'+xS(i).toFixed(1)+'" cy="'+yS(v).toFixed(1)+'" r="2.1" fill="#868e96" opacity=".42"/>');
    });
    p.push('<path d="'+stepPath(values.truth,xS,yS)+'" fill="none" stroke="#2e7d32" stroke-width="2.5" stroke-dasharray="7,4"/>');
    p.push('<path d="'+linePath(values.kernel,xS,yS)+'" fill="none" stroke="#f08c00" stroke-width="2" opacity=".85"/>');
    p.push('<path d="'+linePath(values.local,xS,yS)+'" fill="none" stroke="#e53935" stroke-width="2" opacity=".78"/>');
    p.push('<path d="'+stepPath(values.tv,xS,yS)+'" fill="none" stroke="#1971c2" stroke-width="3.2" stroke-linejoin="round"/>');
    [data.jump1,data.jump2].forEach(function(j){
      p.push('<line x1="'+xS(j)+'" y1="'+margin.top+'" x2="'+xS(j)+'" y2="'+(margin.top+ph)+'" stroke="#2e7d32" stroke-width="1.2" stroke-dasharray="3,3" opacity=".7"/>');
    });
    jumps.forEach(function(j){
      p.push('<line x1="'+xS(j)+'" y1="'+margin.top+'" x2="'+xS(j)+'" y2="'+(margin.top+ph)+'" stroke="#1971c2" stroke-width="1.5" stroke-dasharray="5,3" opacity=".75"/>');
    });
    var grandY=yS(values.grand);
    p.push('<line x1="'+margin.left+'" y1="'+grandY+'" x2="'+(margin.left+pw)+'" y2="'+grandY+'" stroke="#795548" stroke-width="1.2" stroke-dasharray="5,4" opacity=".55"/>');
    p.push('<rect x="'+margin.left+'" y="'+margin.top+'" width="'+pw+'" height="'+ph+'" fill="none" stroke="#adb5bd"/>');
    p.push('<text x="'+(margin.left+pw/2)+'" y="17" text-anchor="middle" font-size="13" font-weight="600" fill="#343a40">'+title+'</text>');
    p.push('<text x="'+(margin.left+pw/2)+'" y="'+(H-5)+'" text-anchor="middle" font-size="11" fill="#495057">Ordered time</text>');
    p.push('<text x="14" y="'+(margin.top+ph/2)+'" text-anchor="middle" font-size="11" fill="#495057" transform="rotate(-90,14,'+(margin.top+ph/2)+')">'+yLabel+'</text>');
    var lx=margin.left+8, ly=31;
    var legend=[
      ["#2e7d32","True","7,4"],["#1971c2","TV fit",""],["#f08c00","Kernel",""],["#e53935","Local linear",""],["#795548","Grand mean","5,4"]
    ];
    legend.forEach(function(item,k){
      var x=lx+k*112;
      p.push('<line x1="'+x+'" y1="'+ly+'" x2="'+(x+19)+'" y2="'+ly+'" stroke="'+item[0]+'" stroke-width="2.2"'+(item[2]?' stroke-dasharray="'+item[2]+'"':'')+'/>');
      p.push('<text x="'+(x+24)+'" y="'+(ly+4)+'" font-size="10" fill="#495057">'+item[1]+'</text>');
    });
    return '<svg viewBox="0 0 '+W+' '+H+'" xmlns="http://www.w3.org/2000/svg">'+p.join("")+'</svg>';
  }

  function update() {
    var n=Number(nControl.value), lambda=Number(lambdaControl.value);
    var noise=Number(noiseControl.value), seed=Number(seedControl.value);
    var data=generateData(n,noise,seed);
    var tvResult=fitGroupTV(data.logs,lambda);
    var kernel=fitKernel(data.logs,data.x,bandwidth);
    var local=fitLocalLinear(data.logs,data.x,bandwidth);
    var grand=meanLog(data.logs);
    var jumps=detectJumps(tvResult.fit);
    function trace(m){return m[0]+m[2];}
    function offdiag(m){return m[1];}
    var traceValues={
      truth:data.truth.map(trace), observed:data.logs.map(trace), tv:tvResult.fit.map(trace),
      kernel:kernel.map(trace), local:local.map(trace), grand:trace(grand)
    };
    var offValues={
      truth:data.truth.map(offdiag), observed:data.logs.map(offdiag), tv:tvResult.fit.map(offdiag),
      kernel:kernel.map(offdiag), local:local.map(offdiag), grand:offdiag(grand)
    };
    traceHost.innerHTML=buildPlot("Trace of log Σ — regime levels","tr(log Σ)",traceValues,data,jumps);
    offHost.innerHTML=buildPlot("Off-diagonal of log Σ — correlation structure","(log Σ)₁₂",offValues,data,jumps);

    var tvMSE=mse(tvResult.fit,data.truth), kernelMSE=mse(kernel,data.truth), localMSE=mse(local,data.truth);
    var tolerance=Math.max(2,Math.round(0.03*n));
    var matched=jumps.length===2 &&
      Math.abs(jumps[0]-data.jump1)<=tolerance && Math.abs(jumps[1]-data.jump2)<=tolerance;
    var jumpMessage=matched
      ? '✓ Both regime boundaries recovered within '+tolerance+' time points.'
      : jumps.length>2
        ? 'Too many jumps: increase λ to fuse short, noisy regimes.'
        : jumps.length<2
          ? 'Too few jumps: decrease λ to retain weaker boundaries.'
          : 'Two jumps were found, but their locations do not match both true boundaries.';
    diagnostics.innerHTML=
      '<table><tr style="border-bottom:1px solid #dee2e6"><td colspan="2"><b>Regime-detection diagnostics</b> '+
      '(λ='+lambda.toFixed(2)+', h='+bandwidth.toFixed(2)+')</td></tr>'+
      '<tr><td>TV fit MSE (Frobenius²)</td><td style="color:#1971c2"><b>'+tvMSE.toFixed(4)+'</b></td></tr>'+
      '<tr><td>Kernel MSE (Frobenius²)</td><td style="color:#f08c00">'+kernelMSE.toFixed(4)+'</td></tr>'+
      '<tr><td>Local-linear MSE (Frobenius²)</td><td style="color:#e53935">'+localMSE.toFixed(4)+'</td></tr>'+
      '<tr><td>Detected jumps</td><td><b>'+jumps.length+'</b> (true: 2)</td></tr>'+
      '<tr><td>Detected indices</td><td style="font-family:monospace">'+(jumps.length?jumps.join(", "):"none")+'</td></tr>'+
      '<tr><td>ADMM iterations</td><td>'+tvResult.iterations+'</td></tr></table>'+
      '<div style="margin-top:8px;padding:7px 10px;background:'+(matched?'#e8f5e9':'#fff3e0')+
      ';border-left:3px solid '+(matched?'#43a047':'#f08c00')+'">'+jumpMessage+'</div>';
  }

  function scheduleUpdate() {
    if (pendingFrame !== null) cancelAnimationFrame(pendingFrame);
    pendingFrame=requestAnimationFrame(function(){pendingFrame=null;update();});
  }
  [nControl,lambdaControl,noiseControl,seedControl].forEach(function(control){
    control.addEventListener("input",scheduleUpdate);
  });
  update();
  return root;
}

Figure 2: Interactive: TV-regularized vs. kernel vs. local-linear Fréchet regression for market regime detection on 2×2 SPD matrices

TipTry these experiments

• Tune \(\lambda\) for jump recovery: Set \(\lambda \approx 0.4\)–\(0.6\) — the TV fit should recover approximately 2 jumps matching the true regime boundaries (30% and 65% of the time axis). At \(\lambda = 0\), every point is its own regime (\(n-1\) jumps). At \(\lambda \to \infty\), the fit collapses to the unconditional Fréchet mean (brown dashed line).
• Compare to smoothing methods: Kernel (orange) and local-linear (red) both produce smooth curves that blur across the true jump boundaries. This is the fundamental tradeoff: smoothing methods assume continuity, TV regularization allows and encourages sharp transitions.
• Increase noise (\(\sigma = 0.3\)–\(0.5\)): The advantage of TV regularization grows — at high noise, smoothing methods spread jump information across many neighboring points, while TV regularization concentrates it at the jump locations.
• Examine the jump locations: The blue vertical dashed lines show where the TV estimator places jumps. Compare these to the green dashed lines (true jumps). With well-chosen \(\lambda\), they should align closely.
• The regime background shading: Green = low-volatility regime, pink = crisis (high vol + high correlation), blue = recovery. The TV estimator’s constant pieces should match these levels.

8.4 When to Use TV Regularization vs. Smoothing for Portfolio Risk

The choice between the three methods depends on the nature of market dynamics:

Scenario	Recommended method	Reason
Gradual market evolution	Local-linear (L8)	Covariance changes continuously with economic conditions — smoothing is appropriate
Known structural breaks (e.g., policy announcements)	TV (L9) with \(\lambda\) tuned to expected jumps	TV preserves the sharp transitions that smoothing would blur
Exploratory regime detection	TV (L9) across a range of \(\lambda\)	The solution path reveals which changes are most persistent — analogous to the dendrogram in hierarchical clustering

9 Asymptotic Theory

The large-sample theory is formulated for mean curves of bounded variation, with the empirical metric

\[ d_n(\hat{\mu}, \mu) = \Biggl\{\frac{1}{n}\sum_{i=1}^n d^2(\hat{\mu}(t_i), \mu(t_i))\Biggr\}^{1/2} \]

quantifying the estimation error.

NoteAssumptions (H0)–(H2)

(H0) Bounded variation of the target: \(\mu \in \mathscr{G}_{\mathcal{M}}(C)\) for some fixed \(C > 0\), where \(\mathscr{G}_{\mathcal{M}}(C)\) is the class of curves \(\gamma : [a, b] \to \mathcal{M}\) with \(\operatorname{TV}(\gamma) \le C\).

(H1) Sub-Gaussian tails around the mean curve: There exist constants \(\beta, \zeta > 0\) such that

\[ \sup_{1 \le i \le n} \mathbb{E}\exp\!\bigl\{\beta d^2(\mu(t_i), Y_i)\bigr\} \le \zeta. \]

(H2) Local entropy control: For a fixed \(R > 0\) there exists \(K > 0\) such that

\[ \log N\bigl(\delta, \mathscr{G}_{\mathcal{M}}^{r}(r), d_n\bigr) \le K\delta^{-1} \]

for all \(\delta > 0\), all \(n \ge 1\), and all \(0 < r \le R\).

Theorem 1 Suppose the response space belongs to a family of Hadamard spaces for which the local entropy condition (H2) holds, the target curve satisfies \(\mu \in \mathscr{G}_{\mathcal{M}}(C)\), and (H1) is satisfied. If the regularization parameter is chosen as

\[ \lambda \asymp n^{-2/3}, \]

then

\[ d_n(\hat{\mu}, \mu) = O_p(n^{-1/3}) \]

uniformly over the model class.

Several consequences are worth emphasizing:

1. \(n^{-1/3}\) is the canonical rate for estimating a one-dimensional bounded-variation signal — the same rate as in classical TV denoising. The theorem shows this rate remains achievable even when the response lives in a nonlinear Hadamard space.

2. General entropy bound: If the entropy bound is replaced by \(\log N(\delta, \cdot) \le K\delta^{-\alpha}\) for \(\alpha \in (0, 2)\), the rate becomes \(O_p(n^{-1/(2+\alpha)})\).

ImportantSpaces covered

The entropy condition (H2) holds for the SPD matrix space (with affine-invariant or log-Euclidean metric), the BHV tree space, and the Wasserstein space \((\mathcal{P}(\mathbb{R}), W_1)\).

10 Key Takeaways

• TV-regularized Fréchet regression is designed for ordered predictors where the regression function is approximately piecewise constant with jumps — it replaces local smoothing with a global variational penalty.
• Total variation of a metric-space curve is defined intrinsically as the supremum of sum of distances along partitions — the natural generalization of \(\int |\gamma'|\).
• The step-function representation theorem shows that the infinite-dimensional problem reduces to optimization over \(n\) points in \(\mathcal{M}\), and the solution is piecewise constant with adaptively chosen jump locations.
• Computation alternates between data-fitting steps (moving toward observations along geodesics) and fusion steps (shrinking adjacent fitted values together).
• The convergence rate is \(O_p(n^{-1/3})\) with \(\lambda \asymp n^{-2/3}\) — the canonical rate for bounded-variation signal estimation, and minimax optimal.

11 Exercises

1. TV of a step function. Show that for a step function \(\gamma\) with \(J\) jumps at locations \(\tau_1, \ldots, \tau_J\) and values \(\omega_0, \omega_1, \ldots, \omega_J\), the total variation is \(\operatorname{TV}(\gamma) = \sum_{j=1}^J d(\omega_j, \omega_{j-1})\). 📝 Show Solution

2. Extreme \(\lambda\) behavior. What is \(\hat{\mu}\) when \(\lambda = 0\)? When \(\lambda \to \infty\)? Interpret both limits. 📝 Show Solution

3. Fusion threshold. In the fusion step, the parameter \(\theta_j = \min\{\alpha_r n\lambda / (2d(p_j, p_{j+1})), 1/2\}\). Explain why the \(1/2\) cap is necessary and what happens when \(\theta_j = 1/2\). 📝 Show Solution

4. Comparison with kernel smoothing. Contrast TV-regularized and kernel Fréchet regression for a regression function that is smooth everywhere except at two known jump points. Which method would you choose and why? 📝 Show Solution

5. Sparsely observed functional data. Extend TV-regularized Fréchet regression to sparsely observed functional data where the functions take values in a general metric space. Outline how you would approach this extension. 📝 Show Solution

Exercise 1: TV of a Step Function

Exercise: Show that for a step function, the total variation equals the sum of jump magnitudes.

Solution:

For a step function \(\gamma\) with jumps at \(\tau_j\), any partition that includes the jump points as partition points will have sum of increments equal to \(\sum d(\omega_j, \omega_{j-1})\). Refining the partition by adding points within constant segments adds zero increment (since \(\gamma\) is constant there). Therefore the supremum over all partitions is exactly \(\sum_j d(\omega_j, \omega_{j-1})\).

Formally: let the jump points be \(\tau_1 < \cdots < \tau_J\). For any partition \(\{s_k\}\), the sum of increments is at most the sum of jumps whose intervals \([\tau_{j-1}, \tau_j]\) are crossed, and the supremum is achieved by any partition that includes all \(\tau_j\).

Exercise 2: Extreme \(\lambda\)

Exercise: What is \(\hat{\mu}\) when \(\lambda = 0\)? When \(\lambda \to \infty\)?

Solution:

\(\lambda = 0\): No penalty on total variation. The objective is \(L_0(\gamma) = \frac{1}{n}\sum_i d^2(\gamma(t_i), Y_i)\). This separates across \(i\), and the minimizer satisfies \(\hat{\mu}(t_i) = Y_i\) for each \(i\). The fit is pure interpolation — every observation is its own Fréchet mean. Total variation is uncontrolled.

\(\lambda \to \infty\): The penalty dominates. Any variation is prohibitively expensive, so the minimizer is a constant function: \(\hat{\mu}(t_i) \equiv \hat{\omega}\) for all \(i\), where \(\hat{\omega}\) is the unconditional sample Fréchet mean of all \(Y_i\). This is the simplest possible model — a single constant.

Between the extremes, \(\lambda\) traces a solution path from interpolation (\(\lambda = 0\)) to the grand mean (\(\lambda \to \infty\)), with the number of constant pieces decreasing monotonically.

Exercise 3: Fusion Threshold

Exercise: Explain why the \(1/2\) cap is necessary in the fusion step.

Solution:

The fusion step moves \(p_j\) and \(p_{j+1}\) toward each other along the geodesic by a fraction \(\theta_j\). Without the \(1/2\) cap, when \(d(p_j, p_{j+1})\) is very small, \(\theta_j = \alpha_r n\lambda / (2 d(p_j, p_{j+1}))\) can exceed 1, meaning the points would overshoot each other — moving past the midpoint and potentially crossing.

The cap \(\theta_j \le 1/2\) ensures: 1. The points move toward each other but never cross. 2. When \(\theta_j = 1/2\), \(p_j\) and \(p_{j+1}\) are fused to the same point (the midpoint of the geodesic). 3. The proximal interpretation is preserved — \(\theta_j = 1/2\) corresponds to evaluating the proximal operator of the distance function. 4. The operator remains firmly nonexpansive, guaranteeing convergence in Hadamard spaces.

Geometric intuition: Moving each point halfway along the geodesic toward the other is exactly the operation that minimizes \(d(p_j, Y_j)^2 + d(p_{j+1}, Y_{j+1})^2 + \text{const} \cdot d(p_j, p_{j+1})\) — it’s the exact solution to the two-point fusion subproblem.

Exercise 4: TV vs. Kernel for Jump Discontinuities

Exercise: Contrast TV-regularized and kernel Fréchet regression for a function with two known jumps.

Solution:

TV-regularized regression wins when jumps are present. Here’s why:

• Kernel regression smooths across jumps — the estimated function shows gradual transitions rather than sharp jumps. The bandwidth \(h\) controls the smoothness; reducing \(h\) sharpens the transition but increases variance everywhere.

• TV-regularized regression preserves sharp jumps — the estimated function is exactly piecewise constant. Jumps are preserved as long as \(\lambda\) is not too large. The jump locations are detected adaptively.

Practical recommendation: - If the regression function is mostly smooth with a few jumps → use TV-regularization. The smooth parts will be well-approximated by steps if \(n\) is large enough. - If the regression function is mostly smooth without jumps → use local Fréchet regression (kernel or local-linear). - If you’re unsure → try both and compare. Cross-validation can select between them.

Exercise 5: Extension to Sparsely Observed Functional Data

Exercise: Outline an approach to extend TV-regularized Fréchet regression to sparsely observed metric-space functional data.

Solution:

Setup: Each subject \(j\) is observed at a sparse, subject-specific set of time points \(t_{j1}, \ldots, t_{jm_j}\) with responses \(Y_{j\ell} \in \mathcal{M}\). The goal is to estimate a common mean function \(\mu(t)\).

Approach outline:

1. Pooling + interpolation: First, estimate individual trajectories using Fréchet interpolation (e.g., geodesic interpolation between observation times), then apply TV-regularized Fréchet regression to the pooled interpolated values.

2. Joint estimation: Alternatively, directly minimize

   \[ L_\lambda(\gamma) = \frac{1}{\sum_j m_j} \sum_{j=1}^N \sum_{\ell=1}^{m_j} d^2(\gamma(t_{j\ell}), Y_{j\ell}) + \lambda \operatorname{TV}(\gamma), \]

   where \(\gamma\) is evaluated at all unique observation times. This avoids pre-smoothing.

3. Computation: The proximal algorithm still applies as long as we can compute geodesics in \(\mathcal{M}\) and the data-fitting step uses only the observations available at each time point.

12 Self-Assessment Quiz

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

👉 Lecture 9 Quiz — 10 Multiple-Choice Questions

13 Further Reading

• Lin and Müller (2021) — The foundational paper on TV-regularized Fréchet regression.
• Weinmann et al. (2014) — Cyclic proximal point algorithms for TV-regularized problems on manifolds.
• Petersen and Müller (2019) — The Fréchet regression framework that TV-regularization builds upon.
• Sturm (2003); Afsari (2011) — Hadamard space theory: convexity of squared distance, uniqueness of Fréchet means.

References

Afsari, Bijan. 2011. “Riemannian \(L^p\) Center of Mass: Existence, Uniqueness, and Convexity.” Proceedings of the American Mathematical Society 139 (2): 655–73. https://doi.org/10.1090/S0002-9939-2010-10541-5.

Lin, Zhenhua, and Hans-Georg Müller. 2021. “Total Variation Regularized Fréchet Regression for Metric-Space Valued Data.” The Annals of Statistics 49 (6): 3510–33. https://doi.org/10.1214/21-AOS2095.

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

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.

Weinmann, Andreas, Laurent Demaret, and Martin Storath. 2014. “Total Variation Regularization for Manifold-Valued Data.” SIAM Journal on Imaging Sciences 7 (4): 2226–57. https://doi.org/10.1137/130951075.