next up previous index
Next: Convergence of the Jacobi Up: Eigensystems Previous: Introduction

   
Jacobi Transformations of a Symmetric Matrix

The orthogonal transformations Ppq annihilate element (p,q) of an object matrix.

Successive transformations undo previously set zeros, but the off-diagonal terms eventually get smaller and smaller, until they get nearly zero, and the only thing that's left is the diagonal with eigenvalues.

Taking the product of all the transformations Ppq yields the matrix of eigenvectors.

The method is absolutely foolproof for symmetric real matrices. But it is slow. Painfully slow for matrices larger than $10\times 10$.

It is a simple and stable algorithm though, and it parallelises well too.

The Jacobi rotation matrix has the form:

\begin{displaymath}\boldsymbol{P}_{pq} = \left(
\begin{array}{ccccccc}
1 &&&&&...
... & c && \\
&&&&& \cdots & \\
&&&&&& 1
\end{array} \right)
\end{displaymath} (3.5)

All diagonal elements are 1 with the exception of (p,p) and (q,q), which are c. The (p,q) element is s and the (q,p) element is -s. All other elements are 0. Furthermore:
c = $\displaystyle \cos{\varphi}$ (3.6)
s = $\displaystyle \sin{\varphi}$ (3.7)

hence

c2 + s2 = 1 (3.8)

Because $\boldsymbol{P}_{pq} \in \hbox{SO}(n)$, Ppq-1 = PpqT and so

 \begin{displaymath}
\boldsymbol{A}' = \boldsymbol{P}_{pq}^{T} \cdot \boldsymbol{A}
\cdot \boldsymbol{P}_{pq}
\end{displaymath} (3.9)

This operation will affect only rows p and q and columns p and qleaving the rest of the matrix unchanged.

It is quite easy to see what the effect of equation (3.9) is going to be on selected terms.

First we need to come up with an expression that describes a generic term of matrix Ppq in terms of Kronecker deltas:

Pij = $\displaystyle \delta_{ij} + \delta_{ip}\delta_{jp} (c - 1)
+ \delta_{iq}\delta_{jq} (c - 1) + \delta_{ip}\delta_{jq} s
- \delta_{jp}\delta_{iq} s$  
  = $\displaystyle \delta_{ij}
+ \left(
\delta_{ip}\delta_{jp} + \delta_{iq}\delta_{...
... - 1 \right)
+ \left(
\delta_{ip}\delta_{jq} - \delta_{iq}\delta_{jp}
\right) s$ (3.10)

Now we can evaluate a'rp for, say, $r \neq p$ and $r \neq q$(assume summation over dummy indexes i and j)

a'rp = P-1ri aij Pjp = Pir aij Pjp  
  = $\displaystyle P_{ir} a_{ij}
\left(
\delta_{jp}
+ \left(
\delta_{jp}\delta_{pp} ...
...ght)
+ \left(
\delta_{jp}\delta_{pq} - \delta_{jq}\delta_{pp}
\right) s
\right)$  
  = $\displaystyle P_{ir} a_{ij}
\left(
\delta_{jp}
+ \delta_{jp} \left( c - 1 \right)
- \delta_{jq}\delta_{pp} s
\right)$  
  = $\displaystyle P_{ir} \left( a_{ip} + (c - 1) a_{ip} - s a_{iq} \right)
= P_{ir} \left( c a_{ip} - s a_{iq} \right)$  
  = $\displaystyle \left(
\delta_{ir}
+ \left(
\delta_{ip}\delta_{rp} + \delta_{iq}\...
...} - \delta_{iq}\delta_{rp}
\right) s
\right)
\left( c a_{ip} - s a_{iq} \right)$  
  = $\displaystyle \delta_{ir} \left( c a_{ip} - s a_{iq} \right)
= \left( c a_{rp} - s a_{rq} \right)$ (3.11)

In summary:

    
a'rp = $\displaystyle c a_{rp} - s a_{rq} \qquad \hbox{for} \qquad r\neq p, r\neq q$ (3.12)
a'rq = $\displaystyle c a_{rq} + s a_{rp} \qquad \hbox{for} \qquad r\neq p, r\neq q$ (3.13)
a'pp = c2 app + s2 aqq - 2 s c apq (3.14)
a'qq = s2 app + c2 aqq + 2 s c apq (3.15)
a'pq = $\displaystyle \left( c^2 - s^2 \right) a_{pq}
+ s c \left( a_{pp} - a_{qq} \right) = a'_{qp}$ (3.16)

The purpose of the Jacobi rotation Ppq is to kill a'pq. Thus:

\begin{displaymath}a'_{pq} = 0 = \left( c^2 - s^2 \right) a_{pq}
+ s c \left( a_{pp} - a_{qq} \right)
\end{displaymath} (3.17)

which implies

 \begin{displaymath}
\frac{c^2 - s^2}
{s c}
= \frac{a_{qq} - a_{pp}}{a_{pq}}
\end{displaymath} (3.18)

Now, observe that:

\begin{eqnarray*}\cos{2\varphi} &=& \cos^2{\varphi} - \sin^2{\varphi} \\
\sin{2\varphi} &=& 2\sin{\varphi}\cos{\varphi}
\end{eqnarray*}


This means that if we divide both sides of equation (3.18) by 2, we'll get:

\begin{displaymath}\cot{2\varphi} = \frac{\cos{2\varphi}}{\sin{2\varphi}}
= \fr...
...^2}
{2 s c}
= \frac{a_{qq} - a_{pp}}{2 a_{pq}} \doteq \theta
\end{displaymath} (3.19)

