4

I'm trying to implement the tiled QR decomposition in MATLAB (in an attempt to understand it), and I'm trying to check that my SGEQRF (upper corner tiles) function is working correctly. I have a similar C implementation that does work correctly (for SGEQRF) but I get a slightly different answer for some matrices, and I'm trying to understand how to reclaim Q from the H vectors.

For example on the matrix A = $$ \left( \begin{array}{ccc} 1 & 1 & 1 & 1\\ 2 & 2 & 2 & 2 \\ 3 & 3 & 3 & 3\\ 4 & 4 & 4 & 4 \end{array}\right) $$

my C implementation retrieves a resulting matrix: $$ \left( \begin{array}{ccc} -5.477 & -5.477 & -5.477 & -5.477\\ 0.309 & 0 & 0 & 0 \\ 0.463 & 0 & 0 & 0\\ 0.618 & 0 & 0 & 0 \end{array}\right) $$

with tau matrix: $$ \left( \begin{array}{ccc} 1.183 \\ 2 \\ 2\\ 0 \end{array}\right) $$

When I used the LAPACK SGEQRF function on the same matrix, it gets different tau and result matrices, however using SOGRQF on my result does return a Q matrix such that Q*R = A, so am satisfied that implementation does work.

What I don't understand however, is how to recover Q from the resulting matrix (ideally I would like to implement this myself in MATLAB so I can understand whats going on).

Presumably you can compute $Q_1$ from $\tau_1 $, $r_{1,1}$ and $v_1$ (where $v_1$ is the lower diagonal first column of the resulting matrix) somehow but I can't see how.

For my example it seems like:

$$ \left( \begin{array} $1 \\ v_0 \\ \end{array} \right)\ $$ $\times \tau_0 \times r_{0,0} +$ $$ \left( \begin{array} $r_{0,0}\\ 0 \\ \end{array} \right)$$

is approximately equal to the first column of A. (Sorry I can't get that latex sum to work for some reason :/ ).

However I'm not sure this is the case always, and the first value doesn't seem to quite be correct (I get 1.0051 instead of 1, which seems to be more than rounding error likely). Further I'm not sure how to propagate this through the matrix.

How does this change when you have the tiled version, where each tile below the corner has its own $\tau$ value?

Thanks

LonelyCat
  • 71
  • 5

1 Answers1

3

Ok I worked this out for the none tiled case. Each matrix $Q_1$ is defined by $Q_1 = I - 2 v v^T$ for the traditional QR decomposition. In this case the vector $v$ is called $w$ and stored in the lower section of the matrix - so for the example above: $$ w =\left( \begin{array}{ccc} 1 \\ 0.309 \\ 0.463\\ 0.618\\ \end{array} \right) $$

We can then recover each of the householder matrices $Q_i'$ by $Q_i' = I - \tau_iw_iw_i^T$ where $Q_i'$ represents the lower square of the matrix, with the identity above.

To regenerate $Q$, we simply do $Q = Q_1^TQ_2^T...Q_n^T$.

Edit: My idea for tiled was wrong so I'm going to think more on that, but removed the incorrect section for now.

LonelyCat
  • 71
  • 5