# Efficient computation of the extension of a linear basis to completion when the basis is almost complete (ideally using LAPACK routines)

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

## 1 Answer

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.

• p is large, so wont that be time consuming? – wjastle May 17 '12 at 20:10
• Yes, it will be. In which form is your $B$ given? As it is orthonormal, it is likely that you got it in a way that makes the orthonormality explicit. Then this could be exploited. – Arnold Neumaier May 17 '12 at 20:17
• but your $P$ is also expensive to compute. – Arnold Neumaier May 17 '12 at 20:22
• Its a wavelet basis matrix obtained from DWT on columns of the identity matrix. The basis is incomplete because there are more wavelets than datapoints due to boundary effects. – wjastle May 17 '12 at 20:23
• Probably adding a few artificial data points outside the boundary and increasing the interval before doing the wavelet transform will then produce a complete orthonormal basis when restricted to the interval of interest. – Arnold Neumaier May 17 '12 at 20:26