161762 Multivariate Analysis for Big Data

function eigen2x2(a, b, d) {
  const trace = a + d;
  const det   = a * d - b * b;
  const disc  = Math.sqrt(Math.max(trace * trace - 4 * det, 0));
  const l1    = (trace + disc) / 2;
  const l2    = (trace - disc) / 2;

  let v1, v2;
  if (Math.abs(b) > 1e-12) {
    v1 = [l1 - d, b];
    v2 = [l2 - d, b];
  } else {
    v1 = a >= d ? [1, 0] : [0, 1];
    v2 = a >= d ? [0, 1] : [1, 0];
  }
  const n1 = Math.hypot(v1[0], v1[1]);
  const n2 = Math.hypot(v2[0], v2[1]);
  v1 = [v1[0] / n1, v1[1] / n1];
  v2 = [v2[0] / n2, v2[1] / n2];

  return { values: [l1, l2], vectors: [v1, v2] };
}

function mean(arr) {
  return arr.reduce((s, x) => s + x, 0) / arr.length;
}

function variance(arr) {
  const m = mean(arr);
  return arr.reduce((s, x) => s + (x - m) ** 2, 0) / (arr.length - 1);
}

function covariance(xArr, yArr) {
  const mx = mean(xArr), my = mean(yArr);
  let s = 0;
  for (let i = 0; i < xArr.length; i++) s += (xArr[i] - mx) * (yArr[i] - my);
  return s / (xArr.length - 1);
}

function correlation(xArr, yArr) {
  return covariance(xArr, yArr) / Math.sqrt(variance(xArr) * variance(yArr));
}

function standardize(arr) {
  const m  = mean(arr);
  const sd = Math.sqrt(variance(arr));
  return arr.map(x => (x - m) / sd);
}

/** 95 % covariance ellipse (exact for χ²₂ because CDF = 1 − e^{−x/2}) */
function ellipsePoints(cov_aa, cov_ab, cov_bb, mx, my, nPts = 120, level = 0.95) {
  const chi2 = -2 * Math.log(1 - level);          // exact for 2 df
  const r    = Math.sqrt(chi2);
  const eig  = eigen2x2(cov_aa, cov_ab, cov_bb);
  const s1   = Math.sqrt(Math.max(eig.values[0], 0));
  const s2   = Math.sqrt(Math.max(eig.values[1], 0));
  const v1   = eig.vectors[0], v2 = eig.vectors[1];

  const pts = [];
  for (let i = 0; i <= nPts; i++) {
    const t  = (2 * Math.PI * i) / nPts;
    const ct = Math.cos(t), st = Math.sin(t);
    pts.push({
      x: mx + r * (s1 * ct * v1[0] + s2 * st * v2[0]),
      y: my + r * (s1 * ct * v1[1] + s2 * st * v2[1])
    });
  }
  return pts;
}

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

/** Generate n standard-normal pairs (Box-Muller, seeded) */
function generateNormals(n, seed) {
  const rng = mulberry32(seed);
  const z1 = [], z2 = [];
  for (let i = 0; i < n; i++) {
    const u1 = rng(), u2 = rng();
    z1.push(Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2));
    z2.push(Math.sqrt(-2 * Math.log(u1)) * Math.sin(2 * Math.PI * u2));
  }
  return { z1, z2 };
}

