The Jacobi eigenvalue algorithm, introduced by Carl Gustav Jacob Jacobi in 1846, is the oldest numerical method for computing the complete eigensystem of a real symmetric matrix \(\bA = \bA^T \in \fR^{n \times n}\). The method systematically drives the matrix toward diagonal form by applying a sequence of orthogonal similarity transformations based on elementary plane rotations.
Although the shifted QR algorithm is computationally faster for general dense matrices, the Jacobi method remains of profound importance in scientific computing. It achieves high relative accuracy for small eigenvalues of graded positive definite matrices that QR cannot match, and its localized transformations allow for massive fine-grained parallelization on modern GPU and systolic architectures.
Givens Plane Rotations and the Zeroing Angle
At each step, the algorithm chooses an off-diagonal target entry \(a_{pq}\) with \(p < q\) and constructs a Givens plane rotation \(\bG(p, q, \theta)\) that annihilates that entry while preserving the overall spectrum of \(\bA\).
The matrix \(\bG(p, q, \theta) \in \fR^{n \times n}\) is an identity matrix except for the four entries at the intersections of rows and columns \(p\) and \(q\): \[
\begin{align}
g_{pp} = c, \qquad g_{pq} = s, \qquad g_{qp} = -s, \qquad g_{qq} = c,
\end{align}
\] where \(c = \cos\theta\) and \(s = \sin\theta\). The transformation \(\bA' = \bG^T \bA \bG\) modifies only rows and columns \(p\) and \(q\) of \(\bA\).
To zero the transformed off-diagonal entry \(a'_{pq} = 0\), we consider the \(2 \times 2\) orthogonal similarity transformation on the submatrix at indices \(\{p, q\}\): \[
\begin{align}
\begin{pmatrix} a'_{pp} & a'_{pq} \\ a'_{pq} & a'_{qq} \end{pmatrix}
=
\begin{pmatrix} c & -s \\ s & c \end{pmatrix}
\begin{pmatrix} a_{pp} & a_{pq} \\ a_{pq} & a_{qq} \end{pmatrix}
\begin{pmatrix} c & s \\ -s & c \end{pmatrix}.
\end{align}
\] Multiplying the matrices yields the off-diagonal formula \[
\begin{align}
a'_{pq} = (c^2 - s^2) a_{pq} + cs(a_{pp} - a_{qq}) = a_{pq} \cos(2\theta) + \frac{1}{2}(a_{pp} - a_{qq}) \sin(2\theta).
\end{align}
\] Setting \(a'_{pq} = 0\) yields the trigonometric condition \[
\begin{align}
\cot(2\theta) = \frac{a_{qq} - a_{pp}}{2a_{pq}} \equiv \tau.
\end{align}
\] Using the identity \(\cot(2\theta) = \frac{1 - \tan^2\theta}{2\tan\theta}\), let \(t = \tan\theta\). This produces the quadratic equation \[
\begin{align}
t^2 + 2\tau t - 1 = 0.
\end{align}
\] To minimize the perturbation to \(\bA\) and prevent catastrophic cancellation, we select the root of smaller absolute value: \[
\begin{align}
t = \frac{\operatorname{sign}(\tau)}{|\tau| + \sqrt{1 + \tau^2}}, \qquad c = \frac{1}{\sqrt{1 + t^2}}, \qquad s = ct.
\end{align}
\]
The standard quadratic formula gives roots \(t = -\tau \pm \sqrt{1 + \tau^2}\). Multiplying the root of smaller magnitude by its conjugate algebraic expression yields \[
\begin{align}
t = \frac{(-\tau \pm \sqrt{1 + \tau^2})(\tau \pm \sqrt{1 + \tau^2})}{\tau \pm \sqrt{1 + \tau^2}} = \frac{1}{\tau + \operatorname{sign}(\tau)\sqrt{1 + \tau^2}} = \frac{\operatorname{sign}(\tau)}{|\tau| + \sqrt{1 + \tau^2}}.
\end{align}
\] This formula satisfies \(|t| \leq 1\), corresponding to a rotation angle \(|\theta| \leq \pi/4\), and avoids subtracting nearly equal quantities.
Monotone Convergence of the Off-Diagonal Norm
Let the off-diagonal squared norm be defined by \[
\begin{align}
\operatorname{off}(\bA)^2 = \sum_{i \neq j} a_{ij}^2 = \|\bA\|_F^2 - \sum_{i=1}^n a_{ii}^2.
\end{align}
\] Applying the Jacobi rotation \(\bA' = \bG(p, q, \theta)^T \bA \bG(p, q, \theta)\) with the optimal angle strictly reduces the off-diagonal norm: \[
\begin{align}
\operatorname{off}(\bA')^2 = \operatorname{off}(\bA)^2 - 2a_{pq}^2.
\end{align}
\]
Because \(\bG\) is orthogonal, the total Frobenius norm is preserved: \(\|\bA'\|_F^2 = \|\bA\|_F^2\). The rotation affects only entries in rows and columns \(p\) and \(q\). For the \(2 \times 2\) submatrix, the sum of squares is invariant under orthogonal similarity: \[
\begin{align}
(a'_{pp})^2 + (a'_{qq})^2 + 2(a'_{pq})^2 = a_{pp}^2 + a_{qq}^2 + 2a_{pq}^2.
\end{align}
\] Because the rotation enforces \(a'_{pq} = 0\), the sum of squared diagonal entries increases: \[
\begin{align}
(a'_{pp})^2 + (a'_{qq})^2 = a_{pp}^2 + a_{qq}^2 + 2a_{pq}^2.
\end{align}
\] All other diagonal entries \(a_{ii}\) (\(i \neq p, q\)) are untouched. Substituting this into \(\operatorname{off}(\bA')^2 = \|\bA'\|_F^2 - \sum (a'_{ii})^2\) proves that \(\operatorname{off}(\bA')^2 = \operatorname{off}(\bA)^2 - 2a_{pq}^2\).
The Cyclic Jacobi algorithm sweeps through all \(N = n(n-1)/2\) off-diagonal index pairs \((p, q)\) with \(1 \leq p < q \leq n\) in a fixed periodic order:
Initialize \(\bA^{(0)} = \bA\) and eigenvector accumulator \(\bV^{(0)} = \bI_n\).
While \(\operatorname{off}(\bA^{(k)}) > `tol` \|\bA^{(0)}\|_F\): \begin{enumerate}
For each pair \(1 \leq p < q \leq n\): \begin{enumerate}
Compute \(\tau = \frac{a_{qq} - a_{pp}}{2a_{pq}}\), \(t = \frac{\operatorname{sign}(\tau)}{|\tau| + \sqrt{1 + \tau^2}}\), \(c = \frac{1}{\sqrt{1+t^2}}\), and \(s = ct\).
Apply the rotation to \(\bA\): \(\bA \leftarrow \bG(p, q, \theta)^T \bA \bG(p, q, \theta)\).
Accumulate the eigenvectors: \(\bV \leftarrow \bV \bG(p, q, \theta)\).
\end{enumerate} \end{enumerate} After a few initial sweeps, the off-diagonal entries converge to zero quadratically: \(\operatorname{off}(\bA^{(k+1)}) \leq C \cdot \operatorname{off}(\bA^{(k)})^2\).
Perform one Jacobi rotation by hand to zero out the off-diagonal entry of \[
\begin{align}
\bA = \begin{pmatrix} 3 & 2 \\ 2 & 6 \end{pmatrix}.
\end{align}
\] Explicitly compute \(\tau, t, c, s\), and the resulting diagonal matrix \(\bA'\).
Verify that the eigenvalues computed on the diagonal of \(\bA'\) match the exact roots of \(\det(\bA - \lambda\bI) = 0\).
Prove that for an \(n \times n\) matrix, if the largest off-diagonal entry is selected at each step (Classical Jacobi), then \(\operatorname{off}(\bA^{(k+1)})^2 \leq \left( 1 - \frac{2}{n(n-1)} \right) \operatorname{off}(\bA^{(k)})^2\).