The algorithm the libraries actually run
The foundation of this site computes eigenvalues with Jacobi rotations, and the choice was deliberate. Jacobi is provably convergent, every step is a two-by-two problem checkable by hand, and it needs no theory to justify — which makes it exactly right for a site that wants its machinery inspectable.
It is also not what any library runs. LAPACK runs the QR algorithm, and the gap between the two is this essay.
The algorithm is three lines. Factorise A = QR. Form A′ = RQ. Repeat. Nothing about that suggests eigenvalues should appear, and they do, because A′ = QᵀAQ is a similarity transformation — so the eigenvalues never change while the matrix is ground towards triangular form, and the diagonal of what remains is the spectrum.
Which is elegant, and is not yet an algorithm.
Why the subdiagonal and not the eigenvalue
The quantity plotted needs justifying, because it is not the obvious one.
The algorithm converges when the subdiagonal entries go to zero: at that point the matrix is triangular and its diagonal is the answer. So the subdiagonal entry is what the convergence theorem is about, and it is the natural thing to watch.
It is also the only thing that shows what is happening. Cubic convergence is invisible in the eigenvalue. Three iterations taking a subdiagonal entry from 10⁻¹ to 10⁻³ to 10⁻⁹ to 10⁻²⁷ appear, on a plot of the eigenvalue’s error, as three points that are all at machine precision — the eigenvalue is correct to sixteen digits after the second step and cannot get more correct. The interesting behaviour happens below the floor of the obvious measurement.
The rate, and the two routes to it
Unshifted, the algorithm converges linearly, and the rate is the ratio of the two eigenvalues the subdiagonal entry separates: .
That is a prediction available before anything runs, since the matrices here are built by choosing the spectrum. And the run supplies the second route.
On a symmetric 4×4 matrix with eigenvalues 8, 4, 2 and 1.8 — so λ₄/λ₃ = 0.9000 — the measured contraction, fitted over the tail of the first deflation, is 0.9000. Four significant figures, two routes, no shared arithmetic: one number comes from the spectrum the matrix was assembled from, the other from watching a subdiagonal entry shrink.
The measurement took two attempts, and the failure is worth recording because it is a real trap. The algorithm’s history of subdiagonal magnitudes runs on past each deflation into the next subdiagonal, which starts large again. Fitting a rate to the concatenated history therefore measures whichever convergence contributed the most points — at λ₄/λ₃ = 0.3 that is the second subdiagonal, and the fit came back 0.5, which is λ₃/λ₂. A perfectly plausible number, for the wrong pair of eigenvalues.
The library now returns the history split per deflation, and the rate is fitted on the first segment alone. A single flat array of numbers was hiding the fact that it described several different convergences end to end.
What linear convergence at that rate actually costs
The rate is a ratio of eigenvalues, which means it can be arbitrarily close to 1, which means the algorithm can be arbitrarily slow. That is not a worst case to be noted and moved past — it is the ordinary situation, since eigenvalues of a real matrix are frequently close together.
Measured on the same matrix at two spectra:
| λ₄/λ₃ | Unshifted iterations | Wilkinson-shifted |
|---|---|---|
| 0.333 | 54 | 7 |
| 0.900 | 248 | 6 |
The first row is the case a textbook draws and it understates the difference badly — a factor of eight looks like a useful optimisation. The second row is a factor of 41, and the gap keeps widening as the ratio approaches 1 while the shifted count stays where it is.
That is the difference between an algorithm whose cost depends on the input and one whose cost does not.
The shift
The idea is a single line and it is the whole of what makes the method practical.
Instead of factorising A, factorise A − μI, and restore the shift afterwards. The result is still a similarity, so the eigenvalues are still preserved. What changes is the rate: the convergence is now governed by the ratio of distances from the shift to the nearest eigenvalues rather than by the ratio of the eigenvalues themselves. Choose μ close to an eigenvalue and that ratio is tiny.
The obvious choice is the corner entry — the Rayleigh shift — which is an estimate of the eigenvalue being converged to, improving as it converges. It gives quadratic convergence and can stall.
Wilkinson’s shift takes the eigenvalue of the trailing 2×2 block closer to the corner entry. It cannot stall on a symmetric matrix, and its convergence order there is cubic. The site fits the order from the recorded history and measures 2.87.
Fitting an order from four or five points spanning fifteen decades is not a precision measurement, and the assertion is written accordingly: the order must exceed 2.2, which distinguishes cubic from quadratic and does not pretend to more resolution than the data supports.
Where the shift comes from, and why the corner entry is a sensible guess
The choice of μ looks circular the first time — the shift should be close to an eigenvalue, and finding eigenvalues is the problem — so it is worth saying why it is not.
The bottom-right entry of a Hessenberg matrix undergoing this iteration is already an estimate of the eigenvalue the algorithm is converging to, and it improves as the iteration proceeds. Using it as the shift is a bootstrapping argument: a mediocre estimate produces a better one, which produces a better shift, and the loop closes on itself. This is the Rayleigh quotient iteration in disguise, and it is where the quadratic convergence comes from — the error in the shift and the error in the eigenvalue are the same quantity, so each step squares it.
Wilkinson’s refinement fixes the case where that bootstrap stalls. The corner entry can sit exactly between two eigenvalues placed symmetrically about it — which is precisely the ±1 matrix below, where the corner entry is 0 and the eigenvalues are ±1 — and then the estimate never improves. Taking instead the eigenvalue of the trailing 2×2 block breaks the symmetry, because a 2×2 eigenvalue problem is solved exactly by a quadratic formula rather than estimated.
So the shift is not guessed and it is not free: each step costs the solution of a 2×2 eigenvalue problem, which is a square root and some arithmetic on numbers already in hand. Against an iteration that costs O(n²), it is nothing.
The refusal, which is the reason shifts exist
The strongest available demonstration is the smallest matrix in this essay.
Take the 2×2 permutation that swaps two coordinates. Its eigenvalues are +1 and −1: same modulus, ratio exactly 1. The unshifted rate is 1, so the theory predicts no convergence at all — not slow convergence, none.
Run for two hundred iterations. The subdiagonal entry moves by 1.09·10⁻¹⁴, which is rounding. It does not decrease, it does not oscillate towards anything, it sits exactly where it started while the algorithm faithfully performs two hundred similarity transformations.
Apply a Wilkinson shift and it converges in one iteration, returning −1 and +1 to 10⁻¹².
The assertion that the unshifted algorithm converged is fed this run and must throw. It does. That refusal is the essay’s centre: the shift is not an optimisation applied to a working algorithm, it is the ingredient that makes it an algorithm.
The other two things the three-line description omits
Shifts are one. The other two are why the algorithm is affordable, and they belong to the form that makes it affordable — but the shape is worth stating here since it completes the account.
Hessenberg reduction first. Reduce the matrix to one nonzero subdiagonal by Householder similarities, once, at O(n³). The form is preserved by a QR step, so every subsequent iteration costs O(n²) rather than O(n³). Skipping it changes no answer and makes the algorithm unusable.
Deflation. When a subdiagonal entry becomes negligible the problem splits, and the algorithm continues on the leading block only. Without it every iteration would cost the full n² however much of the matrix had already converged.
Together with shifts, those three take the method from provably convergent given enough iterations to twelve iterations for six eigenvalues, which is what the site measures — two per eigenvalue, and the number quoted in every reference for a reason.
The non-symmetric case, in brief
Everything measured here is symmetric, and the general case differs in ways worth naming since the libraries have to handle it.
Real eigenvalues are not guaranteed. A real matrix can have complex conjugate pairs, and a real shift cannot converge to one — the subdiagonal entry belonging to a complex pair does not go to zero because there is no real triangular form to reach. The answer is the Francis double shift: apply two conjugate shifts at once, which is a real operation whose product is real, and let the algorithm converge to a quasi-triangular form with 2×2 blocks on the diagonal where the complex pairs live.
The convergence order drops. Wilkinson’s cubic convergence is a symmetric-case result. The general case gives quadratic, which is still excellent and is not the same claim.
And the answer is more fragile. A symmetric matrix’s eigenvalues move by at most the size of a perturbation, which is Weyl’s theorem and is checked in symmetry is worth more than precision with no constant. A defective matrix’s can move by a fractional power of it — the site measures an 8×8 Jordan block whose eigenvalues shift by 10⁻² under a perturbation at the level of storing the matrix. So for a non-symmetric problem the algorithm can be perfect and the answer still wrong, and whose fault that is belongs to the problem rather than to the iteration.
That last point is the site’s spine arriving in the eigenvalue problem. The QR algorithm is backward stable: it returns the exact eigenvalues of a matrix near the one it was given. Whether those are near the eigenvalues of the intended matrix is a question about conditioning, and it has nothing to do with the algorithm at all.
What is asserted here
The unshifted rate is the ratio of the two eigenvalues it separates, to within 0.08, at every position of the slider — the two-routes claim.
The shifted algorithm converges, and in fewer iterations than the unshifted one, everywhere.
Both shift strategies return the eigenvalues the matrix was built with, to 10⁻⁹.
And against Jacobi, in the library check: every eigenvalue agrees between the QR algorithm and the foundation’s Jacobi solver to 10⁻¹¹, which is a comparison between two entirely different algorithms sharing no code.
The refusal: the claim that the unshifted algorithm converges is fed the ±1 matrix and must throw. It does.
The connection to the power method
One more framing, because it explains why the algorithm works at all rather than merely that it does.
The QR algorithm is the power method run on every direction at once. The power method takes a vector, multiplies by A repeatedly, and converges to the dominant eigenvector at a rate |λ₂/λ₁| — the same ratio that governs the unshifted algorithm here, which is not a coincidence.
Doing it for a whole basis at once is subspace iteration: multiply an orthonormal basis by A, re-orthonormalise, repeat. The re-orthonormalisation is a QR factorisation, and unwinding the algebra shows that subspace iteration on the identity produces exactly the sequence of matrices the QR algorithm produces. The three-line recurrence is a compressed form of running n power methods simultaneously and keeping them from collapsing onto each other.
That framing explains the failure on the ±1 matrix immediately. The power method does not converge when the two largest eigenvalues have equal modulus — it oscillates between them forever — and the QR algorithm inherits the failure exactly. The same defect appears elsewhere on this site for the same reason: a rate that is known in advance notes that the Jacobi iteration matrix for the model problem has a spectrum symmetric about zero, so a power iteration for its spectral radius does not converge either, and the closed form has to be used instead.
Re-orthonormalisation being the load-bearing step also connects it back to the orthogonality field. Without it the columns collapse onto the dominant direction within a few steps — the same collapse a bound that holds with probability has to prevent between power iterations of the randomised range finder, and the same fix.
What this changes about the foundation
Not its correctness. Jacobi computes eigenvalues, computes them accurately, and symmetry is worth more than precision uses it for claims that stand.
What it changes is the account of why the eigenvalue problem is affordable. Reading the foundation alone would leave the impression that eigenvalues are computed by a provably convergent iteration that sweeps until the off-diagonal norm is small — true of Jacobi, and true of nothing anybody uses at scale.
The real answer is that they are computed by an algorithm which does not converge at all without a shift, whose convergence order is decided by the shift strategy, and which is affordable only because of a reduction that changes no answer. The elegance is entirely in the parts the three-line description omits, which is a shape this site keeps finding: the derivation is short and correct, and everything that makes it usable is in what the derivation does not mention.
There is a second reason worth naming, and it is not about speed. Jacobi’s accuracy on the small eigenvalues of a graded matrix is better than the QR algorithm’s — it computes each rotation from a 2×2 block without ever forming quantities that mix scales, so a tiny eigenvalue is found to high relative accuracy rather than to high absolute accuracy. For most problems that distinction does not matter and the QR algorithm’s speed decides. For a matrix whose eigenvalues span many decades and whose small ones are the answer, it can decide the other way, and Jacobi survives in LAPACK for exactly that case.
Which is the shape rank is a decision finds at the other end of the spectrum: a method that is right about the large values and wrong about the small ones is answering a question nobody asked, whenever the small ones are the point.