Looking at http://eigen.tuxfamily.org/dox-devel/gr ... stems.html, it seems like SparseQR is the only solver that supports rectangular problems. The short form of my question is does Eigen support a rectangular solver that uses the normal equations backed by a CG solver that doesn't compute A^t A.

Specifically, I have an overconstrained sparse linear solve. The three ways I've looked at are:

(1) SparseQR -> The problem I'm working with has a very challenging and big QR factorization; perhaps it is possible to make the QR factorization work on my problem. I found some good advice here: viewtopic.php?f=74&t=125165 and I could provide my matrix if people are interested. But I think that for many problems (and mine specifically), solving with QR isn't a good idea, and seems to take upwards of hours.

(2) Normal equations by direct multiplication -> When working on smaller systems with Eigen, I typically just did solveCG(A.transpose() * A, A.transpose() * b), which works just fine. However, for my current problem the system is very big and about 95% of the time is spent in the computation of A.transpose() * A (and the app will use upwards of 15GB of memory). The actual CG solve at the end once A^tA is computed is done very fast (less than a minute.)

(3) Normal equations by a manual CG solver -> The solution I've currently gone with is to use a quick CG solver I wrote that replaces any matrix-vector product with two calls to Eigen's matrix-vector multiply routines. This is working fine, but is a bit clunky (I do all the preconditioning etc. myself) and probably much less elegant than Eigen's CG solver. This seems like a pretty common thing to want to do (I think Eigen may support an A^tA view for DenseMatrix?), so I'm wondering if Eigen already has some native support for this sort of thing. The interface seems exactly like the SparseQR solver, just using the normal equations and Eigen::ConjugateGradient on the A^tA view instead.

In order to reduce numerical issues, least-square CG is a bit more subtle than simply replacing A*p by A.adjoint()*(A*p) and b by A.transpose()*b [1]. So we would need a different class for that, but that would definitely be a useful addition to Eigen.

[1] https://web.njit.edu/~jiang/math614/hes ... tiefel.pdf

A patch is in the pipe: http://eigen.tuxfamily.org/bz/show_bug.cgi?id=975, please give it a try and let us known if it's working well for you.

Looks great!

This will be very useful. For many problems I think LSCG is better than QR (especially ones where there is memory pressure or a good variable reordering can't be found.)