I have a $p \times n$ matrix $B$ (where $n < p$) with orthonormal columns and would like to find a numerically efficient way to extend this matrix to get a complete $p$-dimensional orthonormal basis. In other words I want to compute a $p$ by $p-n$ complementary matrix $W$ so that $[B|W]$ is orthonormal. At the moment I am computing the QR factorisation of the rank $p-n$ projection matrix$P=I-BB^T= QR$ using the LAPACK routines dgeqp3 and dorgqr and discarding the trailing $p$ columns of $Q$. I feel sure that there must be a more efficient method because in my application $B$ is almost a complete basis i.e $p-n$ is small.  I have wondered about building a complementary basis $R$ recursively by working through the columns of $P$. At each step adding a column $C$ to the current $R$ only if $\det([R|C]^T[R|C)\neq0$. After I have $p-n$ independent vectors in $R$ I could form $W$ as $W=(R^TR)^{-1/2} R$. However, I'm not sure if this is a numerically stable solution, or if there is another more efficient method.  Ideally I would like a method I can implement using standard LAPACK routines. Originally asked on stackoverflow here
Efficient computation of the extension of a linear basis to completion when the basis is almost complete (ideally using LAPACK routines)
linear algebra;basis set;lapack
Just compute a QR factorization of $B$. The last $p-n$ columns then have he required property of forming a basis complement.For large $n$ and small $p-n$, you can take a random $p\times k$ matrix $X$ with $k=p-n$ or slightly larger (for an increased likelihood of having numerical stability), and compute a QR factorization of $X-B(B^TX)$, which, with probability one, spans the complementary subspace.