fast merging of orthogonal bases 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" :-)
 A: 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$.
A: 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.
