Introduction

Imagine measuring a thousand people and recording their height, weight, arm span, and leg length. Those four numbers are correlated — tall people tend to be heavy and long-limbed too. A lot of that information is redundant.

Principal Component Analysis (PCA) asks a simple question: what are the directions in this data along which points spread out the most? Those directions — the principal components — turn out to carry the bulk of the information, while the remaining directions carry mostly noise.

PCA is not a new discovery. It was introduced independently by Karl Pearson in 1901 and Harold Hotelling in 1933, long before digital computers. Today it underpins face recognition, genome analysis, financial modeling, and nearly every field that deals with high-dimensional data. The algorithm is solved and efficient: a single pass through a matrix decomposition finds all principal components in polynomial time. The challenge is not finding them — it is understanding what they mean.

The formal tool is the eigendecomposition of the covariance matrix. The eigenvectors with the largest eigenvalues are the principal components; the eigenvalues tell you how much variance each direction explains. Project your data onto the top k components and you have reduced dimensionality while preserving as much variance as possible — a provably optimal linear compression.

Try It

Use the sliders to reshape the point cloud. The red arrow is the first principal component — the direction of greatest variance. The blue arrow is the second component, always perpendicular to the first.

<div class="controls">
  <label>Spread (X):: <span id="sx-val">3.0</span>
    <input type="range" id="sx" min="0.5" max="4" step="0.1" value="3.0">
  </label>
  <label>Spread (Y):: <span id="sy-val">1.0</span>
    <input type="range" id="sy" min="0.5" max="4" step="0.1" value="1.0">
  </label>
  <label>Tilt (°):: <span id="tilt-val">30</span>
    <input type="range" id="tilt" min="-90" max="90" step="5" value="30">
  </label>
