bug in matrix solve over inexact coefficients?

[AlexGhitza]

Hi,

I'm observing the following with version 9.1 (but not with 9.0 where the behavior is correct):

sage: m = matrix(SR, [0])

sage: b = vector([1])

sage: m.solve_right(b)

(0)

This should of course raise

ValueError: matrix equation has no solutions

It seems to be specific to inexact coefficients (the same problem occurs for RR and CC, but not for QQ).

Best,

Alex

[Dave Morris]

This is the expected behavior after trac #12406. The OUTPUT section of the docstring of solve_right says "... If the system is not square or does not have full rank, then a solution is attempted via other means. For example, over "RDF" or "CDF" a least-squares solution is returned..."

This makes sense, since testing for the existence of a solution over inexact rings is problematic, but I do think it would be better if the method could have some some reasonable test available, even if it is not the default.

[Michael Orlitzky]

This was changed to "do what MATLAB does" because of the numerical
issues inherent to inexact rings. While

m = matrix(SR, [0])

is singular and the system `m*x == [1]` has no solutions, with

m = matrix(SR, [0.0000000000000000000001]),

the system is easily solvable and we can't tell the difference between
the two. I think the real surprise here is that SR is an inexact ring.

[Nils Bruin]

I think this might need some work:

S=RealBallField(100)

M=Matrix(S,2,1,[1,1])
M.solve_right(vector([1,2]))

There's enough information here to conclude there is no solution; or in a rather deranged way, perhaps it should give a rather large ball back so that the multiplication results in a vector consisting of balls that actually contain the values 1 and 2. Currently, it just returns a wrong answer.

For systems that are even generically overdetermined, I don't think it makes sense to default to a least squares solution in a routine that is supposed to detect inconsistent systems. There the stable numerical behaviour actually IS being inconsistent.

[Michael Orlitzky]

On 5/22/20 3:39 PM, Nils Bruin wrote:
> I think this might need some work:
>
> S=RealBallField(100)
> M=Matrix(S,2,1,[1,1])
> M.solve_right(vector([1,2]))
>
> There's enough information here to conclude there is no solution; or in
> a rather deranged way, perhaps it should give a rather large ball back
> so that the multiplication results in a vector consisting of balls that
> actually contain the values 1 and 2. Currently, it just returns a wrong
> answer.

This "accidentally" gave the right answer (no solution) in the past due
to the "check" parameter, implemented in the base matrix class,

# Have to check that we actually solved the equation.
if self*X != B:
raise ValueError("matrix equation has no solutions")
return X

that is generally nonsensical over an inexact ring (it rejects good
answers). The change you're seeing -- and the change in the original
post -- is that in the base matrix class we no longer try to check the
answer using equality over inexact rings. While removing the check may
have toggled a few (right, wrong) answers, I believe that doing so was
correct in and of itself.

> For systems that are even generically overdetermined, I don't think it
> makes sense to default to a least squares solution in a routine that is
> supposed to detect inconsistent systems. There the stable numerical
> behaviour actually IS being inconsistent.

How should we know if an inexact system is overdetermined? (Does the
term "generically overdetermined" have some technical meaning?)

We're only doing least squares over RDF and CDF. The default solve is
otherwise still totally naive. The real benefit of the refactoring is
that now we are in a position to override _solve_right_general() and
_solve_right_nonsingular_square() in any matrix subclasses where
additional information can be taken advantage of. For example, if the
answer can be checked reliably over RealBallField(100), then it should
be fairly easy to do so without wrecking things for, say, RR.

[AlexGhitza]

Hi,

On Saturday, May 23, 2020 at 3:55:06 AM UTC+10, Michael Orlitzky wrote:

This was changed to "do what MATLAB does" because of the numerical
issues inherent to inexact rings. While

  m = matrix(SR, [0])

is singular and the system `m*x == [1]` has no solutions, with

  m = matrix(SR, [0.0000000000000000000001]),

the system is easily solvable and we can't tell the difference between
the two. I think the real surprise here is that SR is an inexact ring.

It's not surprising that SR is inexact, since it has embedded inexact rings:

sage: SR.has_coerce_map_from(RR)

True

sage: SR.has_coerce_map_from(RDF)

True

However, it has a mechanism for detecting whether a given element is exact:

sage: SR(0).is_exact()

True

sage: SR(0.0).is_exact()

False

So one could envision checking whether the coefficient matrix of the system consists solely of exact elements, in which case Sage could do the right thing as opposed to "whatever MATLAB does".

I would also argue that, despite the validity of the arguments regarding inexact rings, this is a change in behavior that would have benefited from a deprecation warning for a short while.

Another user-friendlier way of doing this (making the check argument irrelevant in the inexact case) would have been to display a warning when the user asks for check=True in the inexact case, rather than simply ignoring it.

Best,

Alex

[Michael Orlitzky]

