5
$\begingroup$

Given two matrices $U_1, U_2$, we can use QR factorization to find orthogonal basis for the subspace spanned by (columns of) $\begin{bmatrix}U_1,U_2\end{bmatrix}$.

Now this generally makes no use of the fact that columns of $U_1, U_2$ might already be orthogonal bases themselves (i.e., unit length, orthogonal).

What is the fastest (or at least fast) way to construct orthogonal basis $U$ such that $\text{span}(U)=\text{span}(\begin{bmatrix}U_1,U_2\end{bmatrix})$, if we know that $U_1$ and $U_2$ are already orthogonal bases themselves?

And as a bonus question ;) -- can such algorithm be expressed on top of BLAS routines over $U_1, U_2$, or do we have to get our hands dirty in explicit loops?

EDIT: consistently adding the word "orthogonal" :-)

Improve this question
$\endgroup$
1
  • 1
    $\begingroup$ might already be orthogonal bases themselves (i.e., unit length, orthogonal) -> might already be orthonormal bases themselves (i.e., unit length, orthogonal) ? $\endgroup$
    – robin girard
    Commented Jul 26, 2010 at 12:08

2 Answers 2

Reset to default
1
$\begingroup$

You may think to the situation as a Gram-Schmidt orthonormalization partially performed (the columns of $U_1$ are already orthonormal), and complete it.

The matrix $V=U_2-U_1U_1^TU_2$ has columns orthogonal to $U_1$; then you compute V=QR and $[U_1,Q]$ should be the basis you're looking for.

This exploits the fact that $U_1$ is already orthonormal, but not the same fact for $U_2$.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Nice, cheers. The way I read that formula, it drops complexity from $O(p(r_1 + r_2)^2)$ (QR of $U^{p \times (r_1+r_2)}$ to $O(p\text{ min}(r_1, r_2)^2)$ (QR of whichever $V$ is smaller), which is not bad but I think we should be able to do better. $\endgroup$
    – RedSnow
    Commented Jul 26, 2010 at 15:34
  • $\begingroup$ You're missing something in your computational cost, the first matrix product costs $O(pr_1r_2)$. I don't think you can save more than that, though. $\endgroup$
    – Federico Poloni
    Commented Jul 27, 2010 at 17:23
  • $\begingroup$ True, so the complexity is $O(\text{matrix product} + \text{QR})=O(pr_1r_2 + p\text{ min}(r_1, r_2)^2=O(pr^2)$ for $r\approx r_1 \approx r_2$. Now that's the theoretical big $O$. But where can the orthonormality of $U_2$ be used, to give at least some boost for practical applications? The QR part is obviously much slower than the multiplication, so never mind the matrix product for now. $\endgroup$
    – RedSnow
    Commented Jul 28, 2010 at 5:55
  • $\begingroup$ If $p\gg \max(r_1,r_2)$, then the coefficient in front of both QR and matrix product is 2, so the matrix products are indeed slower than the QR --- a good "cheatsheet" reference for the computational costs is Appendix C of Higham, "Functions of matrices". Other things matter such as the so-called "level-3 factor", but I think in a good BLAS implementation the costs are of the same order of magnitude. Frankly I don't think the orthogonality of the other matrix can be exploited; if you manage to exploit it, let me know because the same idea could apply to an algorithm I'm working on. :) $\endgroup$
    – Federico Poloni
    Commented Jul 28, 2010 at 8:13
1
$\begingroup$

I have an idea for an improvement -- at least for the application I had in mind, where this merging of bases takes place recursively:

After GEQRF routine (QR decomposition of $V$), keep the basis $Q$ in the form of elementary reflectors -- do not convert to explicit matrix with ORGQR. Then, in forming the subsequent $V=U_2-U_1U_1^TU_2$, perform all matrix multiplications with ORMQR (not GEMM).

I'd love to try and time this, to compare speed of ORMQR vs. ORGQR+GEMM, but NumPy doesn't wrap the necessary LAPACK routines, and I'm not familiar with any other software package :(

If anyone can verify efficiency and attest whether this makes any sense, I'll be happy to accept your answer and be done with this question, as it seems nobody has any other ideas.

Improve this answer
$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .