Some fun uses of the QR decomposition

Computing the QR decomposition
Computing eigenvalues

Computing the QR decomposition

Here's the matrix we played with in class, together with it's QR decomposition.

A = [-1,3;2,6;2,3];
[Q,R] = qr(A)

Q =

   -0.3333   -0.8666    0.3714
    0.6667   -0.4952   -0.5571
    0.6667    0.0619    0.7428


R =

    3.0000    5.0000
         0   -5.3852
         0         0

Note that the first column of Q and the first row of R agree with our computations. qr returns the full QR decomposition. You can get get the reduced QR decomposition as follows.

[Q1,R1] = qr(A,0)

Q1 =

   -0.3333   -0.8666
    0.6667   -0.4952
    0.6667    0.0619


R1 =

    3.0000    5.0000
         0   -5.3852

Cleve Moler's qrsteps function lets us see the evolution of the matrix A.

qrsteps(A);

A =

    -1     3
     2     6
     2     3


A =

     3     5
     0     5
     0     2


A =

    3.0000    5.0000
         0   -5.3852
         0         0

We can also automate the procedure ourselves.

a1 = A(:,1);
v = zeros(length(a1),1);
v(1) = norm(a1);
v = v-a1;
Q1 = eye(length(a1)) - 2*v*v'/(v'*v)
Q1*A

Q1 =

   -0.3333    0.6667    0.6667
    0.6667    0.6667   -0.3333
    0.6667   -0.3333    0.6667


ans =

    3.0000    5.0000
    0.0000    5.0000
    0.0000    2.0000

Computing eigenvalues

It turns out that the QR decomposition is the basis of a very important technique for computing eigenvalues. Here, for example, is a matrix together with some of its eigenvalues.

rand('seed',1);
A = round(10*rand(5)-5);
eig(A)'

ans =

  Columns 1 through 4

   1.1874             0.6528            -0.5414            -5.1494 - 1.4426i

  Column 5

  -5.1494 + 1.4426i

And here's an iterative process using the QR decomposition.

B = A;
for i = 1:100
    [Q,R] = qr(B);
    B = R*Q;
end;
B

B =

   -4.3415    4.1922    1.6496   -0.4477    0.2405
   -0.6521   -5.9573    3.8484   -3.7543   -3.0515
    0.0000    0.0000    1.1874   -2.1253    1.5801
    0.0000    0.0000    0.0000    0.6528    4.9318
   -0.0000   -0.0000   -0.0000   -0.0000   -0.5414

Note how the real eigenvalues of A lie on the diagonal of B. The complex eigenvalues are the eigenvalues of the obvious block.

eig(B(1:2,1:2))'

ans =

  -5.1494 - 1.4426i  -5.1494 + 1.4426i

Our next step in this class is to try to figure out why this works. The first thing we need to come to grips with is the fact that it doesn't always work!

rand('seed',2);
A = round(10*rand(5)-5);
B = A;
for i = 1:1000
    [Q,R] = qr(B);
    B = R*Q;
end;
eig(A)'
B

ans =

  Columns 1 through 4

   4.4338            -3.1369 - 3.1162i  -3.1369 + 3.1162i  -4.5800 - 0.6594i

  Column 5

  -4.5800 + 0.6594i


B =

   -4.7444   -5.5839    1.1236    3.8530    5.7388
    0.0827   -4.4156   -3.2431   -6.4272    1.9030
   -0.0000    0.0000    4.2094    1.7134    1.5066
   -0.0000    0.0000    0.6950   -4.4060    4.9983
   -0.0000    0.0000    0.4548   -2.3675   -1.6434

There are however situations where it works exceedingly well.

rand('seed',2);
A = round(10*rand(5)-5);
A = A*A';  % Ensure that input is symmetric.
B = A;
for i = 1:100
    [Q,R] = qr(B);
    B = R*Q;
end;
eig(A)'
B

ans =

    1.4680    6.8778   29.2699  101.1661  115.2182


B =

  115.2182    0.0000    0.0000    0.0000   -0.0000
    0.0000  101.1661   -0.0000   -0.0000   -0.0000
   -0.0000    0.0000   29.2699    0.0000    0.0000
    0.0000   -0.0000   -0.0000    6.8778    0.0000
   -0.0000    0.0000    0.0000   -0.0000    1.4680

There are also techniques to accelerate the convergence.

rand('seed',1);
A = round(10*rand(5)-5);
eig(A)'

ans =

  Columns 1 through 4

   1.1874             0.6528            -0.5414            -5.1494 - 1.4426i

  Column 5

  -5.1494 + 1.4426i

B = A;
n = length(B);
for i = 1:10
    [Q,R] = qr(B);
    B = R*Q;
end;
B

B =

   -3.3171    3.0790    3.0166   -2.0590   -0.4402
   -1.7674   -6.9796    2.9117   -3.9308   -1.9386
    0.0013    0.0058    1.1869   -1.6333    2.0805
    0.0000    0.0000    0.0003    1.8194    4.2902
   -0.0000   -0.0000   -0.0000   -0.6430   -1.7096

B = A;
n = length(B);
for i = 1:10
    mu = B(n,n);
    [Q,R] = qr(B - eye(n));
    B = R*Q + eye(n);
end;
B

B =

   -4.1354   -4.0850   -1.1381   -1.5944    0.3535
    0.7615   -6.1608    1.3881    3.8827   -4.4980
   -0.0015    0.0105   -0.5440    2.9729    3.2846
    0.0000   -0.0000   -0.0000    0.6513   -3.4234
    0.0000   -0.0000   -0.0000    0.0002    1.1889

Some fun uses of the QR decomposition

Contents

Computing the QR decomposition

Computing eigenvalues