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 ?

  • $\begingroup$ Are you trying to compute $Q_1 Q_1^T$? Because of course $Q Q^T = I$. $\endgroup$
    – vibe
    Jan 12 at 5:59
  • $\begingroup$ That's not the case for m > n .. here A and Q are tall/skinny. $\endgroup$ 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).

  • $\begingroup$ 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. $\endgroup$
    – matthiasbe
    Jan 13 at 8:31
  • $\begingroup$ I encourage you to try it both ways and measure the difference on representative datasets, I could be wrong. $\endgroup$ Jan 13 at 15:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.