How to compute orthonormal basis for rank deficient rectangular matrix

Question

I'm doing a project which needs to compute the orthonormal basis of a rectangular matrix and this matrix may or may not be rank deficient. In matlab, we can just call function orth() which is based on svd to handle this problem, but I need to implement it in C. I tried to use Gram-Schmidt and it works well when the matrix is full-rank. Is there any library in C which can solve this problem? Or some hints on implementing it in C? Thanks a lot.

Answer 1

When a matrix is rank-deficient by k (for instance, if a 4x3 matrix is of rank 2, then k=1), then you need to generate k random vectors and then subtract the components of other orthonormalized vectors previously generated by Gram-Schmidt from the random vectors.

It is necessary to generate vectors yourself as you can see the following example:

A 2x2 matrix

1 1
0 0

By Gram-Schmidt, you will get

1 0
0 0

It is because the original column space only contains a single direction. When you detect such a case (norm(v2) < eps, where eps is a sufficiently small floating-point number, like 1e-10), generate a random vector for the second vector v2, say

1  0.2
0 -0.3

subtract the component of v1

1  0
0 -0.3

and normalize the current v2

1  0
0  1

In such a way, you can construct the directions that don't exist in the original column space.

Note that if computational speed is of your concern, at this point or at your production code, SVD should be avoided when possible. Although SVD and matrix-matrix product are both of O(N^3) complexity, SVD usually has a factor which is larger than 10 times slower compared to matrix-matrix product. In fact, in such a problem, SVD can be useful in your investigation, but a better tool should be QR decomposition, which is essentially Gram-Schidt but in different order of operations to have better numerical stability.

Answer 2

Try the LAPACK library. It has the svd routine (and many others).