SVD: Circle → Ellipse

// SVD: the unit circle maps to an ellipse.
// For any 2x2 matrix A we can write A = U Σ Vᵀ where U, V are rotations
// (or reflections) and Σ = diag(σ1, σ2) with σ1 ≥ σ2 ≥ 0. Geometrically:
//   1. Vᵀ rotates the circle (no change in shape).
//   2. Σ stretches it into an axis-aligned ellipse with semi-axes σ1, σ2.
//   3. U rotates that ellipse to its final position.
// The right singular vectors v1, v2 are the preimages of the ellipse axes;
// the left singular vectors u1, u2 are the directions of the axes.
//
// User can grab either pink tip (σ₁u₁ or σ₂u₂) on the ellipse and drag it
// anywhere. We hold the right singular vectors v₁, v₂ fixed and back-compute
// A from the dragged left columns:
//      A = [σ₁u₁ | σ₂u₂] · [v₁ | v₂]ᵀ
// so every drag is a full edit of the matrix, not just one entry.

const HIT_R = 22;
const N_SAMPLES = 96;

let W = 0, H = 0;
let cx = 0, cy = 0;
let unit = 60;

// matrix A — start with something interesting
let A = [[1.2, 0.7], [0.4, 1.0]];

// which handle is being dragged: null | 1 | 2
let dragHandle = null;

function resize(width, height) {
  W = width; H = height;
  cx = W * 0.5; cy = H * 0.55;
  unit = Math.min(W, H) * 0.18;
}

function worldToScreen(x, y) {
  return [cx + x * unit, cy - y * unit];
}

function screenToWorld(px, py) {
  return [(px - cx) / unit, -(py - cy) / unit];
}

// 2x2 SVD via eigendecomposition of AᵀA.
function svd2(M) {
  const a = M[0][0], b = M[0][1], c = M[1][0], d = M[1][1];
  // AᵀA = [[a²+c², ab+cd], [ab+cd, b²+d²]]
  const e = a * a + c * c;
  const f = b * b + d * d;
  const g = a * b + c * d;
  const tr = e + f;
  const det = e * f - g * g;
  const disc = Math.max(0, tr * tr / 4 - det);
  const s = Math.sqrt(disc);
  const lam1 = tr / 2 + s;
  const lam2 = Math.max(0, tr / 2 - s);
  const sigma1 = Math.sqrt(Math.max(0, lam1));
  const sigma2 = Math.sqrt(Math.max(0, lam2));

  // right singular vectors v1, v2 are eigenvectors of AᵀA.
  let v1x, v1y;
  if (Math.abs(g) > 1e-9) {
    v1x = g;
    v1y = lam1 - e;
  } else if (e >= f) {
    v1x = 1; v1y = 0;
  } else {
    v1x = 0; v1y = 1;
  }
  const n1 = Math.hypot(v1x, v1y) || 1;
  v1x /= n1; v1y /= n1;
  // v2 ⟂ v1
  const v2x = -v1y, v2y = v1x;

  // u_i = A v_i / sigma_i
  let u1x, u1y, u2x, u2y;
  if (sigma1 > 1e-9) {
    u1x = (a * v1x + b * v1y) / sigma1;
    u1y = (c * v1x + d * v1y) / sigma1;
  } else { u1x = 1; u1y = 0; }
  if (sigma2 > 1e-9) {
    u2x = (a * v2x + b * v2y) / sigma2;
    u2y = (c * v2x + d * v2y) / sigma2;
  } else { u2x = -u1y; u2y = u1x; }

  return {
    sigma1, sigma2,
    v1: [v1x, v1y], v2: [v2x, v2y],
    u1: [u1x, u1y], u2: [u2x, u2y],
  };
}

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

