How to solve $Q$ matrix from Householder QR-factorization? - Lapack

Question

I'm using the subroutine sgeqr2 from Lapack. This subroutine solves the QR-factorization

$$A = QR$$

It's easy to find the $R$ matrix, because the in-out argument $A$ of subroutine sgeqr2 will return a matrix, where the upper values from the diagonal is the $R$-values. Eeasy.

But how can I find the $Q$ matrix? According to the "manual".

The matrix Q is represented as a product of elementary reflectors

     Q = H(1) H(2) . . . H(k), where k = min(m,n).

  Each H(i) has the form

     H(i) = I - tau * v * v**T

  where tau is a real scalar, and v is a real vector with
  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
  and tau in TAU(i).

The conclusion here is to take the vector $tau$ and pick the first value of $tau$. Then multiply it with a matrix $vv^T$.

Have I interprent this text correct?

For a matrix $A$

0.674878,   0.151285,   0.875139,   0.150518,
0.828102,   0.150747,   0.934674,   0.474325,
0.476510,   0.914686,   0.740681,   0.060455,
0.792594,   0.471488,   0.529343,   0.743405,
0.084739,   0.475160,   0.419307,   0.628999,
0.674878,   0.151285,   0.875139,   0.150518

I get the $R$ matrix

-1.568159 -0.751743 -1.762103 -0.808132 
0.000000 0.887756 0.233565 0.241302 
0.000000 0.000000 0.500422 -0.142971 
0.000000 0.000000 0.000000 -0.700355 
0.000000 0.000000 0.000000 0.000000 
0.000000 0.000000 0.000000 0.000000

And $tau$ vector

1.430363 1.205732 1.007224 1.655577 

How should I do, to find $Q$ matrix if I know the dimension of $A$ and the values from $tau$?

I wrote som MALAB/Octave code to find the $Q$-matrix from $H$, but it won't work.

tau = [1.430363 1.205732 1.007224 1.655577];

H = eye(4); % Initial 

for i = 1:4
  v(1:i-1) = 0;
  v(i) = 1;
  H = H*(eye(4) - tau(i)*v'*v);
  H
end

Answer 1

Now I got the answer!

Assume that we have a matrix $A$

0.674878,   0.151285,   0.875139,   0.150518,
0.828102,   0.150747,   0.934674,   0.474325,
0.476510,   0.914686,   0.740681,   0.060455,
0.792594,   0.471488,   0.529343,   0.743405,
0.084739,   0.475160,   0.419307,   0.628999,
0.674878,   0.151285,   0.875139,   0.150518

And we implement it into the subroutine.

Then we get the matrix $A$ back, but it has change into this

-1.568159 -0.751743 -1.762103 -0.808132 
0.369188 0.887756 0.233565 0.241302 
0.212440 -0.675308 0.500422 -0.142971 
0.353358 -0.142374 0.875625 -0.700355 
0.037779 -0.412039 -0.411443 0.439228 
0.300877 0.112496 -0.222825 -0.122951 

The $R$ matrix is the upper triangular above the diagonal.

-1.568159 -0.751743 -1.762103 -0.808132 
0.000000 0.887756 0.233565 0.241302 
0.000000 0.000000 0.500422 -0.142971 
0.000000 0.000000 0.000000 -0.700355 
0.000000 0.000000 0.000000 0.000000 
0.000000 0.000000 0.000000 0.000000 

And then the lower traingular under the diagonal, we want to use those in this loop.

Q = eye(6); % Initial 
v = [ 0 0 0 0 0 0];
m = 6;
for i = 1:4
  v(1:i-1) = 0;
  v(i) = 1;
  v(i+1:m) = A(i+1:m,i);
  A(i+1:m,i)
  Q = Q*(eye(6) - tau(i)*v'*v);
end

Where tau is this vector

1.430363 1.205732 1.007224 1.655577

Then our $Q$ matrix will be.

  -0.430363  -0.194015   0.323945   0.148698  -0.562826  -0.577340
  -0.528073  -0.277360   0.137758  -0.191610   0.739060  -0.205691
  -0.303866   0.773024   0.049331   0.520577   0.184471  -0.051341
  -0.505430   0.103108  -0.770067  -0.285532  -0.234751   0.065334
  -0.054038   0.489479   0.419171  -0.752685  -0.118974   0.033112
  -0.430363  -0.194015   0.323945   0.148698  -0.183642   0.785092

I post code instead of matrix math, because it will be easier to read.

Answer 2

Your implementation is very inefficient.

Consider

Q = Q*(eye(6) - tau(i)*v'*v);

This would create the matrix (eye(6) - tau(i)*v'*v) and then multiply it by Q.

If your matrices are n x n, then this would require $ 2 n^3 $ computations for that multiply.

If you instead update Q with

Q = Q - tau(i) * ( Q * v' ) * v then the computation is $ O( n^2 ) $ per iteration.

There are other inefficiencies in what you do: You have a lot of zeroes in your vector v that you could take advantage of.

But, of course, in your example you are only working with a $ 6 x 6 $ matrix, so it doesn't really matter so much.