What should happen when we solve a linear system?

[Michael Orlitzky]

Our documentation for Matrix.solve_right is schizophrenic, and it's
contagious. I'm trying to review #12406 that fixes some old issues with
the solution of linear systems. The problem I have is that solve_right()
doesn't know what it's supposed to do. The documentation says,

If self is a matrix A, then this function returns a vector or matrix X
such that AX=B. If B is a vector then X is a vector and if B is a
matrix, then X is a matrix.

Note: In Sage one can also write A \ B for A.solve_right(B), i.e.,
Sage implements the "the MATLAB/Octave backslash operator.

Those two statements are already at odds, because that's not what the
MATLAB backslash operator does:

https://www.mathworks.com/help/matlab/ref/mldivide.html

At the most basic level, the backslash operator has two cases for square
and non-square matrices:

* non-square gets an approximate solution
* square assumes non-singularity and throws an error otherwise

The docs for solve_right continue...

INPUT:

...
check - bool (default: True) - if False and self is nonsquare, may
not raise an error message even if there is no solution. This is
faster but more dangerous.

Which again, doesn't quite make sense in the context of an approximate
solution, or indeed *any* solution over an inexact ring, where the
answer will be approximate.

Several of these unusual behaviors are doctested, for example in the
tour of linear algebra,

http://doc.sagemath.org/html/en/tutorial/tour_linalg.html

where the very first example expects a solution from a singular square
system:

sage: A = Matrix([[1,2,3],[3,2,1],[1,1,1]])
sage: Y = vector([0, -4, -1])
sage: A.is_square()
True
sage: A.is_singular()
True
sage: A.solve_right(Y)
(-2, 1, 0)

That creates a sticky situation: to fix solve_right() over RDF or other
inexact rings, we basically have to add a _another_ special case that
ignores the documentation for solve_right(), and that doesn't quite
agree with what happens in exact rings.

Do tl;dr, what should solve_right() actually do? Is there a consistent
description of it that works for both exact and inexact rings? Should it
act like the backslash operator, or should we stop saying that we
implement the backslash operator?

The change that requires the least change would be to document that
non-square systems get approximate solutions, and to do away with the
"check" parameter entirely. Or maybe document that "check" only works
over exact rings, and to set it to False at the top of solve_right when
the ring is exact. Then we can tweak the implementation to work like the
MATLAB operator, by e.g. making it an error to solve a singular square
system. But as I mentioned, that breaks some doctests...

We could also just stop saying that we implement the backslash operator,
and put the approximate solvers in different methods. That's probably
cleaner and simpler, but it too will break a lot of stuff. I'm sure
there are dead-tree books out there that use the "\" exclusively.

[Markus Wageringel]

Over inexact rings, it is generally impossible to determine the rank or to check whether a square matrix is singular, so such a check should not be done at all, IMO. This implies that the `check` keyword should be ignored. This is also what https://trac.sagemath.org/ticket/13932 is about. For inexact rings like RDF/CDF, the Matlab behavior is consistent and predictable. In particular, the answer does not depend on numerical noise. In other words, for near-singular square matrices, the method should _not_ attempt to find a least-squares solution, but continue the computation as if the matrix was regular (which may lead to a warning or error). Otherwise, the output of the method would just be very unpredictable, so one would always have to check which type of solution it decided to return.

For inexact rings other than RDF/CDF, other notions of "approximate solution" may be more appropriate, but that does not seem to be a problem as long as it clearly documented.

For exact rings, it is possible to determine the rank and whether an exact solution exists. Therefore, I would expect the method to always return an exact solution if it exists and raise an error otherwise, regardless of whether the matrix is square or not. The `check` keyword is, in my understanding, a detail that helps speed up a computation if one already knows a priori that what it would check is unnecessary. Otherwise, the output may just be wrong.

I will try to reply to your changes on the ticket soon.