PCA on a 2D Point Cloud

// PCA on a 2D point cloud. ~50 draggable points. Each frame we
// compute the sample mean, the 2x2 covariance matrix, and solve the
// eigenproblem analytically (closed-form for symmetric 2x2). Principal
// axes are drawn through the centroid, length scaled by sqrt(eigenvalue).
// An inset shows the 1D scatter of the PC1 projections.

const N = 50;
const HANDLE_R = 7;
const HIT_R = 18;

let xs, ys;
let dragging = -1;
let W = 0, H = 0;
let inited = false;

function seed(width, height) {
  W = width; H = height;
  xs = new Float32Array(N);
  ys = new Float32Array(N);
  // Anisotropic Gaussian-ish blob, tilted ~30deg, centred mid-screen.
  const cx = width * 0.5, cy = height * 0.55;
  const sx = Math.min(width, height) * 0.22;
  const sy = Math.min(width, height) * 0.07;
  const theta = -Math.PI / 7;
  const ct = Math.cos(theta), st = Math.sin(theta);
  for (let i = 0; i < N; i++) {
    // Box-Muller for unit normals.
    const u1 = Math.max(1e-6, Math.random());
    const u2 = Math.random();
    const z1 = Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);
    const z2 = Math.sqrt(-2 * Math.log(u1)) * Math.sin(2 * Math.PI * u2);
    const lx = z1 * sx;
    const ly = z2 * sy;
    xs[i] = cx + ct * lx - st * ly;
    ys[i] = cy + st * lx + ct * ly;
  }
  inited = true;
}

function init({ width, height }) {
  seed(width, height);
}

function eig2x2(a, b, c, d) {
  // Eigenvalues of [[a,b],[b,d]] (symmetric). Already real.
  const tr = a + d;
  const det = a * d - b * b;
  const disc = Math.max(0, tr * tr * 0.25 - det);
  const s = Math.sqrt(disc);
  const l1 = tr * 0.5 + s;
  const l2 = tr * 0.5 - s;
  // Eigenvector for l1: solve (A - l1 I) v = 0 -> [a-l1, b; b, d-l1] v = 0.
  let vx, vy;
  if (Math.abs(b) > 1e-9) {
    vx = l1 - d;
    vy = b;
  } else if (Math.abs(a - l1) > 1e-9 || Math.abs(d - l1) > 1e-9) {
    // diagonal: pick the axis with the larger value.
    if (a >= d) { vx = 1; vy = 0; }
    else { vx = 0; vy = 1; }
  } else {
    vx = 1; vy = 0;
  }
  const n = Math.hypot(vx, vy) || 1;
  vx /= n; vy /= n;
  return { l1, l2, v1x: vx, v1y: vy, v2x: -vy, v2y: vx };
}

