improve solve_right for some inexact rings

[mwageringel]

Since #12406, the check keyword to solve_right is ignored over all inexact rings. As discussed on sage-devel, this is a bit too drastic, as the solution can be checked for correctness in some cases.

Over SR, we can check if all entries are exact and then check whether the solution is correct. This was done in #33159.

Over ball and interval fields, we can check the error bounds to determine whether the result is valid.

sage: matrix(RIF, [[2/3, 1], [2/5, 3/5]]).solve_right(vector(RIF, [1, 1]))  # this solution is acceptable
([-infinity .. +infinity], [-infinity .. +infinity])
sage: matrix(RIF, [[0]]).solve_right(vector(RIF, [1]))  # this is not
(0)

CC: @aghitza @nbruin

Component: linear algebra

Branch/Commit: u/mjo/ticket/29729 @ 22efbc6

Issue created by migration from https://trac.sagemath.org/ticket/29729

[orlitzky]

comment:1

Here's a proof-of-concept (not well tested) for SR.

---

New commits:

adeb978

Trac #29729: add special case to solve_right() for symbolic systems.

[orlitzky]

Branch: u/mjo/ticket/29729

[orlitzky]

Commit: adeb978

[orlitzky]

comment:2

The SR proof-of-concept passes a ptestlong, at least.

As for the ball/interval fields: half of the problem with RBF is that we have a dedicated implementation of _solve_right_nonsingular_square for complex ball matrices, but not for RBF ones. Thus we miss out on the work that the arb library has already done for these systems. Compare:

sage: A = matrix(RBF, [[2/3, 1], [2/5, 3/5]])
sage: b = vector(RBF, [1, 1])
sage: A.solve_right(b)
(nan, nan)

sage: A = A.change_ring(CBF)
sage: b = b.change_ring(CBF)
sage: A.solve_right(b)
...
ValueError: unable to invert this matrix

(Note: we should catch that exception and turn it into a NotFullRankError).

For the generic interval solve/check, as a first approximation, maybe we could multiply A*x and check that each component thereof intersects the corresponding component (interval) of b. That's as opposed to the default "check" implementation using literal equality that would fail. This needs an expert opinion as my understanding of interval systems is only wikipedia-deep.

[orlitzky]

comment:3

It might also be a good idea to factor out the solution check into a separate method, that way we don't have to override all of _solve_right_general for interval matrices just to replace the strict-equality check.

[mwageringel]

comment:4

Replying to @orlitzky:

(Note: we should catch that exception and turn it into a NotFullRankError).

The NotFullRankError is currently only used internally for control flow, so we should only raise it if we want _solve_right_general to be tried when _solve_right_nonsingular_square fails, in the case of CBF.

For the generic interval solve/check, as a first approximation, maybe we could multiply A*x and check that each component thereof intersects the corresponding component (interval) of b. That's as opposed to the default "check" implementation using literal equality that would fail. This needs an expert opinion as my understanding of interval systems is only wikipedia-deep.

The literal equality check seemingly works with ball/interval arithmetics. Though, this is probably just a coincidence. If vector equality is implemented by checking whether not any entries are different (as opposed to all entries being equal), for ball fields this means it is checked that no entries are disjoint balls, hence they must intersect.

sage: A = matrix(RBF, [[2/3, 1], [2/5, 3/5]])
sage: b = vector(RBF, [1, 1])
sage: x = A.solve_right(b)
sage: A * x == b
True

However, it would be more correct to check that every entry of b, as a ball/interval, is fully contained in the corresponding entry of A * x:

sage: all(u in v for u, v in zip(b, A*x))
True

I am not very familiar with interval arithmetics either, though.

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

22efbc6

Trac #29729: new subclass for RealBallField matrices.

[sagetrac-git]

Changed commit from adeb978 to 22efbc6

[orlitzky]

comment:6

That last commit should fix the square/nonsingular issue but probably needs some polish.

[orlitzky]

comment:7

Replying to @mwageringel:

The NotFullRankError is currently only used internally for control flow, so we should only raise it if we want _solve_right_general to be tried when _solve_right_nonsingular_square fails, in the case of CBF.

If we can make check=True work, I think that's what we want. I'm imagining a square system with exact real entries (like powers of two) and then some extra redundant rows tacked onto the end. If the generic solver can eliminate those junk rows and if check=True can make sure we don't wind up with a wrong answer... well first we need to make check=True work, and then I can just try it.

The literal equality check seemingly works with ball/interval arithmetics. Though, this is probably just a coincidence...
However, it would be more correct to check that every entry of b, as a ball/interval, is fully contained in the corresponding entry of A * x:

From what I gather, a solution to an interval system Ax = b should be an interval vector x such that for all x in x, there exists some A in A and b in b with Ax = b. So I don't think either suggestion is quite right; mine because Ax can wind up being larger than necessary, and thus too lenient if we intersect it with b. Yours I would imagine is too strict: if everything is nonnegative, we should be able to pair up the left endpoints of x with the left endpoints of A and b and not insist that the right endpoints of A and b work (for the left endpoints of x) as well.

[mwageringel]

comment:8

Replying to @orlitzky:

From what I gather, a solution to an interval system Ax = b should be an interval vector x such that for all x in x, there exists some A in A and b in b with Ax = b.

This does not seem quite right to me, as then x could always be shrinkened to just a single point, thus reducing the error to 0.

In my naive understanding, as we do not know the exact values of A and b, the solution x of the system Ax = b should satisfy:

∀ *A*∈**A** ∀ *b*∈**b** ∃ *x*∈**x** such that *Ax* = *b*.

Given x and A, let me denote the product by y = Ax. Then y should contain the set

{ *y*∈ℝ<sup>n</sup> | ∃ *A*∈**A** ∃ *x*∈**x** : *Ax* = *y* }.

(In particular, this means that ∀ A∈A ∀ x∈x ∃ y∈y such that Ax = y.)

Since A is non-empty, we can fix an arbitrary A∈A. Now, given any b∈b, we find an x∈x satisfying Ax=b, and so b∈y.

Therefore, under the above assumptions, it would be a necessary condition that b is contained in y.

[mwageringel]

comment:9

Replying to @mwageringel:

In my naive understanding, as we do not know the exact values of A and b, the solution x of the system Ax = b should satisfy:

∀ *A*∈**A** ∀ *b*∈**b** ∃ *x*∈**x** such that *Ax* = *b*.

Come to think of it, this will never hold if A has more rows than columns (and not everything is exact), as then A and b can easily be chosen such that no solution exists.

[orlitzky]

comment:10

It looks like all of these definitions of a solution that we've been making up are valid in different contexts (see e.g. Interval linear systems as a necessary step in fuzzy linear systems by Lodwick and Dubois). However, since RBF is documented to be implemented with Arb and since square systems will be solved by Arb, we should probably use the same interpretation of "solution" that Arb uses. So far I've been unable to figure out what that is. I may have to ask on the mailing list.

[mkoeppe]

comment:12

Sage development has entered the release candidate phase for 9.3. Setting a new milestone for this ticket based on a cursory review.