Jacobi rotation

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 .

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} = s^2 a_{kk} + c^2 a_{\ell\ell} + 2 s c a_{k\ell}. \,\!

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 of 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.

References

Golub, Gene H.; Van Loan, Charles F. (1996), Matrix Computations (3rd ed.), Baltimore: Johns Hopkins University Press, ISBN 978-0-8018-5414-9

Numerical linear algebra

Key concepts	Floating point Numerical stability

Problems	Matrix multiplication (algorithms) Matrix decompositions Linear equations Sparse problems

Hardware	CPU cache TLB Cache-oblivious algorithm SIMD Multiprocessing

Software	BLAS Specialized libraries General purpose software

This article is issued from Wikipedia - version of the 9/11/2014. The text is available under the Creative Commons Attribution/Share Alike but additional terms may apply for the media files.

Jacobi rotation

Numerically stable computation

See also

References