Jacobi eigenvalue algorithm

The Jacobi Eigenvalue - Method for Symmetric Matrices

Description

For $\phi \in R$ and $1\leq k<l\leq n$ let $J(j,k,l)$ denote the matrix with components

$J_{kk}=\cos \phi ,J_{kl}=-\sin \phi ,J_{lk}:=\sin \phi ,J_{ll}=\cos \phi ,\quad J_{ij}=d_{ij}otherwise.$

$J(\phi ,k,l)$ describes a rotation with angle f in the plane spanned by the k. and l. unit vector. We have $J(\phi ,k,l)^{T}=J(\phi ,k,l)^{-1}$, i.e. $J(\phi ,k,l)$ is orthogonal.

For a real symmetric matrix S let $S(\phi ,k,l):=J(\phi ,k,l)^{-1}SJ(\phi ,k,l)$. This matrix differs from S only in rows and columns k and l where ( if $S'=S(\phi ,k,l),c=\cos \phi ,s=\sin \phi$ for short) : Failed to parse (syntax error): {\displaystyle S'_{kj} = S'_{jk} = c · S_{kj} + s · S_{jl} , S'_{lj} = S '_{jl} = - s · S_{kj} + c · S_{jl} \quad for \quad j ¹ k, l } and Failed to parse (syntax error): {\displaystyle S'_{kk} = c^2 S_{kk} - 2 s c S_{kl} + s^2 S_{ll} , S'_{ll} = s^2 S_{kk} - 2 s c S_{kl} + c^2 S_{ll0 , S'_{kl} = S'_{lk} = (c^2 - s^2 ) · S_{kl} + s c · (S_{ll} - S_{kk} ). }

Define $||S||_{F}^{2}:=trace(S^{T}S)$ (the sum of squares of all components), $\sigma ^{2}:=\sum _{i=1}^{n}S_{ii}^{2}$ (the sum of squares of all diagonal components), $\Gamma (S)^{2}:=||S||_{F}^{2}-\sigma ^{2}$ (the sum of squares of all off-diagonal components). $||S||_{F}$ is the Frobenius-Norm of S and "G"(S ) is called off-diag norm since "G"(S ) = 0 if and only if S is diagonal.

The F - Norm does not change under unitary rotations, so we have $||S'||_{F}=||S||_{F}$ . If we choose f according to $\tan 2\phi ={\frac {2S_{kl}}{S_{kk}-S_{ll}}}$ then $S'_{kl}=0$ and one finds $\Gamma (S')^{2}=\Gamma (S)^{2}-2S_{kl}^{2}$. Thus the off-norm decreases and hence S' becomes "more" diagonal than S. In order to optimize this effect, $S_{kl}$ should be the largest off-diagonal component. Such an element is called a pivot and the rotated matrix $S^{J}$ is called a Jacobi - Rotation.

The Jacobi Eigenvalue - Method repeatedly performs Jacobi - Rotations until the matrix becomes almost diagonal. Then the elements in the diagonal are approximations of the (real) eigenvalues of S.

Convergence

If $p=S_{kl}$ is a pivot element, then by definition $|S_{ij}|\leq |p|$ for $1\leq i,j\leq n,i\neq j$ . Since S has exactly N  := ½ n ( n - 1) off-diag elements, we have Failed to parse (syntax error): {\displaystyle p^2 \le \Gamma(S )^2 £ 2 N p^2 } or $2p^{2}\geq \Gamma (S)^{2}/N$ . This implies $\Gamma (S^{J})^{2}\leq (1-1/N)\Gamma (S)^{2}$ or $\Gamma (S^{J})\leq (1-1/N)^{1/2}\Gamma (S)$ , i.e. the sequence of Jacobi - Rotations converges at least linearly by a factor $(1-1/N)^{1/2}$ to a diagonal matrix.

A number of N Jacobi - Rotations is called a sweep; let $S^{\sigma }$ denote the result. The previous estimate yields $\Gamma (S^{\sigma })\leq (1-1/N)^{N/2}\Gamma (S)$, i.e. the sequence of sweeps converges linearly with a factor  » $e^{1/2}$ .

However the following result of Schönhage yields locally quadratic convergence. To this end let S have m distinct eigenvalues Failed to parse (syntax error): {\displaystyle \lambda_1, … , \lambda_m } with multiplicities Failed to parse (syntax error): {\displaystyle \nu_1, … , \nu_m } and let $d\geq 0$ be the smallest distance of two different eigenvalues. Let us call a number of Failed to parse (syntax error): {\displaystyle N_S := \frac^1_2 n (n - 1) - \sum{\mu=1}{m} \frac^1_2 \nu_{\mu} (\nu_{\mu} - 1) \le N } Jacobi - Rotations a Schönhage - sweep. If $S^{s}$ denotes the result then Failed to parse (syntax error): {\displaystyle \Gamma(S^ s ) \le\sqrt{\frac^n_2 - 1} \frac{\gamma^2}{d - 2\gamma}, \quad \gamma := \Gamma(S ) } . Thus convergence becomes quadratic as soon as Failed to parse (syntax error): {\displaystyle \Gamma(S ) < d / (2 + \sqrt{\frac^n_2 - 1}) } .

Cost

Each Jacobi - Rotation can be done in n steps when the pivot element p is known. However the search for p requires inspection of all Failed to parse (syntax error): {\displaystyle N » \frac^1_2 n^2 } off-diag elements. We can reduce this to n steps too if we introduce an additional index array Failed to parse (syntax error): {\displaystyle m_1, … , m_{n - 1} } with the property that $m_{i}$ is the index of the largest element in row i, (i = 1, … , n - 1) of the current S. Then (k, l) must be one of the pairs $(i,m_{i})$ . Since only columns k and l change, only $m_{k}{\text{and}}m_{l}$ must be updated, which again can be done in n steps. Thus each rotation has O (n) cost and one sweep has O (n³) cost which is equivalent to one matrix multiplication. Additionally the $m_{i}$ must be initialized before the process starts, this can be done in n² steps.

Typically the Jacobi - method converges within numerical precision after a small number of sweeps.

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. $e_{i}$ is an eigenvalue and the column $E_{i}$ an orthonormal eigenvector for $ei$ , i = 1, … , n.

procedure jacobi ($S\in IR^{(n,n)}$; out $e\in IR^{n}$; out Failed to parse (unknown function "\im"): {\displaystyle E \im IR^{(n,n)}} )

var $i,k,l,m\in IN$ $s,c,t,p,psum,y\in IR$ $ind\in IN^{n}$

function maxind $(k\in IN)\in IN$ ! index of largest element in row k begin m := k + 1 for i := k + 2 to n do if $|S_{il}|>|S_{im}|$ then m := l endif endfor return m endfunc

procedure rotate ($k,l,i,j\in IN$ ) ! perform rotation of </math>S_{ij} , S_{kl} </math> begin $\left({\begin{array}{c}S_{kl}\\S_{ij}\end{array}}\right)=\left({\begin{array}{cc}c&-s\\s&c\end{array}}\right)\left({\begin{array}{c}S_{kl}\\S_{ij}\end{array}}\right)$ endproc

begin

! init e, E and index array ind E := 1n; psum := 0.0 for k := 1 to n do $ind_{k}:=maxind(k);e_{k}:=S_{kk}$ endfor

repeat ! next rotation

! find index (k, l) of pivot p m := 1 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}''!teststoppingcriteria''''y'':=''psum'';''psum'':=''psum''+|''p''|;'''if'''''y''<math>\geq$ psum then exit endif

! calculate c = cos f, s = sin f y := ( $e_{l}-e_{k}$ ) / 2; t := |y| + sqrt(p ² + y ²); s := sqrt(p ² + t ²); c := t / s; s := p / s; t := p ² / t if y < 0.0 then s := - s; t := - t endif

! update" $S_{kl}$ and e $S_{kl}:=0.0;e_{k}:=e_{k}-t;e_{l}:=e_{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 $\left({\begin{array}{c}E_{ki}\\E_{li}\end{array}}\right)=\left({\begin{array}{cc}c&-s\\s&c\end{array}}\right)\left({\begin{array}{c}E_{ki}\\E_{li}\end{array}}\right)$ endfor

! rows k, l have changed, update $ind_{k},ind_{l}$ $ind_{k}$  := maxind (k); $ind_{l}$  := maxind (l)

loop

endproc

Notes

The absolute values of the pivot are summed in the variable psum. The iteration stops when the current pivot is numerically small compared to psum. It is unlikely that further iterations produce better approximations. The statement psum := psum + |p|; must be executed with rounding to the nearest floating point number.

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 for l := k + 1 to n do $S_{k}l:=S_{l}k$ endfor endfor

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

Example

Let $S:=\left({\begin{array}{cccc}4&,-30&60&-35\\-30&300&-675&420\\60&-675&1620&-1050\\-35&420&-1050&700\end{array}}\right)$

Then jacobi produces the following eigenvalues and eigenvectors after 3 sweeps (19 iterations) : Failed to parse (syntax error): {\displaystyle e_1 = 2585.25381092892231, \\ \left( \begin{array}{cccc} 0.0291933231647860588 & -0.328712055763188997 & 0.791411145833126331 & -0.514552749997152907 \\ \end{array} \right) , \\ e_2 = 37.1014913651276582, \\ \left( \begin{array}{cccc} -0.179186290535454826 & 0.741917790628453435 & -0.100228136947192199 & -0.638282528193614892 \\ \end{array} \right) , \\ e_3 = 1.4780548447781369, \\ \left( \begin{array}{cccc} -0.582075699497237650 & 0.370502185067093058 & 0.509578634501799626 & 0.514048272222164294\\ \end{array} \right) , \\ e_4 = 0.1666428611718905, \\ \left( \begin{array}{cccc} 0.792608291163763585 & 0.451923120901599794 & 0.322416398581824992 & 0.252161169688241933\\ \end{array} \right) , \\ }

References

1. Jacobi, Über ein leichtes Verfahren, die in der Theorie der Säkularstörungen vorkommenden Gleichungen numerisch aufzulösen, Crelle's Journal 30 (1846), 51 - 94 2. Schönhage, Zur quadratischen Konvergenz Jacobi - Verfahrens, Numer. Math. 6 (1964), 410 - 412 3. Rutishauser, The Jacobi Method for Real Symmetric Matrices, Numer. Math. 9 (1966), 1 - 10