Jacobi rotation

In numerical linear algebra, a Jacobi rotation is a rotation, Qkℓ, 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 Qkℓ 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, Qkℓ is an identity matrix except for four entries, two on the diagonal (qkk and qℓℓ, both equal to c) and two symmetrically placed off the diagonal (qkℓ 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. 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 akℓ 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{1-c}{s}, $$

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 Qkℓ by retaining only the three values k, ℓ, and t, with t set to zero for a null rotation.

Tridiagonal example
Some applications may require multiple zero entries in a similarity matrix, possibly in the form of a tridiagonal matrix. Since Jacobian rotations may remove zeros from other cells that were previously zeroed, it is usually not possible to achieve tridiagonalization by simply zeroing each off-tridiagonal cell individually in a medium to large matrix. However, if Jacobian rotations are repeatedly performed on the above-tridiagonal cell with the highest absolute value using an adjacent cell just below or to the left to rotate on, then the all of the off-triangular cells are expected to converge on zero after several iterations. In the example below, $$A$$ is a 5X5 matrix that is to be tridiagonalized into a similar matrix, $$T$$.

$$\begin{align} &A= \begin{bmatrix} 2 & 5 &9 &11 &7\\ 5 & 3 &6 &-2 &5\\ 9 & 6 &7 &3 &1\\ 11 & -2 &3 &1 &3\\ 7 & 5 &1 &3 &4\\ \end{bmatrix} \\

&\text{eigenvalues of }A = \begin{bmatrix}4.1714549     &5.9308002      &-5.3682585      &24.054762      &-11.788758 \end{bmatrix} \\ &|A| = 37662 \end{align}$$

To tridiagonalize matrix $$A$$ into matrix $$T$$, the off-tridiagonal cells [1,3], [1,4], [1,5], [2,4], [2,5], and [3,5], must continue to be iteratively zeroed until the maximum absolute value of those cells is below an acceptable convergence threshold. This example will use 1.e-14.. The cells below the diagonal will be zeroed automatically, due to the symmetric nature of the matrix. The first Jacobian rotation will be on the off-tridiagonal cell with the highest absolute value, which by inspection is [1,4] with a value of 11. To make this entry zero, the condition specified in the above equations must be met for the cell coordinates to be zeroed ($$h=1, l=4$$) and for the selected rotational coordinates of $$S$$ ($$h=1, k=3$$), and are reproduced below for the first iteration.

$$\begin{align} \\ &9 \sin(\theta) + 11\cos(\theta) = 0 \\ & \theta = \tan^{-1}(-11/9) = -0.8850668 \\ & \sin(\theta) = -0.77395730 \\ & \cos(\theta) = 0.66323890 \\ & \\ &\text{More conveniountly:} \\ &r = \sqrt(11^2+9^2) = 14.2126704 \\ &\sin(\theta) = -11/r = -0.77395730 \\ &\cos(\theta) = 9/r = 0.66323890 \\ &\\ &S= \begin{bmatrix} &1 &0 &0 &0 &0 \\ &0 &0.66323890 &0 &-0.77395730 &0  \\ &0 &0 &1 &0 &0  \\ &0 &0.77395730 &0 &0.66323890 &0  \\ &0 &0 &0 &0 &1  \\ \end{bmatrix}  \\ &\\ &T_1=S^TAS = \begin{bmatrix} 2     &12.083046      &9      &0      &7\\ 12.083046      &-0.16438356      &5.2139171      &0.56164384      &4.8001142\\ 9      &5.2139171      &7      &-4.22079      &1\\ 0      &0.56164384      &-4.22079      &4.1643836      &-3.3104236\\ 7      &4.8001142      &1      &-3.3104236      &4\\ \end{bmatrix}  \\ &\text{eigenvalues of }T_1 = \begin{bmatrix}4.1714549     &5.9308002      &-5.3682585      &24.054762      &-11.788758 \end{bmatrix} \\ &|T_1| = 37662 \end{align}$$

The first rotation iteration, $$T_1$$, produces a matrix with cells [1,4] and [4,1] zeroed, as expected. Furthermore, the eigenvalues and determinant of $$T_1$$are identical to those of $$A$$ and T1 is also symmetric, confirming that the Jacobian rotation was performed correctly. The next iteration for $$T_2$$ will select cell [2,5] which contains the highest absolute value, 4.8001142, of all the cells to be zeroed..

After 10 iterations of zeroing the cell with the maximum absolute value using Jacobian rotations on the cell just below it, the maximum absolute value of all off-tridiagonal cells is 2.6e-15. Assuming this convergence criteria is acceptably low for the application it is being performed for, the similar triangularized $$T$$ matrix is shown below.

$$\begin{align} &T= \begin{bmatrix} 2 & 16.613248     &0      &0      &0\\ 16.613248      &10.184783      &5.1109346      &0      &0\\ 0      &5.1109346      &4.0304188      &4.5541372      &0\\ 0      &0       &4.5541372       &-3.3901672       &0.43515335\\ 0      &0      &0      &0.43515335      &4.1749658\\ \end{bmatrix}  \\

&\text{eigenvalues of }T = \begin{bmatrix}4.1714549     &5.9308002      &-5.3682585      &24.054762      &-11.788758 \end{bmatrix} \\ &|T| = 37662 \end{align}$$

Since $$A$$ and $$T$$ have identical eigenvalues and determinants and $$T$$ is also symmetric, $$A$$ and $$T$$ are similar matrices with $$T$$ being tridiagonalized.

Eigenvalues example
Jacobian rotation can be used to extract the eigenvalues in a similar manner as the triangulation example above, but by zeroing all of the cells above the diagonal, instead of the tridiagonal, and performing the Jacobian rotation directly in the cells to be zeroed, instead of an adjacent cell.

Starting with the same matrix $$A$$ as the tridiagonal example,

$$\begin{align} &A= \begin{bmatrix} 2 & 5 &9 &11 &7\\ 5 & 3 &6 &-2 &5\\ 9 & 6 &7 &3 &1\\ 11 & -2 &3 &1 &3\\ 7 & 5 &1 &3 &4\\ \end{bmatrix} \\

&\text{eigenvalues of }A = \begin{bmatrix}4.1714549     &5.9308002      &-5.3682585      &24.054762      &-11.788758 \end{bmatrix} \\ &|A| = 37662 \end{align}$$

The first Jacobian rotation will be on the off-diagonal cell with the with the highest absolute value, which by inspection is [1,4] with a value of 11, and the rotation cell will also be [1,4], $$l = 4, k = 1$$ in the equations above. The rotation angle is the result of a quadratic solution, but it can be seen in the equation that if the matrix is symmetric, then a real solution is assured.

$$\begin{align} \\ &11(cos^2(\theta) - sin^2(\theta)) + (2-1)cos(\theta)sin(\theta) = 0 \\ &\text{dividing by }cos^2(\theta) \\ &tan^2(\theta) - tan(\theta)/11 - 1 = 0 \\ &tan(\theta) = [1/11 + \sqrt{1/11^2 + 4}]/2 = 1.0464871 \\ &\theta = tan^{-1}(1.0464871) = -0.80810980 \\ &sin(\theta) = 0.72298259 \\ &cos(\theta) = 0.69086625 \\ & \\ &\text{More conveniountly:} \\ &r = sqrt(1.0464871^2+1) = 1.4474582 \\ &sin(\theta) = 1.0464871/r = 0.72298259 \\ &cos(\theta) = 1/r = 0.69086625 \\ &\\ &S= \begin{bmatrix} &0.69086625 &0 &0 &0.72298259 &0 \\ &0 &1 &0 &0 &0  \\ &0 &0 &1 &0 &0  \\ &-0.72298259 &0 &0 &0.69086625 &0  \\ &0 &0 &0 &0 &1  \\ \end{bmatrix}  \\ &\\ &T_1=S^TAS = \begin{bmatrix} -9.5113578     &4.9002964      &4.0488484      &0      &2.6671159\\ 4.9002964      &3      &6      &2.2331805      &5\\ 4.0488484      &6      &7      &8.5794421      &1\\ 0      &2.2331805      &8.5794421      &12.511358      &7.1334769\\ 2.6671159      &5      &1      &7.1334769      &4\\ \end{bmatrix}  \\ &\text{eigenvalues of }T_1 = \begin{bmatrix}4.1714549     &5.9308002      &-5.3682585      &24.054762      &-11.788758 \end{bmatrix} \\ &|T_1| = 37662 \end{align}$$

The first rotation iteration, $$T_1$$, produces a matrix with cells [1,4] and [4,1] zeroed, as expected. Furthermore, the eigenvalues and determinant of $$T_1$$are identical to those of $$A$$ and T1 is also symmetric, confirming that the Jacobian rotation was performed correctly. The next iteration for $$T_2$$ will select cell [3,4] which contains the highest absolute value, 8.5794421, of all the cells to be zeroed..

After 25 iterations of zeroing the cell with the maximum absolute value using Jacobian rotations on the cell just below it, the maximum absolute value of all off-diagonal cells is 9.0233029E-11. Assuming this convergence criteria is acceptably low for the application it is being performed for, the similar diagonalized $$T$$ matrix is shown below.

$$\begin{align} &T= \begin{bmatrix} 24.054762  &0      &0      &0      &0\\ 0      &5.9308002      &0      &0      &0\\ 0      &0      &-11.788758      &0      &0\\ 0      &0       &0       &4.1714549       &0\\ 0      &0      &0      &0      &-5.3682585\\ \end{bmatrix}  \\

&\text{eigenvalues of }T = \begin{bmatrix}4.1714549     &5.9308002      &-5.3682585      &24.054762      &-11.788758 \end{bmatrix} \\ &|T| = 37662 \end{align}$$

The eigenvalues are now displayed across the diagonal, and may be directly extracted for use elsewhere.

Since $$A$$ and $$T$$ have identical eigenvalues and determinants and $$T$$ is also symmetric, $$A$$ and $$T$$ are similar matrices with $$T$$ being successfully diagonalized.