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


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 1


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.

  • $\begingroup$ p is large, so wont that be time consuming? $\endgroup$
    – wjastle
    May 17, 2012 at 20:10
  • $\begingroup$ 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. $\endgroup$ May 17, 2012 at 20:17
  • 1
    $\begingroup$ but your $P$ is also expensive to compute. $\endgroup$ May 17, 2012 at 20:22
  • $\begingroup$ 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. $\endgroup$
    – wjastle
    May 17, 2012 at 20:23
  • $\begingroup$ 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. $\endgroup$ May 17, 2012 at 20:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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