/** Eigendecomposition of n×n symmetric matrix — Jacobi rotations */
function eigenSymmetric(mat) {
  const n = mat.length;
  const A = mat.map(r => [...r]);
  const V = Array.from({length: n}, (_, i) =>
    Array.from({length: n}, (_, j) => i === j ? 1 : 0));
  for (let it = 0; it < 200; it++) {
    let mx = 0, p = 0, q = 1;
    for (let i = 0; i < n; i++)
      for (let j = i + 1; j < n; j++)
        if (Math.abs(A[i][j]) > mx) { mx = Math.abs(A[i][j]); p = i; q = j; }
    if (mx < 1e-12) break;
    const apq = A[p][q], d = A[p][p] - A[q][q];
    let t;
    if (Math.abs(apq) < 1e-15 * Math.abs(d)) t = apq / d;
    else {
      const phi = d / (2 * apq);
      t = 1 / (Math.abs(phi) + Math.sqrt(phi * phi + 1));
      if (phi < 0) t = -t;
    }
    const c = 1 / Math.sqrt(t * t + 1), s = t * c, tau = s / (1 + c);
    A[p][p] -= t * apq; A[q][q] += t * apq; A[p][q] = 0; A[q][p] = 0;
    for (let r = 0; r < n; r++) if (r !== p && r !== q) {
      const rp = A[r][p], rq = A[r][q];
      A[r][p] = A[p][r] = rp - s * (rq + tau * rp);
      A[r][q] = A[q][r] = rq + s * (rp - tau * rq);
    }
    for (let r = 0; r < n; r++) {
      const vp = V[r][p], vq = V[r][q];
      V[r][p] = vp - s * (vq + tau * vp);
      V[r][q] = vq + s * (vp - tau * vq);
    }
  }
  const vals = A.map((r, i) => r[i]);
  const idx = vals.map((_, i) => i).sort((a, b) => vals[b] - vals[a]);
  return { values: idx.map(i => vals[i]), vectors: idx.map(i => V.map(r => r[i])) };
}

Plotly = {
  if (window.Plotly) return window.Plotly;
  const s = document.createElement("script");
  s.src = "./plotly.min.js";
  document.head.appendChild(s);
  return new Promise((resolve, reject) => {
    s.onload = () => resolve(window.Plotly);
    s.onerror = () => reject(new Error("Failed to load Plotly from CDN"));
  });
}

Why Linear Algebra?

You will not do matrix arithmetic by hand

Contract for today: We will not compute. We will interpret pictures.

But you need to understand what multivariate tools do under the hood:

Covariance captures the shape of your data (Lecture 1)
Matrices transform point clouds
Singularity signals redundancy (collinearity)
These ideas lead directly to PCA (Lecture 3)

Matrices as operators

Key idea

A matrix is a function that takes a vector in and produces a new vector out.

\[\mathbf{y} = \mathbf{M}\,\mathbf{x}\]

Every multivariate technique — PCA, LDA, regression — applies a matrix operator to your data.

A matrix doesn’t just store numbers. It rotates, stretches, squeezes, or collapses your data cloud.

Business analogy: A matrix is a recipe that mixes inputs into new composite metrics.
Feature engineering: PCA is automated feature engineering using rotations.

Example of a matrix operation

A data observation is a vector — one column of numbers ($n \times 1$):

\[\mathbf{x} = \begin{pmatrix} 230 \\ 37 \\ 22 \end{pmatrix} {\small \begin{array}{l} \leftarrow \text{TV spend (\$k)} \\ \leftarrow \text{radio spend (\$k)} \\ \leftarrow \text{sales (\$k)} \end{array}}\]

Apply a diagonal scaling matrix $\mathbf{M}$ (standardise each variable by its SD):

\[\underbrace{\begin{pmatrix} 1/s_1 & 0 & 0 \\ 0 & 1/s_2 & 0 \\ 0 & 0 & 1/s_3 \end{pmatrix}}_{\mathbf{M}\;(3 \times 3)} \underbrace{\begin{pmatrix} 230 \\ 37 \\ 22 \end{pmatrix}}_{\mathbf{x}\;(3 \times 1)} =\underbrace{\begin{pmatrix} {\approx}1.23 \\ {\approx}0.91 \\ {\approx}{-0.45} \end{pmatrix}}_{\mathbf{y}\;(3 \times 1)\;\text{(z-scores)}}\]