</div>
<canvas id="canvas" width="480" height="280"></canvas>
<div class="stats" id="stats"></div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; margin: 0; color: #222; background: #fff; }
.controls { display: flex; flex-wrap: wrap; gap: .5rem 1rem; margin-bottom: .5rem; }
.controls label { font-size: .85rem; display: flex; flex-direction: column; gap: 2px; }
.controls input[type=range] { width: 120px; accent-color: #e63946; }
canvas { display: block; border: 1px solid #d0d7de; border-radius: 8px; background: #f8fafc; max-width: 100%; }
.stats { font-size: .85rem; margin-top: .4rem; color: #444; min-height: 1.2em; }
.stats b { color: #c92f3c; }
.stats span { color: #1d6db5; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// PCA on a synthetic 2-D point cloud.
// Points are drawn by sampling a standard normal, scaling by (sx, sy),
// then rotating by `tilt` degrees. We recompute PCA each time a slider moves.
//
// PCA steps (exact, no library):
//   1. Center the data (subtract mean).
//   2. Build the 2×2 covariance matrix.
//   3. Solve the 2×2 eigenvalue problem analytically.
//   4. The eigenvector with the larger eigenvalue is PC1.

const N_POINTS = 120;
const canvas = document.getElementById('canvas');
const ctx = canvas.getContext('2d');
const statsEl = document.getElementById('stats');

// Seeded LCG so the point cloud is reproducible
function lcg(seed) {
  let s = seed >>> 0;
  return () => { s = (1664525 * s + 1013904223) >>> 0; return s / 0xffffffff; };
}
// Box-Muller
function gaussPair(rand) {
  const u = rand(), v = rand();
  const r = Math.sqrt(-2 * Math.log(u + 1e-9));
  return [r * Math.cos(2 * Math.PI * v), r * Math.sin(2 * Math.PI * v)];
}

// Generate base points once (standard normal)
const rand = lcg(42);
const base = [];
for (let i = 0; i < N_POINTS; i++) {
  const [a, b] = gaussPair(rand);
  base.push([a, b]);
}

function getParams() {
  return {
    sx: parseFloat(document.getElementById('sx').value),
    sy: parseFloat(document.getElementById('sy').value),
    tilt: parseFloat(document.getElementById('tilt').value) * Math.PI / 180,
  };
}

function computePoints(sx, sy, tilt) {
  const cos = Math.cos(tilt), sin = Math.sin(tilt);
  return base.map(([a, b]) => {
    const x = sx * a, y = sy * b;
    return [cos * x - sin * y, sin * x + cos * y];
  });
}

// 2×2 PCA: returns {pc1:[ex,ey], pc2:[-ey,ex], lambda1, lambda2}
function pca2d(pts) {
  const n = pts.length;
  let mx = 0, my = 0;
  for (const [x, y] of pts) { mx += x; my += y; }
  mx /= n; my /= n;
  let cxx = 0, cxy = 0, cyy = 0;
  for (const [x, y] of pts) {
    const dx = x - mx, dy = y - my;
    cxx += dx * dx; cxy += dx * dy; cyy += dy * dy;
  }
  cxx /= n; cxy /= n; cyy /= n;
  // Eigenvalues of [[cxx,cxy],[cxy,cyy]]
  const tr = cxx + cyy;
  const det = cxx * cyy - cxy * cxy;
  const disc = Math.sqrt(Math.max(0, (tr / 2) * (tr / 2) - det));
  const lam1 = tr / 2 + disc;
  const lam2 = tr / 2 - disc;
  // Eigenvector for lam1
  let ex, ey;
  if (Math.abs(cxy) > 1e-10) {
    const norm = Math.sqrt((lam1 - cyy) ** 2 + cxy ** 2);
    ex = (lam1 - cyy) / norm;
    ey = cxy / norm;
  } else {
    // Diagonal cov
    ex = cxx >= cyy ? 1 : 0;
    ey = cxx >= cyy ? 0 : 1;
  }
  return { pc1: [ex, ey], pc2: [-ey, ex], lambda1: lam1, lambda2: lam2 };
}

function draw() {
  const { sx, sy, tilt } = getParams();
  const pts = computePoints(sx, sy, tilt);
  const { pc1, pc2, lambda1, lambda2 } = pca2d(pts);

  // Viewport
  const W = canvas.width, H = canvas.height;
  const cx = W / 2, cy = H / 2, scale = 40;

  ctx.clearRect(0, 0, W, H);

  // Axes (light)
  ctx.strokeStyle = '#d0d7de'; ctx.lineWidth = 1;
  ctx.beginPath(); ctx.moveTo(0, cy); ctx.lineTo(W, cy); ctx.stroke();
  ctx.beginPath(); ctx.moveTo(cx, 0); ctx.lineTo(cx, H); ctx.stroke();

  // Points
  ctx.fillStyle = 'rgba(29,61,87,0.35)';
  for (const [x, y] of pts) {
    ctx.beginPath();
    ctx.arc(cx + x * scale, cy - y * scale, 2.5, 0, 2 * Math.PI);
    ctx.fill();
  }

  // Arrow helper
  function arrow(ex, ey, len, color) {
    const x2 = cx + ex * len, y2 = cy - ey * len;
    ctx.strokeStyle = color; ctx.fillStyle = color; ctx.lineWidth = 2.2;
    ctx.beginPath(); ctx.moveTo(cx, cy); ctx.lineTo(x2, y2); ctx.stroke();
    // arrowhead
    const angle = Math.atan2(-(y2 - cy), x2 - cx);
    const ah = 10;
    ctx.beginPath();
    ctx.moveTo(x2, y2);
    ctx.lineTo(x2 - ah * Math.cos(angle - 0.4), y2 + ah * Math.sin(angle - 0.4));
    ctx.lineTo(x2 - ah * Math.cos(angle + 0.4), y2 + ah * Math.sin(angle + 0.4));
    ctx.closePath(); ctx.fill();
  }

  const totalVar = lambda1 + lambda2;
  const pct1 = totalVar > 0 ? (lambda1 / totalVar * 100).toFixed(1) : 50;
  const pct2 = totalVar > 0 ? (lambda2 / totalVar * 100).toFixed(1) : 50;
  const arrowLen = Math.sqrt(lambda1) * scale * 1.4;
  const arrowLen2 = Math.sqrt(lambda2) * scale * 1.4;

  arrow(pc2[0], pc2[1], arrowLen2, '#1d6db5');
  arrow(pc1[0], pc1[1], arrowLen, '#e63946');

  // Labels
  ctx.font = '13px system-ui';
  ctx.fillStyle = '#e63946';
  ctx.fillText('PC1 ' + pct1 + '%',
    cx + pc1[0] * (arrowLen + 12) - 10,
    cy - pc1[1] * (arrowLen + 12) + 5);
  ctx.fillStyle = '#1d6db5';
  ctx.fillText('PC2 ' + pct2 + '%',
    cx + pc2[0] * (arrowLen2 + 12) - 10,
    cy - pc2[1] * (arrowLen2 + 12) + 5);

  // Stats
  const angleDeg = (Math.atan2(pc1[1], pc1[0]) * 180 / Math.PI).toFixed(1);
  statsEl.innerHTML =
    fmt('<b>PC1</b> explains %pct1%% of variance (λ₁ = %lambda1%) · <span>PC2</span> explains %pct2%% (λ₂ = %lambda2%) · PC1 angle: %angleDeg%°', { pct1: pct1, lambda1: lambda1.toFixed(2), pct2: pct2, lambda2: lambda2.toFixed(2), angleDeg: angleDeg });

  // Log to console for node verification
  if (typeof console !== 'undefined') {
    console.log('PCA: lambda1=' + lambda1.toFixed(4) + ' lambda2=' + lambda2.toFixed(4) +
      ' pct1=' + pct1 + '% pc1=[' + pc1[0].toFixed(3) + ',' + pc1[1].toFixed(3) + ']');
  }
}

// Update labels and redraw on slider change
['sx','sy','tilt'].forEach(id => {
  const el = document.getElementById(id);
  const valEl = document.getElementById(id + '-val');
  el.addEventListener('input', () => { valEl.textContent = el.value; draw(); });
});

draw();
Result
Console

Notice what happens as you increase the correlation or tilt the cloud: the red arrow tracks the elongation exactly. When the cloud is a perfect circle, both directions carry equal variance and the arrows are arbitrary; as soon as the cloud stretches, PCA locks onto that stretch. The percentage labels show how much variance each component explains — together they always sum to 100%.

The Real Complexity

PCA sits firmly in the solved category — it is not NP-hard or undecidable. Here is the precise picture:

  • Computing PCA is polynomial. For an n × d data matrix, the covariance matrix takes O(∗nd2∗)O(*nd^{2}*) time to form, and eigendecomposition (or the equivalent Singular Value Decomposition) takes O(∗d3∗)O(*d^{3}*). For datasets with millions of points but modest dimensionality, this is fast. Randomized algorithms (Halko, Martinsson & Tropp, 2011) can approximate the top k components in O(∗ndk∗)O(*ndk*) time.
  • The solution is provably optimal. The Eckart–Young–Mirsky theorem (1936) guarantees that projecting onto the top k eigenvectors minimizes the mean squared reconstruction error among all possible rank-k projections — linear or not. There is no better linear compression.
  • The covariance matrix is symmetric and positive semi-definite, so all eigenvalues are real and non-negative and the eigenvectors are orthogonal. No numerical surprises.
  • Open challenges. The hard questions are not computational but statistical: How many components should you keep? Are the components interpretable? Sparse PCA (finding components that use few original features) is NP-hard in general. Robust PCA (separating a low-rank signal from sparse corruption) is a convex program solvable in polynomial time, but at higher constant cost.

PCA is also closely linked to dimensionality reduction and k-means clustering — projecting onto the top principal components often separates clusters that are invisible in the original high-dimensional space.

Where It Matters

PCA is one of the most widely deployed algorithms in science and engineering:

  • Face recognition (Eigenfaces): Turk and Pentland (1991) showed that human faces live in a very low-dimensional subspace. Project a face image onto the top ~50 principal components and you have a compact representation good enough to recognize identity.
  • Genomics: In a genome-wide association study (GWAS), researchers record hundreds of thousands of genetic variants per person. PCA on the genotype matrix separates populations by ancestry in the top two or three components — a striking visualization of human migration history.
  • Finance: Returns on thousands of stocks are highly correlated. The first principal component captures broad market movement (like an index); subsequent components capture sector rotations. Factor models built from PCA are standard in quantitative finance.
  • Image and signal compression: The JPEG standard's DCT is a close relative of PCA — both exploit the fact that natural images have most energy in a few basis directions.
  • Noise filtering: In spectroscopy, experimental noise distributes evenly across all directions; the signal concentrates in a few. Keeping only the top principal components discards the noise layer while preserving the signal.

Anywhere data is high-dimensional and correlated, PCA provides a principled first step. It does not require labeled data, runs in polynomial time, and offers a clear measure of information retained — making it the baseline against which more complex methods are compared. See also k-means clustering for how low-dimensional projections aid unsupervised learning.

Conclusion

PCA is a rare thing in algorithms: a problem that is both theoretically beautiful and practically indispensable. The Eckart–Young–Mirsky theorem guarantees that no linear method can compress data better for a given rank, and eigendecomposition finds that optimal projection in polynomial time.

What makes PCA deep is not its computational status — it is solved — but what it reveals: high-dimensional datasets often live near a much lower-dimensional subspace, and that subspace carries almost all the meaning. Strip away the noisy directions and the structure of the data becomes visible.

The hard part is never computing PCA. It is knowing which components to keep, what they mean, and when the linearity assumption breaks down. For those questions, PCA leads naturally to its nonlinear cousins — kernel PCA, autoencoders, and manifold methods — but that is a story for another article.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/pca/Content licensed under CC BY-NC 4.0.