Jacobi's eigenvalue method

A symmetric matrix $A$ has an orthogonal matrix of eigenvectors $V$ with $V^{T}AV$ diagonal. Jacobi's method (1846) builds $V$ out of plane rotations. A rotation in the $(p,q)$ plane, through an angle set by $a_{pp}$, $a_{qq}$ and $a_{pq}$ (derived under "One rotation"), makes the entry $a_{pq}$ zero; doing that for every pair $p<q$ once is a sweep. A later rotation can make an earlier zero nonzero again, so one sweep does not finish the job, but the size of the off-diagonal part, $\mathrm{off}(A)=({\sum }_{i\ne j}{a}_{ij}^2{)}^{1/2}$, falls with every rotation.

The cells below run it on a fixed 8 by 8 matrix and measure $\mathrm{off}(A)$ after each sweep. The reference is numpy: the same sweeps in Python and numpy.linalg.eigh on the same matrix (posts/refs/jacobi-eigenvalue-sweeps.py, output beside it).

One rotation

Rotating the $p$ and $q$ rows and columns through $\theta$ changes only those two rows and columns. The new $(p,q)$ entry is

which is zero when $\tan 2\theta =2a_{pq}/(a_{pp}-a_{qq})$. atan2 picks the angle with $|\theta |\le \pi /4$ when $a_{pp}>a_{qq}$, the smaller rotation, which disturbs the rest of the matrix least.

A rotation leaves the sum of squares of all entries unchanged, and it moves exactly $2a_{pq}^2$ of it from off the diagonal onto it. That is why $\mathrm{off}(A)$ can only fall. The cell checks both facts on the first rotation, $(p,q)=(0,1)$.

Sweeps

A sweep rotates every pair in order, $(0,1),(0,2),\dots ,(6,7)$: 28 rotations for $N=8$. The cell runs six sweeps and records $\mathrm{off}(A)$ after each.

While the off-diagonal part is large, a sweep removes a fixed fraction of it: the first two take it from 1.33 to 0.453 to 0.0299. Once it is small next to the gaps between the eigenvalues (here the smallest gap is 0.44), each sweep roughly squares it: $0.0299^2\approx 9\times 10^{-4}$ and the next value is $3.2\times 10^{-5}$; $(3.2\times 10^{-5}{)}^2\approx 10^{-9}$ and the next is $8\times 10^{-13}$. The number of correct digits doubles each sweep, as it does for Newton's method, and the fifth sweep hits rounding: under $10^{-15}$, about $4.5\times 10^{-17}$ times the size of the matrix itself. Through the fourth sweep every value agrees with the numpy twin, which does the same rotations in the same order, to $10^{-15}$; after that both are rounding noise, $4\times 10^{-16}$ here and $5\times 10^{-16}$ in numpy.

The answer

After the sweeps the diagonal of $A$ holds the eigenvalues and the columns of $V$ the eigenvectors: $V^T{A}_{0}V$ is diagonal. The eigenvalues are checked against numpy.linalg.eigh, which uses a different algorithm (LAPACK's tridiagonal reduction), and the product $V\mathrm{diag}(\lambda )V^T$ is checked against the matrix we started from.