Let us call this $\cot{2\varphi}$ a $\theta$ for convenience.

Next recall the following simple algebraic identity:

 \begin{displaymath}
\tan^2{\varphi} + 2\tan{\varphi}\cot{2\varphi} - 1 = 0
\end{displaymath} (3.20)

This is easy to see. Since $\tan{\varphi} = s/c$, we have:

\begin{displaymath}\tan^2{\varphi} + 2\tan{\varphi}\cot{2\varphi} - 1 =
\frac{s^...
... s^2}{2 s c} - 1
= \frac{s^2 + c^2 - s^2}{c^2} - 1 = 1 - 1 = 0
\end{displaymath} (3.21)

Now, denote $\tan{\varphi}$ by t, so that the equation looks as follows:

 \begin{displaymath}
t^2 + 2 \theta t - 1 = 0
\end{displaymath} (3.22)

Solving equation (3.22) with respect to t yields the following:

\begin{eqnarray*}\Delta &=& 4\theta^2 - 4\cdot 1\cdot(-1) = 4\theta^2 + 4 \\
t_...
... \pm 2\sqrt{\theta^2 + 1}}{2}
= -\theta \pm \sqrt{\theta^2 + 1}
\end{eqnarray*}


Now, this solution can be rewritten in a computationally more convenient form. Consider the + case:

\begin{eqnarray*}t_+ &=& -\theta + \sqrt{\theta^2 + 1} \\
&=& -1 \frac{\left(\...
...t{\theta^2 + 1}} \\
&=& \frac{1}{\theta + \sqrt{\theta^2 + 1}}
\end{eqnarray*}


The - case yields:

\begin{eqnarray*}t_- &=& -\theta - \sqrt{\theta^2 + 1} \\
&=& -1 \frac{\left(\...
...\theta^2 + 1}} \\
&=& \frac{-1}{-\theta + \sqrt{\theta^2 + 1}}
\end{eqnarray*}


If we get a positive $\theta$ in the + case and a negative $\theta$ in the - case we'll end up with the same smaller t that corresponds to an angle $\varphi < \pi/4$. This will yield the most stable reduction. So we can rewrite the formula for t as follows:

\begin{displaymath}t = \frac{\hbox{sign}(\theta)}{\vert\theta\vert + \sqrt{\theta^2 + 1}}
\end{displaymath} (3.23)

Once we have the t, we get c and s as follows3.1

c = $\displaystyle \frac{1}{\sqrt{t^2 + 1}}$ (3.24)
s = c t (3.25)

Equations (3.12) through (3.16) are now rewritten to minimize the round off error and to make them look like the new quantity is equal to the old one plus a small correction. And so,

 
a'pq = a'qp = 0 = (c2 - s2) apq + sc (app - aqq) (3.26)

by definition. Hence:

\begin{displaymath}a_{qq} = \frac{c^2 - s^2}{sc}a_{pq} + a_{pp}
\end{displaymath} (3.27)

Then we have
app' = c2 app + s2 aqq - 2scapq  
  = $\displaystyle c^2 a_{pp} + s^2 \left( \frac{c^2 - s^2}{sc}a_{pq} + a_{pp} \right)
- 2sca_{pq}$  
  = $\displaystyle (c^2 + s^2) a_{pp} + a_{pq}
\left(
s^2 \frac{c^2 - s^2}{sc} - 2sc
\right)$  
  = $\displaystyle a_{pp} + a_{pq} \frac{s^2c^2 - s^4 - 2s^2c^2}{sc}$  
  = app - t apq (3.28)

Similarly from the same equation (3.26) we get:

\begin{displaymath}a_{pp} = a_{qq} - \frac{c^2 - s^2}{sc}a_{pq}
\end{displaymath} (3.29)

and then
aqq' = s2 app + c2 aqq + 2 s c apq  
  = $\displaystyle s^2\left(
a_{qq} - \frac{c^2 - s^2}{sc}a_{pq}
\right)
+ c^2 a_{qq} + 2 s c a_{pq}$  
  = $\displaystyle (s^2 + c^2) a_{qq} + a_{pq}
\left(
2 s c - s^2 \frac{c^2 - s^2}{sc}
\right)$  
  = $\displaystyle a_{qq} + a_{pq} \frac{2 s^2 c^2 - s^2 c^2 + s^4}{sc}$  
  = aqq + t apq (3.30)

For a'rp and a'rq the computation is more trivial:

a'rp = c arp - s arq  
  = arp + (c - 1) arp - s arq  
  = $\displaystyle a_{rp} - s \left(a_{rq} + \frac{1 - c}{s} a_{rp}\right)$ (3.31)

and
a'rq = c arq + s arp  
  = arq + (c - 1) arq + s arp  
  = $\displaystyle a_{rq} + s \left(a_{rp} - \frac{1 - c}{s} a_{rq}\right)$ (3.32)

where

\begin{displaymath}\frac{1 - c}{s} = \frac{(1-c)(1+c)}{s(1+c)} =
\frac{1 - c^2}{...
...^2}{s(1+c)} =
\frac{s}{1 + c} = \tau = \tan{\frac{\varphi}{2}}
\end{displaymath} (3.33)



 
next up previous index
Next: Convergence of the Jacobi Up: Eigensystems Previous: Introduction
Zdzislaw Meglicki
2001-02-26