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, [1] but only became widely used in the 1950s with the advent of computers. [2]
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 (the square-root sum of squares of all components), 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, 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.
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 [3] 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
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 ≈ 1/2 n2 off-diagonal elements. We can reduce this to O(n) complexity too if we introduce an additional index array with the property that 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 . 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 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(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. Additionally the must be initialized before the process starts, which can be done in n2 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 .
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(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
1. The logical array changed holds the status of each eigenvalue. If the numerical value of or 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 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.
Let
Then jacobi produces the following eigenvalues and eigenvectors after 3 sweeps (19 iterations) :
When the eigenvalues (and eigenvectors) of a symmetric matrix are known, the following values are easily calculated.
The following code is a straight-forward implementation of the mathematical description of the Jacobi eigenvalue algorithm in the Julia programming language.
using LinearAlgebra, Test
function find_pivot(Sprime)
n = size(Sprime,1)
pivot_i = pivot_j = 0
pivot = 0.0
for j = 1:n
for i = 1:(j-1)
if abs(Sprimei,j]) > pivot
pivot_i = i
pivot_j = j
pivot = abs(Sprimei,j])
end
end
end
return (pivot_i, pivot_j, pivot)
end
# in practice one should not instantiate explicitly the Givens rotation matrix
function givens_rotation_matrix(n,i,j,θ)
G = Matrix{Float64}(I,(n,n))
Gi,i = Gj,j = cos(θ)
Gi,j = sin(θ)
Gj,i = -sin(θ)
return G
end
# S is a symmetric n by n matrix
n = 4
sqrtS = randn(n,n);
S = sqrtS*sqrtS';
# the largest allowed off-diagonal element of U' * S * U
# where U are the eigenvectors
tol = 1e-14
Sprime = copy(S)
U = Matrix{Float64}(I,(n,n))
while true
(pivot_i, pivot_j, pivot) = find_pivot(Sprime)
if pivot < tol
break
end
θ = atan(2*Sprimepivot_i,pivot_j/(Sprimepivot_j,pivot_j - Sprimepivot_i,pivot_i )) / 2
G = givens_rotation_matrix(n,pivot_i,pivot_j,θ)
# update Sprime and U
Sprime .= G'*Sprime*G
U .= U * G
end
# Sprime is now (almost) a diagonal matrix
# extract eigenvalues
λ = diag(Sprime)
# sort eigenvalues (and corresponding eigenvectors U) by increasing values
i = sortperm(λ)
λ = λi
U = U:,i
# S should be equal to U * diagm(λ) * U'
@test S ≈ U * diagm(λ) * U'
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 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 with .
The Jacobi Method is also well suited for parallelism.
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, [1] but only became widely used in the 1950s with the advent of computers. [2]
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 (the square-root sum of squares of all components), 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, 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.
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 [3] 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
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 ≈ 1/2 n2 off-diagonal elements. We can reduce this to O(n) complexity too if we introduce an additional index array with the property that 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 . 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 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(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. Additionally the must be initialized before the process starts, which can be done in n2 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 .
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(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
1. The logical array changed holds the status of each eigenvalue. If the numerical value of or 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 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.
Let
Then jacobi produces the following eigenvalues and eigenvectors after 3 sweeps (19 iterations) :
When the eigenvalues (and eigenvectors) of a symmetric matrix are known, the following values are easily calculated.
The following code is a straight-forward implementation of the mathematical description of the Jacobi eigenvalue algorithm in the Julia programming language.
using LinearAlgebra, Test
function find_pivot(Sprime)
n = size(Sprime,1)
pivot_i = pivot_j = 0
pivot = 0.0
for j = 1:n
for i = 1:(j-1)
if abs(Sprimei,j]) > pivot
pivot_i = i
pivot_j = j
pivot = abs(Sprimei,j])
end
end
end
return (pivot_i, pivot_j, pivot)
end
# in practice one should not instantiate explicitly the Givens rotation matrix
function givens_rotation_matrix(n,i,j,θ)
G = Matrix{Float64}(I,(n,n))
Gi,i = Gj,j = cos(θ)
Gi,j = sin(θ)
Gj,i = -sin(θ)
return G
end
# S is a symmetric n by n matrix
n = 4
sqrtS = randn(n,n);
S = sqrtS*sqrtS';
# the largest allowed off-diagonal element of U' * S * U
# where U are the eigenvectors
tol = 1e-14
Sprime = copy(S)
U = Matrix{Float64}(I,(n,n))
while true
(pivot_i, pivot_j, pivot) = find_pivot(Sprime)
if pivot < tol
break
end
θ = atan(2*Sprimepivot_i,pivot_j/(Sprimepivot_j,pivot_j - Sprimepivot_i,pivot_i )) / 2
G = givens_rotation_matrix(n,pivot_i,pivot_j,θ)
# update Sprime and U
Sprime .= G'*Sprime*G
U .= U * G
end
# Sprime is now (almost) a diagonal matrix
# extract eigenvalues
λ = diag(Sprime)
# sort eigenvalues (and corresponding eigenvectors U) by increasing values
i = sortperm(λ)
λ = λi
U = U:,i
# S should be equal to U * diagm(λ) * U'
@test S ≈ U * diagm(λ) * U'
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 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 with .
The Jacobi Method is also well suited for parallelism.