Same idea in PCA, LDA, regression — just a different choice of $\mathbf{M}$ in \[\mathbf{y} = \mathbf{M}\,\mathbf{x}\]

Rotation matrices

A rotation changes our coordinate directions by an angle $\theta$, without stretching the cloud.

\[ \begin{pmatrix} x' \\ y' \end{pmatrix} = \underbrace{ \begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{pmatrix} }_{\text{rotation by }\theta} \begin{pmatrix} x \\ y \end{pmatrix} \]

$\cos\theta$ keeps the “same-axis” contribution
$\sin\theta$ mixes in the perpendicular axis
Signs control the direction of rotation
The columns are the new unit axes expressed in old coordinates

Numbers you will see in the transformation playground

Angle	$\cos\theta$	$\sin\theta$
$30°$	$0.866$	$0.500$
$60°$	$0.500$	$0.866$

Key properties

Rotation preserves shape: all distances unchanged
For a pure rotation matrix, $\det = 1$
That follows from $\cos^2\theta + \sin^2\theta = 1$
$\det = 1$ → area preserved, no reflection
Columns are orthonormal: each has length 1, they are perpendicular

Transformation playground

Select a matrix transformation and see how it reshapes the point cloud. Blue = original (standardised TV vs sales). Red = transformed. Watch: distances, area, and whether the cloud collapses.

viewof tfPreset = {
  const el = Inputs.select(
    [
      "Identity",
      "Rotation 30°",
      "Rotation 60°",
      "Scale (stretch X)",
      "Scale (squash Y)",
      "Shear (mix X with Y)",
      "Near-singular",
      "Singular (rank 1)",
      "Decorrelate"
    ],
    { label: "Transformation", value: "Identity" }
  );

  el.classList.add("tf-preset");
  return el;
}

tfMatrix = {
  // Compute PCA rotation from standardised TV vs sales
  const xRaw = advertising.map(d => d.TV);
  const yRaw = advertising.map(d => d.sales);
  const xs = standardize(xRaw);
  const ys = standardize(yRaw);
  const vx  = variance(xs);
  const vy  = variance(ys);
  const cxy = covariance(xs, ys);
  const eig = eigen2x2(vx, cxy, vy);
  // Eigenvector matrix rows = transpose of eigenvectors (V^T rotates to PC space)
  const pcaRot = [
    eig.vectors[0][0], eig.vectors[0][1],
    eig.vectors[1][0], eig.vectors[1][1]
  ];

  const presets = {
    "Identity":            [1, 0, 0, 1],
    "Rotation 30°":        [0.866, -0.5, 0.5, 0.866],
    "Rotation 60°":        [0.5, -0.866, 0.866, 0.5],
    "Scale (stretch X)":   [1.8, 0, 0, 1],
    "Scale (squash Y)":    [1, 0, 0, 0.3],
    "Shear (mix X with Y)":  [1, 0.8, 0, 1],
    "Near-singular":       [1, 1, 0.99, 1],
    "Singular (rank 1)":   [1, 2, 2, 4],
    "Decorrelate":         pcaRot
  };
  return presets[tfPreset];
}

tfDet = {
  const [a, b, c, d] = tfMatrix;
  return a * d - b * c;
}

tfData = {
  const xRaw = advertising.map(d => d.TV);
  const yRaw = advertising.map(d => d.sales);
  const xs = standardize(xRaw);
  const ys = standardize(yRaw);

  const [a, b, c, d] = tfMatrix;
  const original    = xs.map((x, i) => ({ x, y: ys[i], group: "Original" }));
  const transformed = xs.map((x, i) => ({
    x: a * x + b * ys[i],
    y: c * x + d * ys[i],
    group: "Transformed"
  }));
  return [...original, ...transformed];
}

html`<div style="display:flex;align-items:center;gap:1.2em;margin-bottom:0.2em;flex-wrap:wrap;">
  <div style="background:rgba(0,0,0,0.05);border-radius:6px;padding:0.35em 0.6em;font-family:monospace;font-size:0.88em;">
    <strong style="font-family:sans-serif;">M</strong> =
    <table style="display:inline-table;vertical-align:middle;border-collapse:collapse;margin-left:0.3em;">
      <tr>
        <td style="padding:0.15em 0.45em;border:1px solid #999;text-align:right;">${tfMatrix[0].toFixed(3)}</td>
        <td style="padding:0.15em 0.45em;border:1px solid #999;text-align:right;">${tfMatrix[1].toFixed(3)}</td>
      </tr>
      <tr>
        <td style="padding:0.15em 0.45em;border:1px solid #999;text-align:right;">${tfMatrix[2].toFixed(3)}</td>
        <td style="padding:0.15em 0.45em;border:1px solid #999;text-align:right;">${tfMatrix[3].toFixed(3)}</td>
      </tr>
    </table>
  </div>
  <span style="font-size:0.88em;padding:0.25em 0.65em;border-radius:4px;
    background:${Math.abs(tfDet) < 0.01 ? '#c0392b' : '#27ae60'};color:white;font-weight:bold;">
    det = ${tfDet.toFixed(3)}
  </span>
</div>`

Plot.plot({
  width: 380,
  height: 380,
  grid: true,
  color: {
    domain: ["Original", "Transformed"],
    range: ["steelblue", "#e74c3c"],
    legend: true
  },
  marks: [
    Plot.dot(tfData, {
      x: "x", y: "y",
      fill: "group",
      fillOpacity: 0.45,
      r: 3
    }),
    Plot.frame()
  ],
  x: { domain: [-4, 4], label: "x" },
  y: { domain: [-4, 4], label: "y" }
})

What did you notice?

Quick discussion (turn to a neighbour)

Which transformation preserved distances between points?

→ Rotation — all pairwise distances preserved (det = ±1, orthogonal matrix)

Which one collapsed the cloud to a line?

→ Singular (rank 1) — det = 0, the 2D cloud was squashed onto 1D

What was special about det = 0?

→ Information is destroyed. Different inputs map to the same output — you cannot invert the transformation.

Correlation vs covariance

Toggle standardisation to see the difference between the covariance and correlation matrices.

viewof corrVar1 = Inputs.select(columns, { label: "X variable", value: "TV" })
viewof corrVar2 = Inputs.select(columns, { label: "Y variable", value: "sales" })
viewof useCorrelation = Inputs.toggle({ label: "Standardise → correlation matrix", value: false })

corrData = {
  const xRaw = advertising.map(d => d[corrVar1]);
  const yRaw = advertising.map(d => d[corrVar2]);
  const x = useCorrelation ? standardize(xRaw) : xRaw;
  const y = useCorrelation ? standardize(yRaw) : yRaw;

  const vx  = variance(x);
  const vy  = variance(y);
  const cxy = covariance(x, y);
  const r   = correlation(x, y);
  const mx  = mean(x);
  const my  = mean(y);
  const ell = ellipsePoints(vx, cxy, vy, mx, my);
  const points = x.map((xi, i) => ({ x: xi, y: y[i] }));
  return { vx, vy, cxy, r, mx, my, ell, points };
}

