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

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

• You may use *ORMQR function from Lapack. You should also use *geqrf instead of *geqr2. – Pawel Kowal Feb 17 at 15:17

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.

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.