The QR factorization and alogorithm

First, let $A$ be a real square matrix with column vectors $A_i$. We use the Gram-Schmidt process to create a orthogonal basis: $\{U_i\}$:

$$~ U_1 = A_1; \quad U_2 = A_2 - proj_{U_1}(A_1); \dots; U_n = U_n - \sum_{j < n} proj_{U_j}(A_n). ~$$

An then let $E_i = U_i/\|U_i\|$ be an orthonormal basis. If we express the $A_i$s in terms of the $E_i$s we get $A = QR$ where

$$~ Q = [E_1 E_2 \cdots E_n]; \quad R = E_i \cdot A_j \text{ when } i \leq j \text{ otherwise } 0. ~$$

Note that $Q$ has a special property, namely $Q^TQ = I$, so $Q^T = Q^{-1}$. Such $Q$ are called orthogonal and for complex valued $Q$, unitary. This gives:

The QR factorization: Any real square matrix can be decomposed as $QR$ where $Q$ is orthogonal and $R$ is upper triangular.

In Julia, the function qr will return the two matrices in the factorization:

using LinearAlgebra
A = [1 2; 3 5]
Q,R = qr(A)
Q

2×2 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:
 -0.316227766016838   -0.9486832980505138
 -0.9486832980505138   0.316227766016838

2×2 Array{Float64,2}:
 -3.1622776601683795  -5.375872022286245 
  0.0                 -0.3162277660168371

Q*R - A

2×2 Array{Float64,2}:
 4.440892098500626e-16  0.0
 0.0                    0.0

The QR algorithm

The $QR$ algorithm use the $QR$ factorization to find an upper triangular matrix similar to $A$. Similar matrices have the same eigenvalues and the eigenvalues of an upper triangular matrix are the diagonal entries, so this can be used to iteratively find the eigenvalues of $A$. Let's see how.

Question 1:

Let $A$, $P$, and $B$ be $n\times n$ matrices; and suppose a) $P$ has an inverse; b) $P^{-1}AP = B$. Show using algebra that if $\lambda$ is an eigen vector of $A$ that $\lambda$ is an eigen vector of $B$.

Question 2:

Suppose $A = QR$, as above. Show that $B =RQ$ is similar to $A$.

The QR algorithm

The QR algorithm is an iterative algorithm that use the QR factorization to find an upper triangular matrix similar to $A$ and hence reveal the eigenvalues of $A$.

Here we define $A$ to be a $10 \times 10$ matrix with real eigenvalues $1,2, \dots, 10$.

D = diagm(0 => 1:10)
S = (rand(10,10) .- 0.5)*2
A = S * D * inv(S)

10×10 Array{Float64,2}:
 -22.71983836667061    33.56871064290637     -76.70137962765264     -69.15429630758342    55.79781663212105    -28.146582986294245   -64.95757792474657      72.1352395127158      -44.987653611175745    20.744466921183914 
  11.823353097621512  -11.393519972020245     40.23724575331813      38.02396409450759   -29.814708191143485    16.543713105746992    31.424320324608807    -38.75614473769216      23.129227317296753   -11.690194846457155 
  22.34240067109741   -30.65659440125156      69.26584448494376      57.89145740529617   -43.329444361391       21.895784748803862    55.75119080904199     -60.6188645231603       34.52388111890334    -14.116183464288634 
  47.73794134643261   -61.324431340893234    126.70849835680573     117.154121105631     -87.20848027633124     40.5047405386003     111.92441357113532    -124.03658579100116      68.7005674792452     -30.862083902777083 
   4.69882995681666    -6.222994118741145     12.266686709553255      9.334431455658436   -2.4434716520614685    1.222905664006495    10.277114058552485    -10.964700638936531      7.174778604717872    -3.362972325352626 
 -51.59146617597165    64.74650483923645    -134.2918761022453     -119.24053404604481    92.65362371993683    -38.60604129219455   -119.74910545908486     130.58643086212228     -74.17018540312591     32.465491620706246 
   6.924098738510573   -8.801617860265281     25.33525949523694      22.76995717742598   -20.86362758940716     10.356259097836812    22.731400971545057    -21.931633094339446     17.378869636166343    -9.389576123014688 
  24.569473871225647  -32.95220516298983      69.40613273119935      60.20115332432265   -48.12022321783706     21.919636968981102    59.88923743819131     -60.47606658626965      38.96828468482843    -17.34719045683823  
   1.750362141557467   -0.8193636071887695     0.1750071986963846    -3.176179026597268   -0.3676067509441836   -3.42382213490154      1.3344610941520614    -0.36544859743420943    6.7786491979497185   -0.9448225119599609
  49.437101360639474  -64.57572954175458     131.63723404025848     116.55874689623136   -92.58456977259061     43.25551676992922    116.15992651835697    -129.39088100404456      70.81024642253944    -25.291077890853003

Here $A$ is similar to a diagonal matrix $D$ by construction, the diagonal values of $D$ are as specified.

The $QR$ algorithm starts with $A_1=A$. We form $A_2$ by

First finding $Q_1 R_1 = A_1$ through the QR-factorization
forming $A_2 = R_1 Q_1$. By the second problem above, $A_2$ is similar to $A_1$, so will share the same eigen values:

A1 = A
Q1, R1 = qr(A1)
A2 = R1 * Q1
sort(eigvals(A1))

10-element Array{Float64,1}:
  0.9999999999991012
  2.000000000000222 
  3.000000000000071 
  4.0000000000006155
  4.999999999999847 
  6.000000000000937 
  6.999999999999327 
  7.999999999999745 
  9.000000000000053 
 10.000000000000108

The magic of the algorithm is that if we repeat this to form a sequence $A_1, A_2, A_3, \dots$ this will converge to an upper triangluar matrix. This animation shows every 5 steps of 100 steps of the above thresholding values at 1e-3.

Question 3

Using Julia (or MATLAB if you want) do the above calculation 25 times to form $A_{25}$. Print out your matrix. How many values in the lower triangular part have size less than 1e-3?