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 (a process known as diagonalization). It is named after Carl Gustav Jacob Jacobi, who first proposed the method in 1846, but only became widely used in the 1950s with the advent of computers.
Description.
Let formula_1 be a symmetric matrix, and formula_2 be a Givens rotation matrix. Then:
is symmetric and similar to formula_1.
Furthermore, formula_5 has entries:
where formula_7 and formula_8.
Since formula_9 is orthogonal, formula_1 and formula_5 have the same Frobenius norm formula_12 (the square-root sum of squares of all components), however we can choose formula_13 such that formula_14, in which case formula_5 has a larger sum of squares on the diagonal:
Set this equal to 0, and rearrange:
if formula_18
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 (real) eigenvalues of "S".
Convergence.
If formula_20 is a pivot element, then by definition formula_21 for formula_22 . Let formula_23 denote the sum of squares of all off-diagonal entries of formula_1. Since formula_1 has exactly formula_26 off-diagonal elements, we have formula_27 or formula_28 . Now formula_29. This implies
 formula_30 or formula_31,
i.e. the sequence of Jacobi rotations converges at least linearly by a factor formula_32 to a diagonal matrix.
A number of formula_33 Jacobi rotations is called a sweep; let formula_34 denote the result. The previous estimate yields
i.e. the sequence of sweeps converges at least linearly with a factor ≈ formula_36 .
However the following result of Schönhage yields locally quadratic convergence. To this end let "S" have "m" distinct eigenvalues formula_37 with multiplicities formula_38 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 formula_40 denotes the result then
Thus convergence becomes quadratic as soon as 
formula_42
Cost.
Each Jacobi rotation can be done in O("n") steps when the pivot element "p" is known. However the search for "p" requires inspection of all "N" ≈  "n"2 off-diagonal elements. We can reduce this to O("n") complexity too if we introduce an additional index array formula_43 with the property that formula_44 is the index of the largest element in row "i", ("i" = 1, ..., "n" − 1) of the current "S". Then the indices of the pivot ("k", "l") must be one of the pairs formula_45. Also the updating of the index array can be done in O("n") average-case complexity: First, the maximum entry in the updated rows "k" and "l" can be found in O("n") steps. In the other rows "i", only the entries in columns "k" and "l" change. Looping over these rows, if formula_46 is neither "k" nor "l", it suffices to compare the old maximum at formula_46 to the new entries and update formula_46 if necessary. If formula_46 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("n") complexity. However, this will happen on average only once per rotation. Thus, each rotation has O("n") and one sweep O("n"3) average-case complexity, which is equivalent to one matrix multiplication. Additionally the formula_44 must be initialized before the process starts, which can be done in "n"2 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 formula_51.
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, i.e. formula_52 is an eigenvalue and the column formula_53 an orthonormal eigenvector for formula_52, "i" = 1, ..., "n".
 procedure jacobi("S" ∈ R"n"×"n"; out "e" ∈ R"n"; out "E" ∈ R"n"×"n")
 var
 "i", "k", "l", "m", "state" ∈ N
 "s", "c", "t", "p", "y", "d", "r" ∈ R
 "ind" ∈ N"n"
 "changed" ∈ L"n"
 function maxind("k" ∈ N) ∈ 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("k" ∈ N; "t" ∈ R) ! "update ek and its status"
 "y" := "e""k"; "e""k" := "y"+"t"
 if "changed""k" and ("y"="e""k") then "changed""k" := false; "state" := "state"−1
 elsif (not "changed""k") and ("y"≠"e""k") then "changed""k" := true; "state" := "state"+1
 endif
 endproc
 procedure rotate("k","l","i","j" ∈ N) ! "perform rotation of Sij, Skl"
 │"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("k"); "e""k" := "S""kk"; "changed""k" := true endfor
 while "state"≠0 do ! "next rotation"
 "m" := 1 ! "find index (k,l) of pivot p"
 for "k" := 2 to "n"−1 do
 if │"S""k" "ind""k"│ > │"S""m" "ind""m"│ then "m" := "k" endif
 endfor
 "k" := "m"; "l" := "ind""m"; "p" := "S""kl"
 ! "calculate c = cos φ, s = sin φ"
 "y" := ("e""l"−"e""k")/2; "d" := │"y"│+√("p"2+"y"2)
 "r" := √("p"2+"d"2); "c" := "d"/"r"; "s" := "p"/"r"; "t" := "p"2/"d"
 if "y"<0 then "s" := −"s"; "t" := −"t" endif
 "S""kl" := 0.0; update("k",−"t"); update("l","t")
 ! "rotate rows and columns k and l
 for "i" := 1 to "k"−1 do rotate("i","k","i","l") endfor
 for "i" := "k"+1 to "l"−1 do rotate("k","i","i","l") endfor
 for "i" := "l"+1 to "n" do rotate("k","i","l","i") endfor
 ! "rotate eigenvectors"
 for "i" := 1 to "n" do
 │"E""ik"│ │"c" −"s"││"E""ik"│
 │"E""il"│ │"s" "c"││"E""il"│
 endfor
 ! "rows k, l have changed, update rows indk, indl"
 "ind""k" := maxind("k"); "ind""l" := maxind("l")
 loop
 endproc
Notes.
1. The logical array "changed" holds the status of each eigenvalue. If the numerical value of formula_55 or formula_56 changes during an iteration, the corresponding component of "changed" is set to "true", otherwise to "false". The integer "state" counts the number of components of "changed" which have the value "true". Iteration stops as soon as "state" = 0. This means that none of the approximations formula_57 has recently changed its value and thus it is not very likely that this will happen if iteration continues. Here it is assumed that floating point operations are optimally rounded to the nearest floating point number.
2. The upper triangle of the matrix "S" is destroyed while the lower triangle and the diagonal are unchanged. Thus it is possible to restore "S" if necessary according to
 for "k" := 1 to "n"−1 do ! "restore matrix S"
 for "l" := "k"+1 to "n" do
 "S""kl" := "S""lk"
 endfor
 endfor
3. The eigenvalues are not necessarily in descending order. This can be achieved by a simple sorting algorithm.
 for "k" := 1 to "n"−1 do
 "m" := "k"
 for "l" := "k"+1 to "n" do
 if "e""l" > "e""m" then
 "m" := "l"
 endif
 endfor
 if "k" ≠ "m" then
 swap "e""m","e""k"
 swap "E""m","E""k"
 endif
 endfor
4. The algorithm is written using matrix notation (1 based arrays instead of 0 based).
5. When implementing the algorithm, the part specified using matrix notation must be performed simultaneously.
6. This implementation does not correctly account for the case in which one dimension is an independent subspace. For example, if given a diagonal matrix, the above implementation will never terminate, as none of the eigenvalues will change. Hence, in real implementations, extra logic must be added to account for this case.
Example.
Let 
formula_58
Then "jacobi" produces the following eigenvalues and eigenvectors after 3 sweeps (19 iterations) :
formula_59
formula_60
formula_61
formula_62
formula_63
formula_64
formula_65
formula_66
Applications for real symmetric matrices.
When the eigenvalues (and eigenvectors) of a symmetric matrix are known, the following
values are easily calculated.
Generalizations.
The Jacobi Method has been generalized to complex Hermitian matrices, general nonsymmetric real and complex matrices as well as block matrices.
Since singular values of a real matrix are the square roots of the eigenvalues of the symmetric matrix formula_101 it can also be used for the calculation of these values. For this case, the method is modified in such a way that "S" must not be explicitly calculated which reduces the danger of round-off errors. Note that formula_102 with formula_103 .
The Jacobi Method is also well suited for parallelism.