QR algorithm: finding eigenvalues

Motivation

A lot of iterative algorithms exist in linear algebra to find eigenvalues of symmetric matrices. QR algorithm addresses the issue of finding eigenvalues of nonsymmetric matrices as well. More accurately this algorithm can find real eigenvalues but complex ones are not treated well.

How does it work?

Assume we have a matrix A with real eigenvalues. Then it is possible to find an orthogonal Q and an upper-triangular R matrix such that:

\begin{equation} A = Q \cdot R. \end{equation}

For clarity:

  • orthogonality means: $Q \cdot Q^T = E$,
  • upper-triangular a matrix if each element under its diagonal is zero.

Because of the orthogonality:

\begin{equation} Q^T \cdot A = R. \end{equation}

Then after applying Q from the right side we get another matrix. Repeating this process more times known as the QR-iteration.

\begin{equation} A_{m+1} = Q_m^T \cdot A_m \cdot Q_m \end{equation}

Again, due to the fact that Q is an orthogonal matrix the spectra of A remains unchanged after one step (eq. 3) of the iteration. After enough step the series of $A_m$-s will converge to an upper-triangular matrix for which the eigenvalues are given as the diagonal elements (think of the determinant). The next question is how to calculate a right Q? This can be achieved by a rotation-like transformation. Each transformation makes an element become zero under the diagonal. The picture below shows such a matrix ($S_{ij}$):

smatrix

In order to calculate $\theta_{ij}$:

\begin{equation} -sin \theta_{ij} \cdot a_{jj} + cos \theta_{ij} \cdot a_{ij} = 0 \rightarrow tg \theta_{ij} = \frac{a_{ij}}{a_{jj}}. \end{equation}

It can be possible that $a_{jj}$ is zero, then $\theta_{ij} = \frac{\pi}{2}$.

Overview of the algorithm

  1. m = 1 and $A_m = A$.
  2. Determine an $S_{ij}$ to create zero under the diagonal of $A_m$ in row i and column j.
  3. Repeat step 2 until all of the elements become zero under the diagonal.
  4. Multiply the $S_{ij}$ matrices to get the $Q_m$ matrix.
  5. Apply equation 3.
  6. Repeat from step 2 until $A_m$ is close enough to an upper-triangular matrix.
  7. The eigenvalues are given as the elements of the diagonal in the last $A_m$ matrix.

For the implementational details and for specific examples see the corresponding jupyter notebook on GitHub.

Code on GitHub