Advanced Linear Algebra: Foundations to Frontiers

Subsection 11.3.2 Jacobi's method for computing the Spectral Decomposition

The oldest algorithm for computing the eigenvalues and eigenvectors of a (symmetric) matrix is due to Jacobi and dates back to 1846.

• [23] C. G. J. Jacobi, Über ein leichtes Verfahren, die in der Theorie der Säkular-störungen vorkommenden Gleichungen numerisch aufzulösen, Crelle's Journal 30, 51-94 (1846).

If we recall correctly (it has been 30 years since we read the paper in German), the paper thanks Jacobi's student Seidel for performing the calculations for a \(5 \times 5 \) matrix, related to the orbits of the planets, by hand...

This is a method that keeps resurfacing, since it parallelizes easily. The operation count tends to be higher (by a constant factor) than that of reduction to tridiagonal form followed by the tridiagonal QR algorithm.

Jacobi's original idea went as follows:

• We start with a symmetric matrix:

  \begin{equation*} A = \left( \begin{array}{c c c c} \alpha_{0,0} \amp \alpha_{0,1} \amp \alpha_{0,2} \amp \alpha_{0,3} \\ \alpha_{1,0} \amp \alpha_{1,1} \amp \alpha_{1,2} \amp \alpha_{1,3} \\ \alpha_{2,0} \amp \alpha_{2,1} \amp \alpha_{2,2} \amp \alpha_{2,3} \\ \alpha_{3,0} \amp \alpha_{3,1} \amp \alpha_{3,2} \amp \alpha_{3,3} \\ \end{array} \right). \end{equation*}

• We find the off-diagonal entry with largest magnitude. Let's say it is \(\alpha_{3,1} \text{.}\)

• We compute a Jacobi rotation so that

  \begin{equation*} \left( \begin{array}{c c} \gamma_{3,1} \amp \sigma_{3,1} \\ - \sigma_{3,1} \amp \gamma_{3,1} \end{array} \right) \left( \begin{array}{c c} \alpha_{1,1} \amp \alpha_{1,3} \\ \alpha_{3,1} \amp \alpha_{3,3} \end{array} \right) \left( \begin{array}{c c} \gamma_{3,1} \amp - \sigma_{3,1} \\ \sigma_{3,1} \amp \gamma_{3,1} \end{array} \right) = \left( \begin{array}{c c} \times \amp 0 \\ 0 \amp \times \end{array} \right), \end{equation*}

  where the \(\times\)s denote nonzero entries.

• We now apply the rotation as a unitary similarity transformation from the left to the rows of \(A \) indexed with \(1 \) and \(3 \text{,}\) and from the right to columns \(1\) and \(3 \text{:}\)

  \begin{equation*} \begin{array}{l} \left( \begin{array}{c c c c} \alpha_{0,0} \amp \times \amp \alpha_{0,2} \amp \times \\ \times \amp \times \amp \times \amp 0 \\ \alpha_{2,0} \amp \times \amp \alpha_{2,2} \amp \times \\ \times \amp 0 \amp \times \amp \times \end{array} \right) = \\ \left( \begin{array}{c c c c} 1 \amp 0 \amp 0 \amp 0 \\ 0 \amp \gamma_{3,1} \amp 0 \amp \sigma_{3,1} \\ 0 \amp 0 \amp 1 \amp 0 \\ 0 \amp -\sigma_{3,1} \amp 0 \amp \gamma_{3,1} \end{array} \right) \left( \begin{array}{c c c c} \alpha_{0,0} \amp \alpha_{0,1} \amp \alpha_{0,2} \amp \alpha_{0,3} \\ \alpha_{1,0} \amp \alpha_{1,1} \amp \alpha_{1,2} \amp \alpha_{1,3} \\ \alpha_{2,0} \amp \alpha_{2,1} \amp \alpha_{2,2} \amp \alpha_{2,3} \\ \alpha_{3,0} \amp \alpha_{3,1} \amp \alpha_{3,2} \amp \alpha_{3,3} \\ \end{array} \right) \left( \begin{array}{c c c c} 1 \amp 0 \amp 0 \amp 0 \\ 0 \amp \gamma_{3,1} \amp 0 \amp -\sigma_{3,1} \\ 0 \amp 0 \amp 1 \amp 0 \\ 0 \amp \sigma_{3,1} \amp 0 \amp \gamma_{3,1} \end{array} \right). \end{array} \end{equation*}

  The \(\times\)s here denote elements of the matrix that are changed by the application of the Jacobi rotation.

• This process repeats, reducing the off-diagonal element that is largest in magnitude to zero in each iteration.

Notice that each application of the Jacobi rotation is a unitary similarity transformation, and hence preserves the eigenvalues of the matrix. If this method eventually yields a diagonal matrix, then the eigenvalues can be found on the diagonal of that matrix. We do not give a proof of convergence here.

( 2,1 ) - ( 3,1 ) - ( 4,1 ) - ( 5,1 ) -
          ( 3,2 ) - ( 4,2 ) - ( 5,2 ) - 
                    ( 4,3 ) - ( 5,3 ) - 
                              ( 5,4 )

The key insight is that applying a Jacobi rotation to zero an element, \(\alpha_{i,j} \text{,}\) reduces the square of the Frobenius norm of the off-diagonal elements of the matrix by \(\alpha_{i,j}^2 \text{.}\) In other words, let \({\rm off}( A )\) equal the matrix \(A \) but with its diagonal elements set to zero. If \(J_{i,j} \) zeroes out \(\alpha_{i,j} \) (and \(\alpha_{j,i} \)), then

From this exercise, we learn:

• The good news: every time a Jacobi rotation is used to zero an off-diagonal element, \({\rm off} ( A )\) decreases by twice the square of that element.

• The bad news: a previously introduced zero may become nonzero in the process.

The original algorithm developed by Jacobi searched for the largest (in absolute value) off-diagonal element and zeroed it, repeating this processess until all off-diagonal elements were small. The problem with this is that searching for the largest off-diagonal element in an \(m \times m \) matrix requires \(O( m^2) \) comparisons. Computing and applying one Jacobi rotation as a similarity transformation requires \(O(m) \) flops. For large \(m \) this is not practical. Instead, it can be shown that zeroing the off-diagonal elements by columns (or rows) also converges to a diagonal matrix. This is known as the column-cyclic Jacobi algorithm. Zeroing out every pair of off-diagonal elements once is called a sweep. We illustrate this in Figure 11.3.2.1. Typically only a few sweeps (on the order of five) are needed to converge sufficiently.

We conclude by noting that the matrix \(Q \) such that \(A = Q \Lambda Q^H \) can be computed by accumulating all the Jacobi rotations (applying them to the identity matrix).