function tick({ ctx, width, height, input }) {
  if (!inited || width !== W || height !== H) seed(width, height);
  const mx = input.mouseX, my = input.mouseY;

  if (input.mouseDown && dragging === -1) {
    let best = -1, bestD = HIT_R;
    for (let i = 0; i < N; i++) {
      const d = Math.hypot(mx - xs[i], my - ys[i]);
      if (d < bestD) { bestD = d; best = i; }
    }
    dragging = best;
  }
  if (!input.mouseDown) dragging = -1;
  if (dragging !== -1) {
    const m = 8;
    xs[dragging] = Math.max(m, Math.min(width - m, mx));
    ys[dragging] = Math.max(m, Math.min(height - m, my));
  }

  // Mean.
  let mxs = 0, mys = 0;
  for (let i = 0; i < N; i++) { mxs += xs[i]; mys += ys[i]; }
  mxs /= N; mys /= N;

  // Covariance.
  let sxx = 0, syy = 0, sxy = 0;
  for (let i = 0; i < N; i++) {
    const dx = xs[i] - mxs, dy = ys[i] - mys;
    sxx += dx * dx; syy += dy * dy; sxy += dx * dy;
  }
  sxx /= (N - 1); syy /= (N - 1); sxy /= (N - 1);

  const { l1, l2, v1x, v1y, v2x, v2y } = eig2x2(sxx, sxy, sxy, syy);
  const s1 = Math.sqrt(Math.max(0, l1));
  const s2 = Math.sqrt(Math.max(0, l2));

  // Background.
  ctx.fillStyle = "#0a0c14";
  ctx.fillRect(0, 0, width, height);

  // Grid.
  ctx.strokeStyle = "rgba(120,140,180,0.06)";
  ctx.lineWidth = 1;
  ctx.beginPath();
  for (let x = 0; x < width; x += 40) { ctx.moveTo(x, 0); ctx.lineTo(x, height); }
  for (let y = 0; y < height; y += 40) { ctx.moveTo(0, y); ctx.lineTo(width, y); }
  ctx.stroke();

  // 1-sigma ellipse (axes-aligned in the PC frame).
  ctx.save();
  ctx.translate(mxs, mys);
  ctx.rotate(Math.atan2(v1y, v1x));
  ctx.strokeStyle = "rgba(180,200,255,0.18)";
  ctx.lineWidth = 1;
  ctx.beginPath();
  ctx.ellipse(0, 0, s1, s2, 0, 0, Math.PI * 2);
  ctx.stroke();
  ctx.strokeStyle = "rgba(180,200,255,0.10)";
  ctx.beginPath();
  ctx.ellipse(0, 0, 2 * s1, 2 * s2, 0, 0, Math.PI * 2);
  ctx.stroke();
  ctx.restore();

  // Principal axes through centroid, length = 2.5 * sigma.
  const k = 2.5;
  ctx.lineWidth = 2.5;
  ctx.lineCap = "round";
  ctx.strokeStyle = "rgba(255,160,90,0.95)";
  ctx.beginPath();
  ctx.moveTo(mxs - v1x * s1 * k, mys - v1y * s1 * k);
  ctx.lineTo(mxs + v1x * s1 * k, mys + v1y * s1 * k);
  ctx.stroke();
  ctx.strokeStyle = "rgba(120,220,180,0.95)";
  ctx.beginPath();
  ctx.moveTo(mxs - v2x * s2 * k, mys - v2y * s2 * k);
  ctx.lineTo(mxs + v2x * s2 * k, mys + v2y * s2 * k);
  ctx.stroke();

  // Points.
  for (let i = 0; i < N; i++) {
    const isDrag = i === dragging;
    ctx.fillStyle = isDrag ? "rgba(255,230,140,0.98)" : "rgba(160,200,255,0.88)";
    ctx.strokeStyle = "rgba(255,255,255,0.6)";
    ctx.lineWidth = 1;
    ctx.beginPath();
    ctx.arc(xs[i], ys[i], isDrag ? HANDLE_R + 2 : HANDLE_R, 0, Math.PI * 2);
    ctx.fill();
    if (isDrag) ctx.stroke();
  }

  // Centroid marker.
  ctx.fillStyle = "rgba(255,255,255,0.95)";
  ctx.beginPath();
  ctx.arc(mxs, mys, 3.5, 0, Math.PI * 2);
  ctx.fill();

  // Inset: 1D scatter of PC1 projections.
  const insetW = Math.min(260, width - 24);
  const insetH = 56;
  const ix = 12;
  const iy = height - insetH - 12;
  ctx.fillStyle = "rgba(20,24,36,0.92)";
  ctx.fillRect(ix, iy, insetW, insetH);
  ctx.strokeStyle = "rgba(255,160,90,0.5)";
  ctx.lineWidth = 1;
  ctx.strokeRect(ix + 0.5, iy + 0.5, insetW - 1, insetH - 1);

  // Compute projections onto PC1.
  let pmin = Infinity, pmax = -Infinity;
  const proj = new Float32Array(N);
  for (let i = 0; i < N; i++) {
    const p = (xs[i] - mxs) * v1x + (ys[i] - mys) * v1y;
    proj[i] = p;
    if (p < pmin) pmin = p;
    if (p > pmax) pmax = p;
  }
  const span = Math.max(1e-3, pmax - pmin);
  const axisY = iy + insetH * 0.62;
  ctx.strokeStyle = "rgba(255,160,90,0.55)";
  ctx.beginPath();
  ctx.moveTo(ix + 10, axisY); ctx.lineTo(ix + insetW - 10, axisY);
  ctx.stroke();
  // Tick at centroid (projection = 0).
  const z0 = ix + 10 + ((0 - pmin) / span) * (insetW - 20);
  ctx.strokeStyle = "rgba(255,255,255,0.7)";
  ctx.beginPath();
  ctx.moveTo(z0, axisY - 6); ctx.lineTo(z0, axisY + 6); ctx.stroke();
  ctx.fillStyle = "rgba(255,200,120,0.95)";
  for (let i = 0; i < N; i++) {
    const px = ix + 10 + ((proj[i] - pmin) / span) * (insetW - 20);
    ctx.beginPath();
    ctx.arc(px, axisY, 2.5, 0, Math.PI * 2);
    ctx.fill();
  }
  ctx.fillStyle = "rgba(220,230,250,0.85)";
  ctx.font = "11px ui-monospace, monospace";
  ctx.fillText("PC1 projection", ix + 8, iy + 14);

  // HUD.
  const total = l1 + l2;
  const ratio = total > 1e-6 ? (l1 / total) : 1;
  const angleDeg = (Math.atan2(v1y, v1x) * 180 / Math.PI).toFixed(1);
  ctx.fillStyle = "rgba(220,230,250,0.95)";
  ctx.font = "13px ui-monospace, monospace";
  ctx.fillText(`λ₁ = ${l1.toFixed(0)}   σ₁ = ${s1.toFixed(1)}`, 12, 20);
  ctx.fillText(`λ₂ = ${l2.toFixed(0)}   σ₂ = ${s2.toFixed(1)}`, 12, 38);
  ctx.fillText(`var explained: ${(ratio * 100).toFixed(1)}%`, 12, 56);
  ctx.fillText(`PC1 angle: ${angleDeg}°`, 12, 74);
  ctx.fillStyle = "rgba(255,160,90,0.9)";
  ctx.fillText("— PC1", width - 132, 20);
  ctx.fillStyle = "rgba(120,220,180,0.9)";
  ctx.fillText("— PC2", width - 80, 20);
  ctx.fillStyle = "rgba(200,210,230,0.55)";
  ctx.fillText("drag points to reshape", width - 180, height - 14);
}