QR factorization and eigenvectors

Question

I really can't find an answer to this. I am writing a C code to solve the eigenvalues problem for a real symmetric matrix A. The procedure is the following.

I have a real symmetric matrix A. I put it in tridiagonalized form, P being the matrix that tridiagonalizes A. Then I write the QR factorization of the tridiagonal matrix and I find eigenvalues and eigenvectors of the final (upper) triangular matrix thus obtained. Up to now everything is fine. But from this point, how can I relate these eigenvalues and eigenvectors to those of the original matrix A?

I read of the QR algorithm but I am looking for a more direct approach. I solved the eigensystem of the triangular matrix yet and I am looking for a relation (involving matrices P and Q) that gives the same for the matrix A. How can I proceed?

Answer 1

Edit: Per the discussion in the comments, I misunderstood the question. The intention of the question is to know whether it is possible to easily recover the eigenvalue decomposition of $A$ if we know the QR decomposition of $A$, and we know the eigenvalue decomposition of $R$.

Basically, the answer is no, it is not possible to recover the eigenvalues of $A$ from this information using any finite algebraic process. To see why, note that general-purpose eigenvalue algorithms always require an infinite iterative process. This is because the solutions of a polynomial equation are the eigenvalues of companion matrix for the polynomial, and by the Abel-Ruffini theorem there is no algebraic solution of 5th order polynomials in general.

But here everything is (in principle) a finite algebraic process:

1. Taking a QR factorization of $A$,
2. Getting the eigenvalues, $\lambda_R$ of $R$ by looking at its diagonal.
3. Finding the eigenvectors of $R$, by finding vectors in the null space of $A - \lambda_R I$

If one could use $Q$ and the eigenvectors and eigenvalues of $R$ in an algebraic process to get the eigenvalues of $A$, this would imply that the eigenvalues of $A$ are algebraic in the entries of $A$, which would imply the existence of algebraic solutions to any polynomial, which would contradict the Abel-Ruffini theorem.

---

Old answer:

The QR algorithm for eigenvalues uses the QR decomposition repeatedly in an iterative fashion. At each step you compute the QR decomposition of the current matrix, then multiply the factors in reverse order to get the next matrix. \begin{align} &A_1 = A \\ Q_1 R_1 = A_1, \quad &A_2 = R_1 Q_1 \\ Q_2 R_2 = A_2, \quad &A_3 = R_2 Q_2 \\ Q_3 R_3 = A_3, \quad &A_4 = R_3 Q_3 \\ \dots \end{align} Eventually this process will converge to a triangular matrix with the eigenvalues on the diagonal. You don't just do the QR decomposition one time. This is a common misconception caused by bad naming. For more details, see here:

"Lecture Notes on Solving Large Scale Eigenvalue Problems." Peter Arbenz. Chapter 4. https://people.inf.ethz.ch/arbenz/ewp/Lnotes/chapter4.pdf

Here is an example in python/numpy where I do the QR algorithm for 20 iterations on a 4x4 symmetric matrix:

import numpy as np

A = np.random.randn(4,4)
A = A + A.T
print('Initial A:')
print(A)
true_eigs, _ = np.linalg.eigh(A)
for k in range(20):
    Q, R = np.linalg.qr(A)
    A = R @ Q

print('Final A:')
print(A)
qr_eigs = A.diagonal()
print('true_eigs=', true_eigs)
print('qr_eigs=', qr_eigs)

This outputs the following:

Initial A:
[[ 3.69583037  0.57095943  1.25982728 -0.91138766]
 [ 0.57095943 -5.72767862  1.40997728  3.09590731]
 [ 1.25982728  1.40997728  0.84352098 -0.00810195]
 [-0.91138766  3.09590731 -0.00810195  2.71813063]]
Final A:
[[-7.00604910e+00  3.73685766e-02 -1.74171532e-05  3.95274785e-16]
 [ 3.73685766e-02  4.60446175e+00 -3.51557676e-03 -6.65769624e-17]
 [-1.74171532e-05 -3.51557676e-03  3.57531666e+00  4.29698141e-16]
 [-9.44786542e-25 -2.34167099e-22  1.01641538e-19  3.56074064e-01]]
true_eigs= [-7.00616937  0.35607406  3.57530465  4.60459402]
qr_eigs= [-7.0060491   4.60446175  3.57531666  0.35607406]