On 5/22/20 6:40 PM, AlexGhitza wrote:
> I would also argue that, despite the validity of the arguments regarding
> inexact rings, this is a change in behavior that would have benefited
> from a deprecation warning for a short while.

We were pretty careful not to break anything in the sage library. For
typical inexact rings, the risk of breakage was minimal, because not
much worked to begin with. What we overlooked was that not all inexact
rings are typical, and that in cases like RealBallField and exact-SR,
the change could turn a right answer to a wrong one *in a situation that
wasn't ridiculous to begin with.* None of those cases were doctested, so
they were easy to overlook.

> Another user-friendlier way of doing this (making the check argument
> irrelevant in the inexact case) would have been to display a warning
> when the user asks for check=True in the inexact case, rather than
> simply ignoring it.

The base class defaults to check=True, and that can't be changed, so you
would get a warning on every inexact solve in that case. Moreover for
the "usual" inexact rings, things were unilaterally improved and not
deserving of a warning.

It will probably take a bit of domain-specific knowledge to make RBF/CBF
work as Nils expects, but reverting SR back to the old behavior when all
of its elements are exact shouldn't be hard if it doesn't cause any new
unforeseen problems. Let's just fix it instead of throwing warnings.

[Dave Morris]

From the documentation at https://www.mathworks.com/help/matlab/ref/mldivide.html, it appears to me that MATLAB gives a warning: "Warning: Matrix is close to singular or badly scaled. Results may be inaccurate."  That seems to me to be a better default behavior than what sage is doing now, but I think it would make sense to let users choose to turn off the warning (and maybe also let them upgrade it to an error).

[Markus Wageringel]

Indeed, instead of ignoring the `check` keyword altogether over inexact fields, it would be useful to verify the solution when that is possible, such as for exact input over SR and for ball and intervall fields. I have opened #29729 for this. For general inexact rings though, the result cannot be verified, so it would not make sense to do the check by default.

Am Samstag, 23. Mai 2020 07:15:45 UTC+2 schrieb Dave Morris:

From the documentation at https://www.mathworks.com/help/matlab/ref/mldivide.html, it appears to me that MATLAB gives a warning: "Warning: Matrix is close to singular or badly scaled. Results may be inaccurate."  That seems to me to be a better default behavior than what sage is doing now, but I think it would make sense to let users choose to turn off the warning (and maybe also let them upgrade it to an error).

 The generic implementation of `solve_right` does not provide this. I would suggest to use a type that has a specialized implementation for this, such as RDF, CDF, CBF. For example over RDF and CDF, a warning is issued if the input is ill-conditioned and an error is raised if the matrix is singular, which is a reasonable behavior:

sage: matrix(RDF, [[2/3, 1], [0.4, 0.6]]).solve_right(vector(RDF, [1, 1]))
/Applications/SageMath/src/bin/sage-ipython:1: LinAlgWarning: Ill-conditioned matrix (rcond=2.77556e-17): result may not be accurate.
  #!/usr/bin/env sage-python
(5404319552844596.0, -3602879701896396.0)
sage: matrix(RDF, [[0]]).solve_right(vector(RDF, [1]))
...
LinAlgError: Matrix is singular.

[Dave Morris]

Although both behaviors may be reasonable, I don't like the fact that each line of the following code gives a completely different result depending on whether F = RR (returns a result that may be garbage, without warning) or F = RDF (gives an error, even though there may be a valid solution).  My own preference for default behavior (for F = RR, at least, which is where I think a naive user is likely to be working) is intermediate between these two, and is what I think MATLAB does: return a result but give a warning.

sage: matrix(F, [[0]]).solve_right(vector(F, [0]))
sage: matrix(F, [[0]]).solve_right(vector(F, [1]))

Along these lines, I note that the code for `Matrix_generic_dense.solve_right` has

try:
    X = self._solve_right_nonsingular_square(C, check_rank=True)
except NotFullRankError:
    X = self._solve_right_general(C, check=check)

It seems to me that the except clause would be a good place to put a warning message over inexact rings. However, I know nothing about numerical analysis (and not much about sage), so I am certainly not giving an expert opinion on anything in this discussion, just the viewpoint of this sage user, and I apologize if it is noise.

[Markus Wageringel]

Don't get me wrong. I would greatly prefer if RR would detect the ill-conditioned cases as well, but, as it stands, this is not implemented. Neither LU nor QR decompositions are supported for RR, so it is not easy to change this. The default implementation for solve_right just computes an echelonization.

Putting a warning in the except clause might be a good first step indeed, as it should at least detect the case of a matrix that is singular to machine precision.

For the record, this is what Matlab does in these cases:

>> linsolve([0], [1])
Warning: Matrix is singular to working precision.

ans =

Inf

>> linsolve([0], [0])
Warning: Matrix is singular to working precision.

ans =

NaN