Jacobi Transformations of a Symmetric Matrix

The orthogonal transformations P_pq 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 P_pq 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:

$$\boldsymbol{P}_{pq} = \left( \begin{array}{ccccccc} 1 &&&&&... ... & c && \\ &&&&& \cdots & \\ &&&&&& 1 \end{array} \right)$$ (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:

$$\cos{\varphi}$$ c = (3.6)

$$\sin{\varphi}$$ s = (3.7)

hence

c² + s² = 1 (3.8)

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

$$\boldsymbol{A}' = \boldsymbol{P}_{pq}^{T} \cdot \boldsymbol{A} \cdot \boldsymbol{P}_{pq}$$ (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 P_pq in terms of Kronecker deltas:

$$\delta_{ij} + \delta_{ip}\delta_{jp} (c - 1) + \delta_{iq}\delta_{jq} (c - 1) + \delta_{ip}\delta_{jq} s - \delta_{jp}\delta_{iq} s$$ P_ij =

$$\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^-1_ri a_ij P_jp = P_ir a_ij P_jp

$$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)$$ =

$$P_{ir} a_{ij} \left( \delta_{jp} + \delta_{jp} \left( c - 1 \right) - \delta_{jq}\delta_{pp} s \right)$$ =

$$P_{ir} \left( a_{ip} + (c - 1) a_{ip} - s a_{iq} \right) = P_{ir} \left( c a_{ip} - s a_{iq} \right)$$ =

$$\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)$$ =

$$\delta_{ir} \left( c a_{ip} - s a_{iq} \right) = \left( c a_{rp} - s a_{rq} \right)$$ = (3.11)

In summary:

    

$$c a_{rp} - s a_{rq} \qquad \hbox{for} \qquad r\neq p, r\neq q$$ a'_rp = (3.12)

$$c a_{rq} + s a_{rp} \qquad \hbox{for} \qquad r\neq p, r\neq q$$ a'_rq = (3.13)

a'_pp = c² a_pp + s² a_qq - 2 s c a_pq (3.14)

a'_qq = s² a_pp + c² a_qq + 2 s c a_pq (3.15)

$$\left( c^2 - s^2 \right) a_{pq} + s c \left( a_{pp} - a_{qq} \right) = a'_{qp}$$ a'_pq = (3.16)

The purpose of the Jacobi rotation P_pq is to kill a'_pq. Thus:

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

which implies

$$\frac{c^2 - s^2} {s c} = \frac{a_{qq} - a_{pp}}{a_{pq}}$$ (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:

$$\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$$ (3.19)

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

Next recall the following simple algebraic identity:

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

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

$$\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$$ (3.21)

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

$$t^2 + 2 \theta t - 1 = 0$$ (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:

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

Once we have the t, we get c and s as follows^3.1

$$\frac{1}{\sqrt{t^2 + 1}}$$ c = (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 = (c² - s²) a_pq + sc (a_pp - a_qq) (3.26)

by definition. Hence:

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

Then we have

a_pp' = c² a_pp + s² a_qq - 2sca_pq

$$c^2 a_{pp} + s^2 \left( \frac{c^2 - s^2}{sc}a_{pq} + a_{pp} \right) - 2sca_{pq}$$ =

$$(c^2 + s^2) a_{pp} + a_{pq} \left( s^2 \frac{c^2 - s^2}{sc} - 2sc \right)$$ =

$$a_{pp} + a_{pq} \frac{s^2c^2 - s^4 - 2s^2c^2}{sc}$$ =

= a_pp - t a_pq (3.28)

Similarly from the same equation (3.26) we get:

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

and then

a_qq' = s² a_pp + c² a_qq + 2 s c a_pq

$$s^2\left( a_{qq} - \frac{c^2 - s^2}{sc}a_{pq} \right) + c^2 a_{qq} + 2 s c a_{pq}$$ =

$$(s^2 + c^2) a_{qq} + a_{pq} \left( 2 s c - s^2 \frac{c^2 - s^2}{sc} \right)$$ =

$$a_{qq} + a_{pq} \frac{2 s^2 c^2 - s^2 c^2 + s^4}{sc}$$ =

= a_qq + t a_pq (3.30)

For a'_rp and a'_rq the computation is more trivial:

a'_rp = c a_rp - s a_rq

= a_rp + (c - 1) a_rp - s a_rq

$$a_{rp} - s \left(a_{rq} + \frac{1 - c}{s} a_{rp}\right)$$ = (3.31)

and

a'_rq = c a_rq + s a_rp

= a_rq + (c - 1) a_rq + s a_rp

$$a_{rq} + s \left(a_{rp} - \frac{1 - c}{s} a_{rq}\right)$$ = (3.32)

where

$$\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}}$$ (3.33)