Jacobi eigenvalue algorithm

In numerical linear algebra, the Jacobi eigenvalue algorithm is an iterative method for the calculation of the eigenvalues and eigenvectors of a real symmetric matrix. It is named after Carl Gustav Jacob Jacobi, who first proposed the method in 1846, but it only became widely used in the 1950s with the advent of computers.
This algorithm is inherently a dense matrix algorithm: it draws little or no advantage from being applied to a sparse matrix, and it will destroy sparseness by creating fill-in. Similarly, it will not preserve structures such as being banded of the matrix on which it operates.

Description

Let be a symmetric matrix, and be a Givens rotation matrix. Then:
is symmetric and similar to.
Furthermore, has entries:
where and.
Since is orthogonal, and have the same Frobenius norm , however we can choose such that, in which case has a larger sum of squares on the diagonal:
Set this equal to 0, and rearrange:
if
In order to optimize this effect, S_ij should be the off-diagonal element with the largest absolute value, called the pivot.
The Jacobi eigenvalue method repeatedly performs rotations until the matrix becomes almost diagonal. Then the elements in the diagonal are approximations of the eigenvalues of S.

Convergence

If is a pivot element, then by definition for . Let denote the sum of squares of all off-diagonal entries of. Since has exactly off-diagonal elements, we have or . Now. This implies
or ;
that is, the sequence of Jacobi rotations converges at least linearly by a factor to a diagonal matrix.
A number of Jacobi rotations is called a sweep; let denote the result. The previous estimate yields
that is, the sequence of sweeps converges at least linearly with a factor ≈ .
However the following result of Schönhage yields locally quadratic convergence. To this end let S have m distinct eigenvalues with multiplicities and let d > 0 be the smallest distance of two different eigenvalues. Let us call a number of
Jacobi rotations a Schönhage-sweep. If denotes the result then
Thus convergence becomes quadratic as soon as

Cost

Each Givens rotation can be done in steps when the pivot element p is known. However the search for p requires inspection of all N ≈ n² off-diagonal elements, which means this search dominates the overall complexity and pushes the computational complexity of a sweep in the classical Jacobi algorithm to. Competing algorithms attain complexity for a full diagonalisation.

Caching row maximums

We can reduce the complexity of finding the pivot element from O to O if we introduce an additional index array with the property that is the index of the largest element in row i, of the current S. Then the indices of the pivot must be one of the pairs. Also the updating of the index array can be done in O average-case complexity: First, the maximum entry in the updated rows k and l can be found in O steps. In the other rows i, only the entries in columns k and l change. Looping over these rows, if is neither k nor l, it suffices to compare the old maximum at to the new entries and update if necessary. If should be equal to k or l and the corresponding entry decreased during the update, the maximum over row i has to be found from scratch in O complexity. However, this will happen on average only once per rotation. Thus, each rotation has O and one sweep O average-case complexity, which is equivalent to one matrix multiplication. Additionally the must be initialized before the process starts, which can be done in n² steps.
Typically the Jacobi method converges within numerical precision after a small number of sweeps. Note that multiple eigenvalues reduce the number of iterations since.

Cyclic and parallel Jacobi

An alternative approach is to forego the search entirely, and simply have each sweep pivot every off-diagonal element once, in some predetermined order. It has been shown that this cyclic Jacobi attains quadratic convergence, just like the classical Jacobi.
The opportunity for parallelisation that is particular to Jacobi is based on combining cyclic Jacobi with the observation that Givens rotations for disjoint sets of indices commute, so that several can be applied in parallel. Concretely, if pivots between indices and pivots between indices, then from follows because in computing or the rotation only needs to access rows and the rotation only needs to access rows. Two processors can perform both rotations in parallel, because no matrix element is accessed for both.
Partitioning the set of index pairs of a sweep into classes that are pairwise disjoint is equivalent to partitioning the edge set of a complete graph into matchings, which is the same thing as edge colouring it; each colour class then becomes a round within the sweep. The minimal number of rounds is the chromatic index of the complete graph, and equals for odd but for even. A simple rule for odd is to handle the pairs and in the same round if. For even one may create rounds where a pair for goes into round and additionally a pair for goes into round. This brings the time complexity of a sweep down from to, if processors are available.
A round would consist of each processor first calculating for its rotation, and then applying the rotation from the left. Next, the processors synchronise before applying the transpose rotation from the right, and finally synchronising again. A matrix element may be accessed by two processors during a round, but not by both during the same half of this round.
Further parallelisation is possible by dividing the work for a single rotation between several processors, but that might be getting too fine-grained to be practical.

Algorithm

The following algorithm is a description of the Jacobi method in math-like notation.
It calculates a vector e which contains the eigenvalues and a matrix E which contains the corresponding eigenvectors; that is, is an eigenvalue and the column an orthonormal eigenvector for, i = 1,..., n.
procedure jacobi
var
i, k, l, m, state ∈ N
s, c, t, p, y, d, r ∈ R
ind ∈ Nⁿ
changed ∈ Lⁿ
function maxind ∈ N ! index of largest off-diagonal element in row k
m := k+1
for i := k+2 to n do
if │S_ki│ > │S_km│ then m := i endif
endfor
return m
endfunc
procedure update ! update e_k and its status
y := e_k; e_k := y+''t
if changedk'' and then changed_k := false; state := state−1
elsif and then changed_k := true; state := state+1
endif
endproc
procedure rotate ! perform rotation of S_ij, S_kl
┌ ┐ ┌ ┐┌ ┐
│S_kl│ │c −s││S_kl│
│ │ := │ ││ │
│S_ij│ │s c││S_ij│
└ ┘ └ ┘└ ┘
endproc
! init e, E, and arrays ind, changed
E := I; state := n
for k := 1 to n do ind_k := maxind; e_k := S_kk; changed_k := true endfor
while state≠0 do ! next rotation
m := 1 ! find index of pivot p
for k := 2 to n−1 do
if │S_{k ''indk''}│ > │S_{m ''indm''}│ then m := k endif
endfor
k := m; l := ind_m; p := S_kl
! calculate c = cos φ, s = sin φ
y := /2; d := │y│+√
r := √; c := d/''r; s'' := p/''r; t'' := p²/d
if y<0 then s := −s; t := −t endif
S_kl := 0.0; update; update
! rotate rows and columns k and l
for i := 1 to k−1 do rotate endfor
for i := k+1 to l−1 do rotate endfor
for i := l+1 to n do rotate endfor
! rotate eigenvectors
for i := 1 to n do
┌ ┐ ┌ ┐┌ ┐
│E_ik│ │c −s││E_ik│
│ │ := │ ││ │
│E_il│ │s c││E_il│
└ ┘ └ ┘└ ┘
endfor
! update all potentially changed ind_i
for i := 1 to n do ind_i := maxind endfor
loop
'''endproc'''

Example

Let
Then jacobi produces the following eigenvalues and eigenvectors after 3 sweeps :