aliquote.org

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 eignevalues on its diagonal.

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

♪ Clan of Xymox • A Day

See Also

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