L1 vs. L2 Regression

// VM is the shared utility library, loaded globally via _includes/head-scripts.html (js/**). This page uses VM.numerical.linearRegression (L2), VM.numerical.l1Regression (L1), VM.plotting.paddedRange, VM.ui.renderTable.
VM = window.VM

// severityMax: the slider's bound in each direction (range is
// [-severityMax, severityMax]). Used both for the slider itself and (in
// setup, below) to fix the plot's y-axis to the slider's full reach — not
// the current severity value — so the axis never rescales as the slider
// moves. Without that, an autoscaling axis would shrink the L2 line's
// drift back down to looking small at every step, which defeats the
// point of the demo.
severityMax = 15

// basePoints redraws whenever regenerateTrigger changes (initial load
// counts as one draw), seeded from the current timestamp so each draw is a
// fresh, non-repeating set of points and outlier locations. Points carry
// their pre-outlier (noisy, on-the-line) y-value; outlierSeverity is
// applied on top of this in `points`, below, so dragging the slider never
// needs to re-draw the noise.
//
// Kept outside the ojs-row div below: an output:false cell still renders
// an empty .cell div, and if it sat inside the row it would become the
// row's first child, stealing the `flex: 1 1 200px` fill rule meant for
// the slider.
basePoints = {
  regenerateTrigger
  // Seeded PRNG (mulberry32), seeded from the current time so each
  // "Regenerate" click produces a different draw.
  let state = Date.now()
  const rand = () => {
    state = (state + 0x6D2B79F5) | 0
    let t = Math.imul(state ^ (state >>> 15), 1 | state)
    t = (t + Math.imul(t ^ (t >>> 7), 61 | t)) ^ t
    return ((t ^ (t >>> 14)) >>> 0) / 4294967296
  }
  // Box-Muller transform: turns two uniform samples into one standard-normal sample.
  const gaussian = () => {
    const u1 = Math.max(rand(), 1e-12)
    const u2 = rand()
    return Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2)
  }

  const totalPoints = 24
  const outlierCount = 4
  // Slope and intercept are redrawn on every regenerate too (not just the
  // noise), so the underlying line itself is different each time, not just
  // the scatter around it. The sign is randomized too, so the line slants
  // down as often as it slants up.
  const slope = (rand() < 0.5 ? -1 : 1) * (0.5 + rand() * 2)
  const intercept = rand() * 5
  const noiseStd = 0.5

  const pts = []
  for (let i = 0; i < totalPoints; i++) {
    const x = i * 10 / (totalPoints - 1)
    const y = slope * x + intercept + noiseStd * gaussian()
    pts.push({x, y})
  }

  // Pick outlierCount distinct indices, away from the very ends so the
  // outliers sit inside the visible span of the line rather than right at
  // its edge.
  const outlierIndices = []
  while (outlierIndices.length < outlierCount) {
    const i = 2 + Math.floor(rand() * (totalPoints - 4))
    if (!outlierIndices.includes(i)) outlierIndices.push(i)
  }

  return {points: pts, outlierIndices, slope, intercept}
}

viewof outlierSeverity = {
  const input = Inputs.range([-severityMax, severityMax], {step: 0.5, value: 5, label: "Outlier severity"})
  input.classList.add("ojs-fill")
  return input
}

// regenerateTrigger: clicking "Regenerate points" bumps this, which is the
// only thing basePoints depends on below — that's what makes a fresh draw
// happen on click instead of only once at page load.
viewof regenerateTrigger = {
  const button = Inputs.button("Regenerate points", {value: 0, reduce: v => v + 1})
  button.classList.add("ojs-auto")
  return button
}