function tick({ ctx, width, height, input }) {
  if (width !== W || height !== H) resize(width, height);

  // Compute SVD of the current matrix.
  let { sigma1, sigma2, v1, v2, u1, u2 } = svd2(A);

  // -- input: drag σ₁u₁ or σ₂u₂ tip ---------------------------------
  // Tip world-coords for both σ-vectors:
  let tip1x = u1[0] * sigma1, tip1y = u1[1] * sigma1;
  let tip2x = u2[0] * sigma2, tip2y = u2[1] * sigma2;
  const [s1px, s1py] = worldToScreen(tip1x, tip1y);
  const [s2px, s2py] = worldToScreen(tip2x, tip2y);

  if (input.mouseDown && dragHandle === null) {
    const d1 = Math.hypot(input.mouseX - s1px, input.mouseY - s1py);
    const d2 = Math.hypot(input.mouseX - s2px, input.mouseY - s2py);
    if (d1 < HIT_R && d1 <= d2) dragHandle = 1;
    else if (d2 < HIT_R) dragHandle = 2;
  }
  if (!input.mouseDown) dragHandle = null;

  if (dragHandle !== null) {
    // The dragged tip becomes the new σ_i u_i column. Hold v1, v2 fixed
    // and rebuild A as [σ₁u₁ | σ₂u₂] · [v₁ | v₂]ᵀ.
    let [wx, wy] = screenToWorld(input.mouseX, input.mouseY);
    // clamp so the user can't make the ellipse leave the visible world
    const lim = 5;
    wx = Math.max(-lim, Math.min(lim, wx));
    wy = Math.max(-lim, Math.min(lim, wy));

    if (dragHandle === 1) { tip1x = wx; tip1y = wy; }
    else                  { tip2x = wx; tip2y = wy; }

    // Build A = [tip1 tip2] · [v1 v2]ᵀ. With V orthonormal, this is
    // equivalent to assigning column i of (A V) to tip_i.
    const a00 = tip1x * v1[0] + tip2x * v2[0];
    const a01 = tip1x * v1[1] + tip2x * v2[1];
    const a10 = tip1y * v1[0] + tip2y * v2[0];
    const a11 = tip1y * v1[1] + tip2y * v2[1];
    A = [[a00, a01], [a10, a11]];

    // Recompute SVD for the freshly edited matrix so labels, ellipse,
    // and HUD all reflect the dragged state in the same frame.
    ({ sigma1, sigma2, v1, v2, u1, u2 } = svd2(A));
  }

  // -- background ----------------------------------------------------
  ctx.fillStyle = "#0a0c14";
  ctx.fillRect(0, 0, W, H);

  // grid
  ctx.strokeStyle = "rgba(120,140,180,0.10)";
  ctx.lineWidth = 1;
  ctx.beginPath();
  for (let i = -6; i <= 6; i++) {
    const [x1, y1] = worldToScreen(i, -6);
    const [x2, y2] = worldToScreen(i, 6);
    ctx.moveTo(x1, y1); ctx.lineTo(x2, y2);
    const [x3, y3] = worldToScreen(-6, i);
    const [x4, y4] = worldToScreen(6, i);
    ctx.moveTo(x3, y3); ctx.lineTo(x4, y4);
  }
  ctx.stroke();

  // axes
  ctx.strokeStyle = "rgba(255,255,255,0.30)";
  ctx.lineWidth = 1;
  ctx.beginPath();
  ctx.moveTo(0, cy); ctx.lineTo(W, cy);
  ctx.moveTo(cx, 0); ctx.lineTo(cx, H);
  ctx.stroke();

  // -- unit circle ---------------------------------------------------
  ctx.strokeStyle = "rgba(180,200,230,0.55)";
  ctx.lineWidth = 1.5;
  ctx.beginPath();
  ctx.arc(cx, cy, unit, 0, Math.PI * 2);
  ctx.stroke();

  // -- ellipse (image of circle) -------------------------------------
  ctx.strokeStyle = "rgba(120,220,255,0.95)";
  ctx.fillStyle = "rgba(120,220,255,0.10)";
  ctx.lineWidth = 2;
  ctx.beginPath();
  for (let i = 0; i <= N_SAMPLES; i++) {
    const th = (i / N_SAMPLES) * Math.PI * 2;
    const x = Math.cos(th), y = Math.sin(th);
    const ax = A[0][0] * x + A[0][1] * y;
    const ay = A[1][0] * x + A[1][1] * y;
    const [px, py] = worldToScreen(ax, ay);
    if (i === 0) ctx.moveTo(px, py); else ctx.lineTo(px, py);
  }
  ctx.closePath();
  ctx.fill();
  ctx.stroke();

  // -- right singular vectors v1, v2 (on the circle) -----------------
  drawArrow(ctx, ...worldToScreen(0, 0), ...worldToScreen(v1[0], v1[1]),
            "rgba(255,180,90,0.95)", 2);
  drawArrow(ctx, ...worldToScreen(0, 0), ...worldToScreen(v2[0], v2[1]),
            "rgba(255,180,90,0.55)", 1.7);

  // -- left singular vectors scaled by sigma_i (on the ellipse) ------
  drawArrow(ctx, ...worldToScreen(0, 0),
            ...worldToScreen(u1[0] * sigma1, u1[1] * sigma1),
            "rgba(255,120,200,1)", 2.4);
  drawArrow(ctx, ...worldToScreen(0, 0),
            ...worldToScreen(u2[0] * sigma2, u2[1] * sigma2),
            "rgba(255,120,200,0.7)", 2);

  // -- drag handles at the σ-vector tips -----------------------------
  const tip1 = worldToScreen(u1[0] * sigma1, u1[1] * sigma1);
  const tip2 = worldToScreen(u2[0] * sigma2, u2[1] * sigma2);
  drawHandle(ctx, tip1[0], tip1[1], dragHandle === 1, "rgba(255,150,210,1)");
  drawHandle(ctx, tip2[0], tip2[1], dragHandle === 2, "rgba(255,150,210,0.85)");

  // -- labels near tips ----------------------------------------------
  ctx.font = "12px ui-monospace, monospace";
  ctx.fillStyle = "rgba(255,180,90,0.95)";
  {
    const [px, py] = worldToScreen(v1[0] * 1.15, v1[1] * 1.15);
    ctx.fillText("v₁", px, py);
    const [px2, py2] = worldToScreen(v2[0] * 1.15, v2[1] * 1.15);
    ctx.fillText("v₂", px2, py2);
  }
  ctx.fillStyle = "rgba(255,150,210,1)";
  {
    const [px, py] = worldToScreen(u1[0] * sigma1 * 1.18, u1[1] * sigma1 * 1.18);
    ctx.fillText(`σ₁u₁`, px, py);
    const [px2, py2] = worldToScreen(u2[0] * sigma2 * 1.18, u2[1] * sigma2 * 1.18);
    ctx.fillText(`σ₂u₂`, px2, py2);
  }

  // -- HUD -----------------------------------------------------------
  ctx.fillStyle = "rgba(220,230,250,0.92)";
  ctx.font = "13px ui-monospace, monospace";
  ctx.fillText("A =", 12, 22);
  ctx.fillText(`[ ${pad(A[0][0])}  ${pad(A[0][1])} ]`, 48, 22);
  ctx.fillText(`[ ${pad(A[1][0])}  ${pad(A[1][1])} ]`, 48, 40);

  ctx.fillStyle = "rgba(180,220,255,0.92)";
  ctx.fillText("Σ = diag(σ₁, σ₂)", 12, 68);
  ctx.fillStyle = "rgba(255,150,210,0.95)";
  ctx.fillText(`σ₁ = ${sigma1.toFixed(3)}   (semi-major)`, 12, 88);
  ctx.fillText(`σ₂ = ${sigma2.toFixed(3)}   (semi-minor)`, 12, 106);
  ctx.fillStyle = "rgba(220,230,250,0.92)";
  ctx.fillText(`det(A) = σ₁ σ₂ = ${(sigma1 * sigma2).toFixed(3)}`, 12, 126);
  ctx.fillText(`cond = σ₁/σ₂ = ${sigma2 > 1e-6 ? (sigma1 / sigma2).toFixed(2) : "∞"}`, 12, 144);

  ctx.fillStyle = "rgba(200,210,230,0.55)";
  ctx.fillText("drag a pink tip to reshape A", W - 240, H - 12);
}

