FullPivLU and QR refuse to solve for non square matrix

[nestordemeure]

The doc says that LU will only solve for square matrix but that both FullPivLU and QR should work for non-square matrices.

In practice there is an assert in both implementations of solve to insures the methods fails on non-square matrices.

Is it because the implementation of FullPivLU and QR is not able to deal with non-square matrix yet or is it a copy-paste error (in which case I can submit a PR) ?

[sebcrozet]

It appears the documentation is not accurate. The construction of LU, FullPivLU and QR for non-square matrices are all implemented. However their methods related to inverse/solve are only implemented for square matrices.

[nestordemeure]

This means that, if I am not mistaken, the only way to currenlty solve for a non-square matrix is to use a SVD decomposition.

It feels problematic to have only one costly algorithm for this use case.

[Andlon]

@nestordemeure: Can you elaborate on what exactly you mean by solving a non-square system? Do you mean in a least squares sense?

I believe that if you want to use e.g. QR decomposition or LU decomposition to solve least squares problem, you need certain guarantees about the factors that you get back. For example, if I remember correctly, for the QR decomposition you would need a rank-revealing QR decomposition. The reason is that the R factor is also non-square, and so the usual back-substitution procedure used for solving linear systems after factorization is no longer directly applicable. Only if the R factor has certain structure (which a rank-revealing QR decomposition would give) would you be able to actually easily solve the (non-square) upper triangular problem.

To summarize, if you want to solve non-square systems a specialized algorithm for this purpose may be necessary.

[nestordemeure]

Yes, I mean it in the least square sense.

Reading the documentation I believed that the method implemented deals with non-squares matrix (which might require using one of two underlying implementations as, indeed, the algorithm is different). As sebcrozet said, it is an inacuracy in the documentation.

On a larger scale, it appears that there is currently only one algorithm to solve non-square linear systems in nalgebra and, to the best of my knowledge, it is a costly one. That feels like a hole in the API.

We could probably take inspiration from Eigen and implement some of the algorithms from the Solving linear least squares systems
and Solving Sparse Linear Systems
pages (I am currently busy but I should be able to implement some of those algorithms in nalgebra).

[kmsquire]

Just wanted to say that we just ran into this issue. We were porting an algorithm from C++ which used QR composition to solve a non-linear system of equations in the least-squares sense, and had to change it to use SVD.