{
  const container = document.createElement("div");
  container.style.display = "flex";
  container.style.gap = "1em";
  container.style.alignItems = "flex-start";

  // ── Left panel: matrix + r ──
  const left = document.createElement("div");
  left.style.flex = "0 0 240px";

  const matLabel = useCorrelation ? "Correlation Matrix  R" : "Covariance Matrix  S";
  left.innerHTML = `
    <div style="background:rgba(0,0,0,0.05);border-radius:6px;padding:0.5em 0.8em;
      font-family:monospace;font-size:0.85em;">
      <div style="margin-bottom:0.3em;font-weight:bold;font-family:sans-serif;font-size:0.9em;">
        ${matLabel}
      </div>
      <table style="border-collapse:collapse;text-align:right;width:100%;">
        <tr>
          <td></td>
          <td style="padding:0.15em 0.4em;font-family:sans-serif;font-weight:bold;font-size:0.78em;">${useCorrelation ? "z("+corrVar1+")" : corrVar1}</td>
          <td style="padding:0.15em 0.4em;font-family:sans-serif;font-weight:bold;font-size:0.78em;">${useCorrelation ? "z("+corrVar2+")" : corrVar2}</td>
        </tr>
        <tr>
          <td style="padding:0.15em 0.4em;font-family:sans-serif;font-weight:bold;font-size:0.78em;">${useCorrelation ? "z("+corrVar1+")" : corrVar1}</td>
          <td style="padding:0.2em 0.4em;border:1px solid #999;">${corrData.vx.toFixed(3)}</td>
          <td style="padding:0.2em 0.4em;border:1px solid #999;">${corrData.cxy.toFixed(3)}</td>
        </tr>
        <tr>
          <td style="padding:0.15em 0.4em;font-family:sans-serif;font-weight:bold;font-size:0.78em;">${useCorrelation ? "z("+corrVar2+")" : corrVar2}</td>
          <td style="padding:0.2em 0.4em;border:1px solid #999;">${corrData.cxy.toFixed(3)}</td>
          <td style="padding:0.2em 0.4em;border:1px solid #999;">${corrData.vy.toFixed(3)}</td>
        </tr>
      </table>
    </div>
    <div style="background:#c0392b;color:white;border-radius:6px;padding:0.45em 0.8em;
      font-size:1.05em;font-weight:bold;margin-top:0.6em;text-align:center;">
      Pearson <em>r</em> = ${corrData.r.toFixed(3)}
    </div>
  `;

  // ── Right panel: plot ──
  const right = document.createElement("div");
  right.style.flex = "1";
  const plot = Plot.plot({
    width: 460,
    height: 310,
    grid: true,
    marks: [
      Plot.dot(corrData.points, {
        x: "x", y: "y",
        fill: "steelblue", fillOpacity: 0.5, r: 3
      }),
      Plot.line(corrData.ell, {
        x: "x", y: "y",
        stroke: "#c0392b", strokeWidth: 2.5
      }),
      Plot.frame()
    ],
    x: { label: useCorrelation ? `z(${corrVar1})` : corrVar1 },
    y: { label: useCorrelation ? `z(${corrVar2})` : corrVar2 }
  });
  right.appendChild(plot);

  container.appendChild(left);
  container.appendChild(right);
  return container;
}

The ellipse tilt stays the same — correlation is scale-free, but the axes become in standard deviation units.
S for scale-dependent covariance; R for scale-free correlation.

What is an eigenvector?

A direction the matrix only stretches — never rotates

For a square matrix S, an eigenvector $\mathbf{v}$ satisfies:

\[\mathbf{S}\,\mathbf{v} = \lambda\,\mathbf{v}\]

$\mathbf{v}$ is the direction (eigenvector) — the matrix doesn’t rotate it
$\lambda$ is the eigenvalue — how much the matrix stretches along that direction
A 2×2 covariance matrix has two eigenvectors → the axes of the ellipse
Large $\lambda$ = lots of spread; tiny $\lambda$ = almost no information that way

Business intuition: eigenvectors are the natural “axes” of your data cloud. They become the principal components in Lecture 3.

When variables are redundant

Drag the slider towards 1 and watch the ellipse collapse → singularity. Toggle eigenvectors to see the axes shrink.

viewof corrStrength = Inputs.range([0, 0.99], {
  label: "Correlation strength (r)",
  value: 0.30,
  step: 0.01
})
viewof showSingEigen = Inputs.toggle({ label: "Show eigenvectors", value: false })

