4
$\begingroup$

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

$\endgroup$
2
$\begingroup$

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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