As a segue to my [prior post on this topic](http://stats.stackexchange.com/a/206345/67822) I want to share some tentative (albeit incomplete and flawed) exploration of the functions behind the linear algebra and related R functions. This is supposed to be a work in progress.

Part of the opaqueness of the functions has to do with the "compact" form of the Householder $\text {QR}$ decomposition, which is meant to increase computational efficiency. In general the Householder decomposition would call for normalization of the vector $v$ meant to represent the difference between the column vectors $x$ in the matrix $A$ that we want to decompose and the vectors $y$ corresponding to the reflection across the subspace or "mirror" determined by $v$. That normalized norm-2 $1$ vector $u$ would then be used to compute the different Householder transformations $ I - 2 uu^T x$. Here is an illustration:

[![enter image description here][1]][1]

Instead, the method used by LAPACK eliminates the matrix multiplications involving the computation of $Q$, and instead liberates the first entry in the Householder reflectors by turning them into $1$ so that the final $u$ vector is not necessarily normalized to $\lVert u\rVert= 1$, but it can still be used as a directional vector. The beauty of the method is that given that $R$ needs to populate the replaced diagonal of $A$ throughout the $QR$ decomposition, we can actually take advantage of the $0$ elements in $R$ below the diagonal to fill them in with these "dwarf" $v$ reflectors, which, having their heading elements all equal to $1$ will have now one less dimension to be represented.

Although these computations are meant to facilitate working with block matrices, the intuition can be gained with a simple random matrix. All the code that follows is [here](https://github.com/RInterested/SIMULATIONS_and_PROOFS/blob/master/HOUSEHOLDER%20QR).

In the "manual" calculation in the code attached we can reproduce the "compact QR" matrix as the addition of the $R$ matrix and the lower triangular "storage" matrix for the reflectors (excluding their $1$ lead elements), calculated as:

$u[2:\text{length(u)}]/\tau$ with $u = \frac{\text{sign}(x_i=x_1)\times \lVert x \rVert \begin{bmatrix}1\\0\\0\\\vdots\\0\end{bmatrix}+\begin{bmatrix}x_1\\x_2\\x_3\\\vdots\\x_m\end{bmatrix}}{\text{sign}(x_i=x_1)\times \lVert x \rVert + x[1]}$ and $\tau = \frac{1 + u[2:\text{length(u)}]^T\,u[2:\text{length(u)}]}{2}$.

In this way $tau$ or other derived scalars, such as $\rho$ (in linked code page) can be easily stored to reproduce $Q$, while in the lower $0$ part of $R$ these "decapitated" $u$ vectors are neatly stored:
[![enter image description here][2]][2]

The equivalence between the manual and R functions are probably self-explanatory in the link to the code attached in Github.


  [1]: https://i.sstatic.net/rrarE.png
  [2]: https://i.sstatic.net/LmydM.png