singData = {
  const r = corrStrength;
  const { z1, z2 } = baseNormals;
  const x = z1;
  const y = z1.map((z, i) => r * z + Math.sqrt(Math.max(1 - r * r, 0)) * z2[i]);

  const vx  = variance(x);
  const vy  = variance(y);
  const cxy = covariance(x, y);
  const mx  = mean(x);
  const my  = mean(y);
  const ell = ellipsePoints(vx, cxy, vy, mx, my);
  const det = vx * vy - cxy * cxy;
  const eig = eigen2x2(vx, cxy, vy);
  const condNumber = eig.values[0] / Math.max(eig.values[1], 1e-10);
  const rSample = correlation(x, y);
  const points = x.map((xi, i) => ({ x: xi, y: y[i] }));
  return { det, condNumber, rSample, ell, points, eig, vx, vy, cxy, mx, my };
}

html`<div style="display:flex;gap:0.8em;margin-bottom:0.2em;flex-wrap:wrap;">
  <span style="font-size:0.82em;padding:0.2em 0.6em;border-radius:4px;
    background:${singData.det < 0.05 ? '#c0392b' : '#27ae60'};color:white;">
    det(S) = ${singData.det.toFixed(4)}
  </span>
  <span style="font-size:0.82em;padding:0.2em 0.6em;border-radius:4px;background:rgba(0,0,0,0.07);">
    r = ${singData.rSample.toFixed(3)}
  </span>
  <span style="font-size:0.82em;padding:0.2em 0.6em;border-radius:4px;background:rgba(0,0,0,0.07);">
    condition (λ₁/λ₂) = ${singData.condNumber.toFixed(1)}
  </span>
  ${singData.det < 0.05
      ? html`<span style="font-size:0.82em;padding:0.2em 0.6em;border-radius:4px;
              background:#c0392b;color:white;font-weight:bold;">⚠ Redundant variables warning <span style="font-weight:normal;font-size:0.9em;">(near-singular)</span></span>`
      : html``}
