# Efficiently compute a projection matrix from Householders reflectors

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 ?

• Are you trying to compute $Q_1 Q_1^T$? Because of course $Q Q^T = I$. – vibe Jan 12 at 5:59
• That's not the case for m > n .. here A and Q are tall/skinny. – rchilton1980 Jan 12 at 17:17

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).

• Interesting, thank you for your answer. I wasn't aware the matrix-matrix multiplication could be faster than using the compact QR format with the householders reflectors. – matthiasbe Jan 13 at 8:31
• I encourage you to try it both ways and measure the difference on representative datasets, I could be wrong. – rchilton1980 Jan 13 at 15:09