function pad(v) {
  const s = v.toFixed(2);
  return (s.startsWith("-") ? "" : " ") + s;
}

function drawHandle(ctx, x, y, active, color) {
  ctx.fillStyle = active ? "rgba(255,240,200,1)" : color;
  ctx.strokeStyle = "rgba(255,255,255,0.9)";
  ctx.lineWidth = 1.4;
  ctx.beginPath();
  ctx.arc(x, y, active ? 7 : 5.5, 0, Math.PI * 2);
  ctx.fill();
  ctx.stroke();
}

function drawArrow(ctx, x0, y0, x1, y1, color, lw) {
  const dx = x1 - x0, dy = y1 - y0;
  const len = Math.hypot(dx, dy);
  if (len < 1e-3) return;
  const ux = dx / len, uy = dy / len;
  const head = Math.min(12, len * 0.35);

  ctx.strokeStyle = color;
  ctx.lineWidth = lw;
  ctx.lineCap = "round";
  ctx.beginPath();
  ctx.moveTo(x0, y0);
  ctx.lineTo(x1 - ux * head * 0.5, y1 - uy * head * 0.5);
  ctx.stroke();

  ctx.fillStyle = color;
  ctx.beginPath();
  ctx.moveTo(x1, y1);
  ctx.lineTo(x1 - ux * head - uy * head * 0.5, y1 - uy * head + ux * head * 0.5);
  ctx.lineTo(x1 - ux * head + uy * head * 0.5, y1 - uy * head - ux * head * 0.5);
  ctx.closePath();
  ctx.fill();
}