C++: Solve underdetermined linear system with lapack

Question

Let's put it simple: I have an under-determined linear system of equations

Ax = b

and I want to get one valid solution, no matter which one of the infinite solutions for the system. And I want to get it as efficiently as possible.

I have checked general LAPACK routines and it seems that they cannot handle the under-determined case. For example, dgesv(), whose documentation is found here, will return and integer larger than 1 in INFO if the factor U, from PLU factorization, is singular, and it will not solve the system if that is the case.

I have also checked some routines for Linear Least Squares problems (LLS) (documentation here), which does exactly solve my problem, just not as efficiently as I wished. If the LLS problems I provide is under-determined, the LLS routine will return the vector that minimizes

||Ax-b||

Which is a valid solution. However, it is calculated as the solution to an optimization problem, and I was wondering if there is a more efficient way of obtaining a valid solution for my under-determined problem.

A similar question was made here, but I believe that my question is more concrete than that: I am using LAPACK, and I want to solve an under-determined system of linear equations as efficiently as possible.

Answer 1

For an under-determined system of equations:

The correct approach is to use singular value decomposition (SVD). Lapack offers singular value decomposition in the form of dgesvd.

To perform the SVD you will have to homogenize your problem to turn it into a matrix problem of the form: My = 0. This is easy to do by introducing another degree of freedom (another variable). This will transform the vector x -> y and matrix A -> M. When performing the SVD on the matrix M, the smallest singular vector will be the solution to your under-determined least squares problem.

I would recommend using matlab or octave to experiment before wasting time coding anything up.