Jacobi rotation

From Wikipedia, the free encyclopedia

In numerical linear algebra, a Jacobi rotation is a rotation, Q_kℓ, of a 2-dimensional linear subspace of an n-dimensional inner product space, chosen to zero a symmetric pair of off-diagonal entries of an n×n real symmetric matrix, A, when applied as a similarity transformation:

$A \mapsto Q_{k\ell}^T A Q_{k\ell} = A' . \,\!$

$\begin{bmatrix} {*} & & & \cdots & & & * \\ & \ddots & & & & & \\ & & a_{kk} & \cdots & a_{k\ell} & & \\ \vdots & & \vdots & \ddots & \vdots & & \vdots \\ & & a_{\ell k} & \cdots & a_{\ell\ell} & & \\ & & & & & \ddots & \\ {*} & & & \cdots & & & * \end{bmatrix} \to \begin{bmatrix} {*} & & & \cdots & & & * \\ & \ddots & & & & & \\ & & a'_{kk} & \cdots & 0 & & \\ \vdots & & \vdots & \ddots & \vdots & & \vdots \\ & & 0 & \cdots & a'_{\ell\ell} & & \\ & & & & & \ddots & \\ {*} & & & \cdots & & & * \end{bmatrix}$

It is the core operation in the Jacobi eigenvalue algorithm, which is numerically stable and well-suited to implementation on parallel processors.

It is important to note that only rows k and ℓ and columns k and ℓ of A will be affected, and that A′ will remain symmetric. Also, an explicit matrix for Q_kℓ is rarely computed; instead, auxiliary values are computed and A is updated in an efficient and numerically stable way. However, for reference, we may write the matrix as

$Q_{k\ell} = \begin{bmatrix} 1 & & & & & & \\ & \ddots & & & & 0 & \\ & & c & \cdots & s & & \\ & & \vdots & \ddots & \vdots & & \\ & & -s & \cdots & c & & \\ & 0 & & & & \ddots & \\ & & & & & & 1 \end{bmatrix} .$

That is, Q_kℓ is an identity matrix except for four entries, two on the diagonal (q_kk and q_ℓℓ, both equal to c) and two symmetrically placed off the diagonal (q_kℓ and Q_ℓk, equal to s and −s, respectively). Here c = cos ϑ and s = sin ϑ for some angle ϑ; but to apply the rotation, the angle itself is not required. Using Kronecker delta notation, the matrix entries can be written

$q_{ij} = \delta_{ij} + (\delta_{ik}\delta_{jk} + \delta_{i\ell}\delta_{j\ell})(c-1) + (\delta_{ik}\delta_{j\ell} - \delta_{i\ell}\delta_{jk})s . \,\!$

Suppose h is an index other than k or ℓ (which must themselves be distinct). Then the similarity update produces, algebraically,

$a'_{hk} = a'_{kh} = c a_{hk} - s a_{h\ell} \,\!$

$a'_{h\ell} = a'_{\ell h} = c a_{h\ell} + s a_{hk} \,\!$

$a'_{k\ell} = a'_{\ell k} = (c^2-s^2)a_{k\ell} + sc (a_{kk} - a_{\ell\ell}) = 0 \,\!$

$a'_{kk} = c^2 a_{kk} + s^2 a_{\ell\ell} - 2 s c a_{k\ell} \,\!$

$a'_{\ell\ell} = c^2 a_{kk} + s^2 a_{\ell\ell} + 2 s c a_{k\ell} \,\!$

[edit] Numerically stable computation

To determine the quantities needed for the update, we must solve the off-diagonal equation for zero (Golub & Van Loan 1996, §8.4). This implies that

$\frac{c^2-s^2}{sc} = \frac{a_{\ell\ell} - a_{kk}}{a_{k\ell}} .$

Set β to half this quantity,

$\beta = \frac{a_{\ell\ell} - a_{kk}}{2 a_{k\ell}} .$

If a_kℓ is zero we can stop without performing an update, thus we never divide by zero. Let t be tan ϑ. Then with a few trigonometric identities we reduce the equation to

$t^2 + 2\beta t - 1 = 0 . \,\!$

For stability we choose the solution

$t = \frac{\sgn(\beta)}{|\beta|+\sqrt{\beta^2+1}} .$

From this we may obtain c and s as

$c = \frac{1}{\sqrt{t^2+1}} \,\!$

$s = c t \,\!$

Although we now could use the algebraic update equations given previously, it may be preferable to rewrite them. Let

$\rho= \frac{s}{1+c} ,$

so that ρ = tan(ϑ/2). Then the revised update equations are

$a'_{hk} = a'_{kh} = a_{hk} - s (a_{h\ell} + \rho a_{hk}) \,\!$

$a'_{h\ell} = a'_{\ell h} = a_{h\ell} + s (a_{hk} - \rho a_{h\ell}) \,\!$

$a'_{k\ell} = a'_{\ell k} = 0 \,\!$

$a'_{kk} = a_{kk} - t a_{k \ell} \,\!$

$a'_{\ell\ell} = a_{\ell\ell} + t a_{k \ell} \,\!$

As previously remarked, we need never explicitly compute the rotation angle ϑ. In fact, we can reproduce the symmetric update determined by Q_kℓ by retaining only the three values k, ℓ, and t, with t set to zero for a null rotation.