Rectangular matrix support for Full Pivot LU decomposition

[ghost]

Currently, both LU decomposition algorithms require the input matrix to be square. However, neither Matlab nor Eigen require square matrices for their full-pivot LU decompositions, and Golub and Van Loan seem to (somewhat obliquely) indicate that full pivoting works for non-square matrices.

There's some discussion about numerical (in)stability of partial pivoting with rectangular matrices in the Eigen docs, but that seems to be contradicted by Golub and Van Loan (see PR #160). Also, Matlab seems to only have partial-pivot LU for dense matrices (at least, it only returns a row-interchange permutation matrix: "[L,U,P] = lu(A) returns an upper triangular matrix in U, a lower triangular matrix L with a unit diagonal, and a permutation matrix P, such that LU = PA"), but Matlab supports rectangular LU decomposition. Given the difference, I would think that you all would prefer err on the side of caution and restrict rectangular support to just FullPivLu, but I'm happy to look into it more!

As far as I can tell, the actual work should mostly involve removing m.rows() == m.cols() assertions from the decompose methods of FullPivLu, and adding them to functions where only square matrices make sense, like inverse and det.

[Andlon]

Thanks for writing up this issue!

Huh, it's very strange to me that Matlab does not have full pivoting support for dense LU. In any case, I agree with your proposal that we only implement rectangular matrix support for FullPivLu. Perhaps we could even change the panic message for PartialPivLu so that it tells the user to use FullPivLu if they try to use PartialPivLu with a rectangular matrix?

If you do come across anything more conclusive that suggests that partial pivoting is stable for rectangular matrix, I'd be happy to also support this for PartialPivLu!

[AtheMathmo]

Thanks for writing this up clearly.

I agree with your assessment. Let's stick to FullPivLU for now with an appropriate error message on PartialPivLu as @Andlon suggests.

If you do end up looking into this more thoroughly and believe that we can support this for PartialPivLU then please feel free to say so. I don't want to dictate the development of rulinalg via other existing libraries - though of course they provide a strong signal that we should not ignore!