// points: applies the current outlier severity on top of basePoints. Only
// the outlier indices' y-values change with the slider — each keeps its
// own noisy y-value from basePoints and is shifted by exactly
// `outlierSeverity` on top of it (positive pushes it up, negative pushes
// it down), so severity 0 means "not actually an outlier yet" (identical
// to its original noisy position) and larger magnitude means "further
// from where it should be." Shifting each point's own noisy y — rather
// than recomputing it from the line equation — is what keeps the 4
// outliers from landing exactly on a second, perfectly straight line as
// the slider moves. Every other point stays fixed as the slider moves.
points = {
  const pts = []
  for (let i = 0; i < basePoints.points.length; i++) {
    const p = basePoints.points[i]
    if (basePoints.outlierIndices.includes(i)) {
      pts.push({x: p.x, y: p.y + outlierSeverity})
    } else {
      pts.push(p)
    }
  }
  return pts
}

// result: fits both an L2 (ordinary least-squares) and an L1 (least
// absolute deviations) line to the current points, plus each fit's total
// squared/absolute residual and how long it took to compute. A single fit
// is sub-millisecond and too noisy to time meaningfully, so each fit runs
// `reps` times and we report the average — this is also what makes the
// timing comparison fair between the two methods (L2 is a closed-form
// formula; L1 is iterative and does more work per call).
result = {
  const reps = 200

  const t0 = performance.now()
  let l2
  for (let i = 0; i < reps; i++) l2 = VM.numerical.linearRegression(points)
  const t1 = performance.now()
  let l1
  for (let i = 0; i < reps; i++) l1 = VM.numerical.l1Regression(points)
  const t2 = performance.now()

  const residualStats = fit => {
    let sumSquared = 0, sumAbsolute = 0
    for (const p of points) {
      const residual = p.y - (fit.slope * p.x + fit.intercept)
      sumSquared += residual * residual
      sumAbsolute += Math.abs(residual)
    }
    return {sumSquared, sumAbsolute}
  }

  return {
    points,
    l2: {...l2, ...residualStats(l2), timeMs: (t1 - t0) / reps},
    l1: {...l1, ...residualStats(l1), timeMs: (t2 - t1) / reps}
  }
}

// setup: axis ranges and line endpoints shared by mainPlot.
setup = {
  const xs = []
  for (const p of basePoints.points) xs.push(p.x)
  const xRange = VM.plotting.paddedRange(xs, {emptyRange: [-1, 11], relativePadding: 0.08, minPadding: 0.5})

  // Fixed y-range, sized from the slider's full [-severityMax, severityMax]
  // reach — not the current outlierSeverity value — so the axis never
  // rescales as the slider moves (see the note on severityMax above).
  const ys = []
  for (let i = 0; i < basePoints.points.length; i++) {
    const p = basePoints.points[i]
    ys.push(p.y)
    if (basePoints.outlierIndices.includes(i)) {
      ys.push(p.y + severityMax)
      ys.push(p.y - severityMax)
    }
  }
  const yRange = VM.plotting.paddedRange(ys, {emptyRange: [-2, 20], relativePadding: 0.12, minPadding: 1})

  const lineYs = fit => [fit.slope * xRange.lo + fit.intercept, fit.slope * xRange.hi + fit.intercept]

  return {
    points: result.points,
    outlierIndices: basePoints.outlierIndices,
    xRange, yRange,
    lineXs: [xRange.lo, xRange.hi],
    l2LineYs: lineYs(result.l2),
    l1LineYs: lineYs(result.l1)
  }
}

