Jacobi eigenvalue algorithm on arvau

Wikipedia / 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 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

be a Givens rotation matrix. Then:

is symmetric and similar.

has entries:

have the same Frobenius norm has a larger sum of squares on the diagonal:

Set this equal to 0, and rearrange:

In order to optimize this effect, Sij 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

is a pivot element to a diagonal matrix.

denote the result. The previous estimate yields

;

However the following result of Schönhage yields locally quadratic convergence. To this end let S have m and let d > 0 be the smallest distance of two different eigenvalues. Let us call a number of

denotes the result then

Cost

steps when the pivot element p is known. However the search for p requires inspection of all N≈n2 complexity for a full diagonalisation.

Caching row maximums

We can reduce the complexity of finding the pivot element from O(N) to O(n is the index of the largest element in row i, (i = 1, ..., n1) of the current S. Then the indices of the pivot (k, l. 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 is neither k nor l 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(n3) average-case complexity, which is equivalent to one matrix multiplication must be initialized before the process starts, which can be done in n2 steps.

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. 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 matchingss, 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 processors are available.

for its rotation, and then applying the rotation from the left (rotating between rows). Next, the processors synchronise before applying the transpose rotation from the right (rotating between columns), 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, i = 1, ..., n.

procedure jacobi(S ∈ Rn×n; out e ∈ Rn; out E ∈ Rn×n)
  var
    i, k, l, m, state ∈ N
    s, c, t, p, y, d, r ∈ R
    ind ∈ Nn
    changed ∈ Ln

  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 │Ski│ > │Skm│ then m := i endif
    endfor
    return m
  endfunc

  procedure update(k ∈ N; t ∈ R) ! update ek and its status
    y := ek; ek := y+t
    if changedk and (y=ek) then changedk := false; state := state−1
    elsif (not changedk) and (y≠ek) then changedk := true; state := state+1
    endif
  endproc

  procedure rotate(k,l,i,j ∈ N) ! perform rotation of Sij, Skl
    ┌   ┐    ┌     ┐┌   ┐
    │Skl│    │c  −s││Skl│
    │   │ := │     ││   │
    │Sij│    │s   c││Sij│
    └   ┘    └     ┘└   ┘
  endproc

  ! init e, E, and arrays ind, changed
  E := I; state := n
  for k := 1 to n do indk := maxind(k); ek := Skk; changedk := 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 │Sk indk│ > │Sm indm│ then m := k endif
    endfor
    k := m; l := indm; p := Skl
    ! calculate c = cos φ, s = sin φ
    y := (el−ek)/2; d := │y│+√(p2+y2)
    r := √(p2+d2); c := d/r; s := p/r; t := p2/d
    if y<0 then s := −s; t := −t endif
    Skl := 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
      ┌   ┐    ┌     ┐┌   ┐
      │Eik│    │c  −s││Eik│
      │   │ := │     ││   │
      │Eil│    │s   c││Eil│
      └   ┘    └     ┘└   ┘
    endfor
    ! update all potentially changed indi
    for i := 1 to n do indi := maxind(i) endfor
  loop
endproc

Notes

1. The logical array changed 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 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
        Skl := Slk
    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 el > em then
            m := l
        endif
    endfor
    if k ≠ m then
        swap em,ek
        swap Em,Ek
    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

Then jacobi produces the following eigenvalues and eigenvectors after 3 sweeps (19 iterations) :

Applications for real symmetric matrices

When the eigenvalues (and eigenvectors) of a symmetric matrix are known, the following values are easily calculated.

Singular values
.
2-norm and spectral radius
The 2-norm of a matrix A. In case of a symmetric matrix it is the largest absolute value of its eigenvectors and thus equal to its spectral radius.
Condition number
The condition number. In case of a symmetric matrix it is the absolute value of the quotient of the largest and smallest eigenvalue. Matrices with large condition numbers can cause numerically unstable results: small perturbation can result in large errors. Hilbert matrices are the most famous ill-conditioned matrices. For example, the fourth-order Hilbert matrix has a condition of 15514, while for order 8 it is 2.7×108.
Rank
is the dimension of the range of. Furthermore it is the number of nonzero singular values.
In case of a symmetric matrix r is the number of nonzero eigenvalues. Unfortunately because of rounding errors numerical approximations of zero eigenvalues may not be zero (it may also happen that a numerical approximation is zero while the true value is not). Thus one can only calculate the numerical rank by making a decision which of the eigenvalues are close enough to zero.
Pseudo-inverse
.
holds where Diag(e) denotes the diagonal matrix with vector e is (numerically close to) zero. Since matrix E.
Least squares solution
. In case of a symmetric matrix S.
Matrix exponential
where exp.
Linear differential equations
.
are called stable and unstable.

Julia implementation

The following code is a straight-forward implementation of the mathematical description of the Jacobi eigenvalue algorithm in the Julia programming language.

Generalizations

The Jacobi Method has been generalized to complex Hermitian matrices, general nonsymmetric real and complex matrices as well as block matrices.

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 .