Eigensolvers: Why Finding Eigenvectors Is a Whole Job

Writing down what an eigenvector is takes one line. Computing one for a matrix with a hundred thousand rows is a completely different sport: the neat formulas from eigenvectors and eigendecomposition stop being something you solve and turn into something you approximate. The tool that does the approximating is an eigensolver.

Why there's no formula. Above a small size, there's simply no clean formula to plug the matrix into — and the one route that looks like a formula is numerically shaky anyway, so real solvers don't chase it. (The reason it's a dead end: eigenvalues are the roots of a degree-N polynomial, the matrix's characteristic polynomial, and beyond degree 4 there's no general solution in radicals — the Abel–Ruffini result.) So eigensolvers don't "solve the formula." They approximate the eigenvalues and eigenvectors iteratively.

The simplest one, power iteration, is almost silly: take a random arrow, multiply it by the matrix, multiply again, and again. When one eigenvalue dominates and your starting arrow has even a little component in that direction, each multiplication exaggerates that strongest pull, so the arrow drifts and locks onto the dominant eigenvector. It's like tuning a radio dial — you're not calculating the station's frequency, you're homing in until it comes through clear. (PageRank is essentially this idea applied to a damped link-transition matrix: repeat the update until the importance scores stop changing.)

Real solvers are smarter versions of that idea, and they split by the kind of matrix:

• Dense, symmetric, and you need all of them → LAPACK-style dense solvers: first tidy the matrix into a simpler shape, then grind it down with a family of battle-tested routines. (The shape is called tridiagonal; the routines have names like QR, divide-and-conquer, or MRRR — you don't need to know them to use them.) This is the world behind numpy.linalg.eigh.
• Huge and sparse, and you need only a few → a smarter cousin of power iteration: it keeps multiplying by the matrix and homes in only on the handful of eigenvectors you asked for, never building the whole spectrum. Ask for "the 10 smallest" and that's all it computes. (These methods are called Lanczos for symmetric matrices, Arnoldi for general ones; in practice you reach them through scipy.sparse.linalg.eigsh / ARPACK.)

That second case is spectral clustering: build a sparse similarity graph, turn it into a graph Laplacian, ask an eigensolver for a few of its smallest non-trivial eigenvectors, and use those eigenvectors as new coordinates. In that space, tangled shapes — rings, crescents, two moons — can separate into blobs a plain K-means can grab.

The price tag: a full dense eigendecomposition costs roughly O(N³) time and O(N²) memory. A 100,000 × 100,000 dense matrix has 10 billion entries — about 80 GB just to store it once in double precision, before you've computed a single eigenvector. That single fact is why spectral clustering is filed under "excellent, but for small-to-medium N".

The moral outlasts any one method: the eigensolver is the quiet workhorse, and the interesting computation often happens before the step that gets the name. Spectral clustering is the clean example — it looks like clustering, but the real work already happened inside the eigensolver.