Computing eigenvalues using QR decomposition

December 28, 2022

Matthew Drury made a great job explaining how eigenvalues are usually computed, via the QR decomposition. This is an iterative algorithm which is available in Numpy, numpy.linalg.qr, while in R it is directly qr. We already discussed the QR decomposition in a previous post. Basically, the idea is as follows: Let $X$ be a symmetric matrix, compute $X_k = Q_kR_k$, for $k=1,\dots,n$, and update $X_{k+1}=R_kQ_k$. The sequence of $X_n$ converges to a diagonal matrix $D$ of eigenvalues, such that the eigenvectors are the columns of $\prod_i Q_i$.

In Racket, the very first iteration would read:

(require math)
(define X0 (matrix [[12 -51   4]
                    [ 6 167 -68]
                    [-4  24 -41]]))
(define X (matrix+ X0 (matrix-transpose X0)))
; (define D (diagonal-matrix '(1 1 1)))
(define D (identity-matrix 3))
(define-values (Q R) (matrix-qr X))
(define D (matrix* D Q))
(define X (matrix* R Q))

Wrap the above code in a function, iterate 20 to 30 times and then the $X$ matrix should contain the eigenvalues on its diagonal.

Of note, the Householder method for computing the QR decomposition in Racket is available on Rosetta.

For other numerical approaches to computing eigenvectors and eigenvalues, see this excellent blog post by Marc Khoury: Numerical Algorithms for Computing Eigenvectors.

♪ Clan of Xymox • A Day

See Also

» Prime palindromes » Gaussian elimination » Armstrong numbers » Perfect and amicable numbers » Decimal numbers