// mainPlot draws the inlier points, the outlier points, and both fitted
// lines.
//
// _plotDiv is kept in a closure so Plotly reuses the same DOM node across
// reactive updates — this preserves zoom and pan when the slider or button
// changes.
mainPlot = {
  const Plotly = window.Plotly
  let _plotDiv = null

  return ({points, outlierIndices, xRange, yRange, lineXs, l2LineYs, l1LineYs}) => {
    const inlierXs = [], inlierYs = [], outlierXs = [], outlierYs = []
    for (let i = 0; i < points.length; i++) {
      if (outlierIndices.includes(i)) {
        outlierXs.push(points[i].x)
        outlierYs.push(points[i].y)
      } else {
        inlierXs.push(points[i].x)
        inlierYs.push(points[i].y)
      }
    }

    const data = [
      {x: lineXs, y: l2LineYs, type: "scatter", mode: "lines", name: "L2 fit (least squares)", line: {color: "#dc2626", width: 2.5}},
      {x: lineXs, y: l1LineYs, type: "scatter", mode: "lines", name: "L1 fit (least absolute deviations)", line: {color: "#2563eb", width: 2.5}},
      {x: inlierXs, y: inlierYs, type: "scatter", mode: "markers", name: "Points", marker: {color: "#111827", size: 8, line: {color: "white", width: 1}}},
      {x: outlierXs, y: outlierYs, type: "scatter", mode: "markers", name: "Outliers", marker: {color: "#f97316", size: 13, symbol: "diamond", line: {color: "white", width: 1.5}}}
    ]

    const layout = {
      margin: {l: 0, r: 0, t: 0, b: 0},
      xaxis: {title: "x", range: [xRange.lo, xRange.hi], zeroline: true},
      yaxis: {title: "y", range: [yRange.lo, yRange.hi], zeroline: true},
      legend: {orientation: "h", y: -0.08},
      hovermode: "closest",
      uirevision: "static",
      autosize: true
    }
    const config = {responsive: true, displaylogo: false}

    if (!_plotDiv) {
      _plotDiv = document.createElement("div")
      _plotDiv.className = "plotly-box-large"
      Plotly.newPlot(_plotDiv, data, layout, config)
      // newPlot fires before the div is in the DOM, so Plotly measures 0 width.
      // A ResizeObserver catches the real size once it's inserted and laid out.
      new ResizeObserver(() => Plotly.Plots.resize(_plotDiv)).observe(_plotDiv)
    } else {
      Plotly.react(_plotDiv, data, layout, config)
    }
    return _plotDiv
  }
}

mainPlot(setup)

NoteFit comparison

// comparisonTable reports each fit's slope/intercept, its total squared
// and absolute residual over the current points, and the average time to
// compute it. L2 always has the smaller "sum of squared residuals" (it's
// exactly what it minimizes) and L1 always has the smaller "sum of
// absolute residuals" (same reason) — the column to watch as the slider
// moves is slope and intercept: L2's drift noticeably further from the
// underlying (pre-outlier) line — shown in the note below the table —
// than L1's does.
comparisonTable = {
  const row = (label, fit) => [
    label,
    fit.slope.toFixed(4),
    fit.intercept.toFixed(4),
    fit.sumSquared.toExponential(3),
    fit.sumAbsolute.toExponential(3),
    (fit.timeMs * 1000).toFixed(2)
  ]

  return VM.ui.renderTable({
    html,
    headers: ["Fit", "Slope", "Intercept", "Sum of squared residuals", "Sum of absolute residuals", "Avg. fit time (μs)"],
    csvHeaders: ["Fit", "Slope", "Intercept", "Sum of squared residuals", "Sum of absolute residuals", "Avg fit time (microseconds)"],
    filename: "l1-vs-l2-comparison.csv",
    rows: [row("L2 (least squares)", result.l2), row("L1 (least absolute deviations)", result.l1)]
  })
}

md`The underlying (pre-outlier) line has slope **${basePoints.slope.toFixed(2)}** and intercept **${basePoints.intercept.toFixed(2)}** — compare each fit's slope/intercept above to that as you move the slider.

L1's fit took **${result.l1.iterations}** reweighting iterations to converge at the current severity, which is why it's consistently slower than L2's single closed-form solve, even on a dataset this small.`

Why this matters

Drag the outlier severity slider (in either direction) to see the two fits diverge. L2 minimizes \(\sum_i (y_i - \hat y_i)^2\); L1 minimizes \(\sum_i |y_i - \hat y_i|\). Squaring a residual makes large residuals count far more than small ones, so L2 bends the line hard to shrink one huge outlier’s residual — that’s why it chases outliers. L1’s penalty grows only linearly with distance, so an outlier costs the fit proportionally and never dominates the way it does for L2, which is why L1 is called a robust estimator. The tradeoff is computational: L2 has a closed-form solution, while L1 has none and must be solved iteratively — this page uses iteratively reweighted least squares (IRLS), which is why L1’s fit takes more iterations and more time than L2’s, as shown in the table above.