It is a conjugate gradient like iterative method where the heavy operation is the product with the sparse matrix A and its transpose.
LSQR solves one of the following cost functions:
This method is appealing from two reasons:
1) It applies to rectangular matrices - both over determined and under determined.
2) It can use a warm start initial vector. This is useful when the system slightly changes over time and you want to still keep track of the solution without fully solving the system from scratch.
The C code is downloadable from here: http://www.stanford.edu/group/SOL/software/lsqr.html
(I strongly advise against using the C++ version which is overly complicated).
0) It is advised to start first with the matlab/octave version which has surprisingly very good performance. Only once you are happy with algorithm performance it is advised to switch to C.
1) You have to implement the matrix vector product operator yourself:
Compute the product
If mode = 1, compute y = y + A*x
If mode = 2, compute x = x + A(transpose)*y.
2) For a warm start, I am using the following suggested solution:
// Note that x is not an input parameter.
// If some initial estimate x0 is known
// one could proceed as follows:
// 1. Compute a residual vector r0 = b - A*x0.
// 2. Use LSQR to solve the system A*dx = r0.
// 3. Add the correction dx to obtain a final solution x = x0 + dx.
3) It is advised to link against cblas to get additional performance gain. Use the
library flags: -lgsl -lgslcblas
For a system of size 500x5000 with 150K non zeros, each iteration of LSQR takes 6 milliseconds!
(Using a single core Intel 2.56GHz).