</div>`

{
  const marks = [
    Plot.dot(singData.points, {
      x: "x", y: "y",
      fill: "steelblue", fillOpacity: 0.5, r: 3
    }),
    Plot.line(singData.ell, {
      x: "x", y: "y",
      stroke: "#c0392b", strokeWidth: 2.5
    })
  ];

  if (showSingEigen) {
    const { eig, mx, my } = singData;
    const scale1 = Math.sqrt(Math.max(eig.values[0], 0)) * 2;
    const scale2 = Math.sqrt(Math.max(eig.values[1], 0)) * 2;

    marks.push(Plot.arrow(
      [{
        x1: mx - scale1 * eig.vectors[0][0],
        y1: my - scale1 * eig.vectors[0][1],
        x2: mx + scale1 * eig.vectors[0][0],
        y2: my + scale1 * eig.vectors[0][1]
      }],
      { x1: "x1", y1: "y1", x2: "x2", y2: "y2",
        stroke: "#27ae60", strokeWidth: 3 }
    ));
    marks.push(Plot.arrow(
      [{
        x1: mx - scale2 * eig.vectors[1][0],
        y1: my - scale2 * eig.vectors[1][1],
        x2: mx + scale2 * eig.vectors[1][0],
        y2: my + scale2 * eig.vectors[1][1]
      }],
      { x1: "x1", y1: "y1", x2: "x2", y2: "y2",
        stroke: "#8e44ad", strokeWidth: 3 }
    ));
  }

  marks.push(Plot.frame());

  return Plot.plot({
    width: 350,
    height: 350,
    grid: true,
    marks,
    x: { domain: [-3.5, 3.5], label: "x" },
    y: { domain: [-3.5, 3.5], label: "y" }
  });
}

showSingEigen
  ? html`<div style="font-size:0.80em;margin-top:0.1em;">
      <span style="color:#27ae60;">■</span> <strong>Eigenvector 1</strong> (λ₁ = ${singData.eig.values[0].toFixed(3)})
      &nbsp;&nbsp;
      <span style="color:#8e44ad;">■</span> <strong>Eigenvector 2</strong> (λ₂ = ${singData.eig.values[1].toFixed(3)})
      &nbsp;&nbsp;—&nbsp; as r → 1, λ₂ → 0 (the purple axis vanishes)
    </div>`
  : html``

What a degenerate matrix looks like

A concrete numerical example

Suppose a dataset records TV spend and TV spend doubled as a second column:

Obs	TV	TV×2	sales
1	100	200	12
2	230	460	22
3	45	90	8
4	180	360	18

Key insight: Column 2 is $\text{TV} \times 2 = 2 \cdot \text{TV}$ — no new information. This forces Row 2 = 2 × Row 1 in the covariance matrix below.

The $(3 \times 3)$ covariance matrix collapses: \[ \mathbf{S} = \begin{pmatrix} 4{,}900 & 9{,}800 & 280 \\ 9{,}800 & 19{,}600 & 560 \\ 280 & 560 & 33 \end{pmatrix} \]

Row 2 = 2 × Row 1: no independent second dimension
$\det(\mathbf{S}) = 0$: singular, cannot be inverted
One eigenvalue = 0: a direction with zero variance

Why singularity matters

When det(S) → 0

The covariance matrix cannot be inverted
Parameter estimates become unstable (huge standard errors)
Effective dimensionality is reduced — redundant variables
Near-zero eigenvalues = directions without variance

Business translation: If two KPIs move in lockstep you are measuring the same thing twice. Drop one — or let a multivariate technique handle it.

Advertising in 3D — the full point cloud (time permitting)

PCA preview: PCA finds a new axis system where variables are uncorrelated, ordered by variance.
“Project” means drop a perpendicular shadow onto a line.

Rotate and zoom to explore how three variables relate simultaneously. Toggle Show eigenvectors to draw the three PC lines — each point projected onto the PC axis, mapped back into original TV/radio/sales space.

viewof show3DEigen = Inputs.toggle({ label: "Show eigenvectors", value: false })

data3D = {
  const x = advertising.map(d => d.TV);
  const y = advertising.map(d => d.radio);
  const z = advertising.map(d => d.sales);
  const vars = [x, y, z];
  const S = vars.map((vi) => vars.map((vj) => covariance(vi, vj)));
  const means = vars.map(v => mean(v));
  const eig = eigenSymmetric(S);
  const total = eig.values.reduce((a, b) => a + b, 0);
  const pctExplained = eig.values.map(v => (100 * v / total));
  return { x, y, z, S, means, eig, pctExplained };
}

{
  const div = document.createElement("div");
  div.style.width = "100%";
  div.style.height = "480px";
  div.style.touchAction = "none";
  div.setAttribute("data-prevent-swipe", "");

  const traces = [{
    x: data3D.x, y: data3D.y, z: data3D.z,
    mode: "markers",
    type: "scatter3d",
    marker: { size: 2.5, color: "steelblue", opacity: 0.6 },
    name: "Observations"
  }];

  // Center mark (always shown)
  traces.push({
    x: [data3D.means[0]], y: [data3D.means[1]], z: [data3D.means[2]],
    mode: "markers",
    type: "scatter3d",
    marker: { size: 6, color: "red", symbol: "cross", opacity: 1 },
    name: "Center",
    showlegend: true
  });

  if (show3DEigen) {
    const colors = ["#D55E00", "#0072B2", "#009E73"];
    const labels = ["PC1", "PC2", "PC3"];
    const widths  = [6, 4, 3];

    const vectors = data3D.eig.vectors;   // already column-extracted & sorted ↓ λ
    const values  = data3D.eig.values;
    const m       = data3D.means;
    const n       = data3D.x.length;
    const dot     = (a, b) => a.reduce((s, ai, i) => s + ai * b[i], 0);

    // Center the raw data once
    const centered = Array.from({ length: n }, (_, i) => [
      data3D.x[i] - m[0],
      data3D.y[i] - m[1],
      data3D.z[i] - m[2]
    ]);

    const total = values.reduce((a, b) => a + b, 0);

    for (let k = 0; k < 3; k++) {
      const v = vectors[k];
      // PC score for each observation: how far along eigenvector k
      const scores = centered.map(p => dot(p, v));
      // Sort by score so Plotly draws a continuous line (not a zigzag)
      const idx = scores.map((s, i) => ({ s, i })).sort((a, b) => a.s - b.s).map(d => d.i);
      // Project back into original 3D space  (R: PC1_DRESS = center + score * loading)
      traces.push({
        x: idx.map(i => m[0] + scores[i] * v[0]),
        y: idx.map(i => m[1] + scores[i] * v[1]),
        z: idx.map(i => m[2] + scores[i] * v[2]),
        mode: "lines",
        type: "scatter3d",
        line: { color: colors[k], width: widths[k] },
        name: `${labels[k]} (λ=${values[k].toFixed(1)}, ${(100 * values[k] / total).toFixed(0)}%)`
      });
    }
  }

  const layout = {
    autosize: true,
    margin: { l: 0, r: 0, t: 10, b: 0 },
    scene: {
      aspectmode: "manual",
      aspectratio: { x: 1.2, y: 1.0, z: 0.8 },
      xaxis: { title: "TV ($k)" },
      yaxis: { title: "radio ($k)" },
      zaxis: { title: "sales" },
      camera: { eye: { x: 1.5, y: 1.5, z: 1.0 } }
    },
    legend: { x: 0, y: 1, font: { size: 11 } },
    paper_bgcolor: "rgba(0,0,0,0)",
    plot_bgcolor: "rgba(0,0,0,0)"
  };

  Plotly.newPlot(div, traces, layout, { responsive: true, displayModeBar: false, scrollZoom: true });
  return div;
}

show3DEigen
  ? html`<div style="font-size:0.80em;margin-top:0.1em;">
      Eigenvalues of 3×3 covariance:
      <span style="color:#D55E00;font-weight:bold;">PC1 λ₁ = ${data3D.eig.values[0].toFixed(1)} (${data3D.pctExplained[0].toFixed(0)}%)</span> |
      <span style="color:#0072B2;font-weight:bold;">PC2 λ₂ = ${data3D.eig.values[1].toFixed(1)} (${data3D.pctExplained[1].toFixed(0)}%)</span> |
      <span style="color:#009E73;font-weight:bold;">PC3 λ₃ = ${data3D.eig.values[2].toFixed(1)} (${data3D.pctExplained[2].toFixed(0)}%)</span>
    </div>`
  : html``

Wrapping Up

When $r \to 1$, what happens geometrically, and what does it imply about information?

→ The ellipse collapses toward a line; one eigenvalue goes to $\approx 0$ — one direction contains almost no independent information.

Which methods become unstable when redundancy is extreme, and why?

→ Anything needing $\mathbf{S}^{-1}$ becomes unstable or undefined: inversion amplifies noise along tiny-variance directions.

What does standardisation change, and what does it preserve?

→ It rescales axes to unit variance, turning covariance into correlation, while preserving the relationship pattern (the ellipse tilt).

Next: Lecture 3 — Principal Component Analysis

What’s coming

Eigendecomposition: $\mathbf{S} = \mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^\top$
The eigenvectors you just saw become the principal components
Covariance or correlation matrix? Which to feed PCA and why
Choosing how many components to keep (scree plot, Kaiser rule)
Business applications: dimensionality reduction, visualisation, feature engineering

Angle	\(\cos\theta\)	\(\sin\theta\)
\(30°\)	\(0.866\)	\(0.500\)
\(60°\)	\(0.500\)	\(0.866\)