Efficiently compute a projection matrix from Householders reflectors

Question

Let $A \in \mathbb{R}^{m \times n}$ where $m \geq n$.

Let $B$ and $\tau$ be the result of applying LAPACK's dgeqrfp method (R on the upper right triangle, and the Householders reflectors below the diagonal).

How would you use the reflectors and $\tau$ to perform a projection of any matrix on the space orthogonal to the space spanned by the columns of A ?

One way is to build all the Householder transformations $H_i = I - \tau v_i v_i^t$ where $v_i[1:i] = 0$ and $v_i[i+1:m] = B[i+1:m]$ in matlab's notation. We can then compute the factor $Q = H_n ... H_1$.

Then the projector is $P_{A^\bot} = QQ^t$.

Do you know an more efficient way ?

Answer 1

Fortunately, LAPACK provides routines to deal with the $\mathbf Q$ factor from the $\mathbf A = \mathbf Q \mathbf R$ decomposition, [dgeqrf]. To find the projection of an arbitrary $\mathbf B$ onto the space orthogonal to $ \mathrm {range}(\mathbf A)$, you want to form $\mathbf C = \left(\mathbf I - \mathbf Q \mathbf Q^T\right) \mathbf B$. Here are two options:

(1) you can tabulate $\mathbf Q$ explicitly using [dorgqr] and then compute $\mathbf C$ using two calls of [dgemm], first with a transpose, then without. You can form $\mathbf Q$ by overwriting $\mathbf A$, but you'll probably need a temporary to represent the intermediate value $\mathbf Q^T \mathbf B$.

(2) use two calls [dormqr], which can apply the action of $\mathbf Q^T$ or $\mathbf Q$ without explictly forming it. You'll probably still need a temporary to represent $\mathbf Q^T \mathbf B$.

I would expect (2) to be a little more accurate since it works with the householder format directly. But I wouldn't be surprised if (1) is faster for large problems (due to it being to easier to optimize [dgemm] than [dormqr]). Especially if you have multiple $\mathbf B$'s, across which you could amortize the cost of tabulating $\mathbf Q$ up front.

EDIT: I should clarify that on a single operation (multiply by $\mathbf Q$ or multiply by $\mathbf Q^T$), I wouldn't expect [dorgqr] followed by [dgemm] to beat a single call to [dormqr]. Practically speaking, [dorgqr] is the same householder-accumulation algorithm as [dormqr], it's just a multiply operation $\mathbf Q \mathbf I$ applied to a particular identity input (perhaps a little bit faster because accumulating $\mathbf Q \mathbf I$ generates fillin in a predictable/exploitable way). In this context, my belief that (1) could be faster stems more from the fact you need to apply $\mathbf Q$ twice, which gives you a chance for amortization/reuse of the effort you spent computing $\mathbf Q$ across multiple [dgemm] calls (which is just about the fastest operation you can find). We are further aided by the fact that we don't need to hang onto $\mathbf R$, which means we may reuse/clobber $\mathbf A$ to store $\mathbf Q$ (whereas in the general case, you'd probably prefer [dormqr] because it will leave $\mathbf R$ undisturbed).