Opened 13 months ago
Last modified 2 months ago
#29729 new enhancement
improve solve_right for some inexact rings
Reported by:  ghmwageringel  Owned by:  

Priority:  major  Milestone:  sage9.4 
Component:  linear algebra  Keywords:  
Cc:  AlexGhitza, nbruin  Merged in:  
Authors:  Reviewers:  
Report Upstream:  N/A  Work issues:  
Branch:  u/mjo/ticket/29729 (Commits, GitHub, GitLab)  Commit:  22efbc6a930103d86deb53a6e75cc15bac3cc2b2 
Dependencies:  Stopgaps: 
Description
Since #12406, the check
keyword to solve_right
is ignored over all inexact rings. As discussed on sagedevel, 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, [[0]]).solve_right(vector(SR, [1])) # 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, [[0]]).solve_right(vector(RIF, [1])) # this is not (0)
Change History (12)
comment:1 Changed 13 months ago by
 Branch set to u/mjo/ticket/29729
 Cc AlexGhitza added
 Commit set to adeb978ba1166eddd2dd22191fb7bcae694450b4
comment:2 followup: ↓ 4 Changed 13 months ago by
The SR proofofconcept 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 wikipediadeep.
comment:3 Changed 13 months ago by
 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 strictequality check.
comment:4 in reply to: ↑ 2 ; followup: ↓ 7 Changed 13 months ago by
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) ofb
. 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 wikipediadeep.
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
 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
That last commit should fix the square/nonsingular issue but probably needs some polish.
comment:7 in reply to: ↑ 4 ; followup: ↓ 8 Changed 13 months ago by
Replying to ghmwageringel:
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 ofCBF
.
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 ofA * 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 ; followup: ↓ 9 Changed 13 months ago by
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:
∀ 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∈ℝ^{n}  ∃ A∈A ∃ x∈x : Ax = y }.
(In particular, this means that ∀ A∈A ∀ x∈x ∃ y∈y such that Ax = y.)
Since A is nonempty, 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.
comment:9 in reply to: ↑ 8 Changed 13 months ago by
Replying to ghmwageringel:
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.
comment:10 Changed 13 months ago by
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
 Milestone changed from sage9.2 to sage9.3
comment:12 Changed 2 months ago by
 Milestone changed from sage9.3 to sage9.4
Sage development has entered the release candidate phase for 9.3. Setting a new milestone for this ticket based on a cursory review.
Here's a proofofconcept (not well tested) for SR.
New commits:
Trac #29729: add special case to solve_right() for symbolic systems.