Opened 13 months ago

Last modified 2 months ago

# improve solve_right for some inexact rings

Reported by: Owned by: gh-mwageringel major sage-9.4 linear algebra AlexGhitza, nbruin N/A u/mjo/ticket/29729 22efbc6a930103d86deb53a6e75cc15bac3cc2b2

### Description

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.

```sage: matrix(SR, []).solve_right(vector(SR, ))  # should raise an error
(0)
```

Over ball and intervall 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, []).solve_right(vector(RIF, ))  # this is not
(0)
```

### comment:1 Changed 13 months ago by mjo

• Branch set to u/mjo/ticket/29729
• Cc AlexGhitza added
• Commit set to adeb978ba1166eddd2dd22191fb7bcae694450b4

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.`

### comment:2 follow-up: ↓ 4 Changed 13 months ago by mjo

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.

### comment:3 Changed 13 months ago by mjo

• Cc nbruin added

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.

### comment:4 in reply to: ↑ 2 ; follow-up: ↓ 7 Changed 13 months ago by gh-mwageringel

Replying to mjo:

(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.

### comment:5 Changed 13 months ago by git

• Commit changed from adeb978ba1166eddd2dd22191fb7bcae694450b4 to 22efbc6a930103d86deb53a6e75cc15bac3cc2b2

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

 ​22efbc6 `Trac #29729: new subclass for RealBallField matrices.`

### comment:6 Changed 13 months ago by mjo

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

### comment:7 in reply to: ↑ 4 ; follow-up: ↓ 8 Changed 13 months ago by mjo

Replying to gh-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.

### comment:8 in reply to: ↑ 7 ; follow-up: ↓ 9 Changed 13 months ago by gh-mwageringel

Replying to mjo:

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:

AAbbxx such that Ax = b.

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

{ y∈ℝn | ∃ AAxx : Ax = y }.

(In particular, this means that ∀ AAxxyy such that Ax = y.)

Since A is non-empty, we can fix an arbitrary AA. Now, given any bb, we find an xx satisfying Ax=b, and so by.

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

### comment:9 in reply to: ↑ 8 Changed 13 months ago by gh-mwageringel

Replying to gh-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:

AAbbxx 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.

### comment:10 Changed 13 months ago by mjo

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.

### comment:11 Changed 8 months ago by mkoeppe

• Milestone changed from sage-9.2 to sage-9.3

### comment:12 Changed 2 months ago by mkoeppe

• Milestone changed from sage-9.3 to sage-9.4

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

Note: See TracTickets for help on using tickets.