#7852 closed defect (fixed)
solve_left for RDF matrices is WRONG
Reported by: | was | Owned by: | was |
---|---|---|---|
Priority: | blocker | Milestone: | sage-4.7.2 |
Component: | linear algebra | Keywords: | |
Cc: | jason | Merged in: | sage-4.7.2.alpha3 |
Authors: | Rob Beezer, Leif Leonhardy | Reviewers: | Martin Raum, Leif Leonhardy, Rob Beezer |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | #11848 | Stopgaps: |
Description (last modified by )
Observe the docstring for solve_left for an RDF matrix:
sage: A = random_matrix(RDF,3) sage: A.solve_left? Solve the equation A*x = b, where EXAMPLES: sage: A = matrix(RDF, 3,3, [1,2,5,7.6,2.3,1,1,2,-1]); A [ 1.0 2.0 5.0] [ 7.6 2.3 1.0] [ 1.0 2.0 -1.0] sage: b = vector(RDF,[1,2,3]) sage: x = A.solve_left(b); x (-0.113695090439, 1.39018087855, -0.333333333333) sage: A*x (1.0, 2.0, 3.0)
But that is solve_right.
This was evidently introduced by maybe Grout's "Switch the RDF and CDF matrices to a numpy 1.2 backend; factor out common functionality to matrix_double_dense.pyx.".
Reported by Stephanie Dietzel
Apply
- trac_7852-solve-linear-systems-CDF.patch
- trac_7852-fix_noise_errors_in_preparser_examples.reviewer.patch
- trac_7852-fix_noise_errors_in_polys.reviewer.patch
- trac_7852-fix_noisy_zero_error_in_matrix_double_dense.reviewer.patch
- trac_7852-adjust_noisy_zero_term_threshold_for_polys.reviewer.patch
to the Sage library.
Attachments (5)
Change History (25)
comment:1 Changed 11 years ago by
Changed 10 years ago by
comment:2 Changed 10 years ago by
- Status changed from new to needs_review
Patch corrects the naming problem, and creates upgraded solve_left
and
solve_right
.
This required just a few touch-ups where the name had been wrong, or the improved accuracy changed results.
comment:3 Changed 10 years ago by
comment:4 follow-up: ↓ 5 Changed 10 years ago by
Without the patch, "solve_right" seems to work for nonsquare matrices over RDF, I guess because it's using the generic "solve_right" method from matrix2.pyx. With the patch, it looks like this won't work any more; is that right? Is there a good reason for not allowing square matrices? We could instead add a method _solve_right_nonsingular_square
to matrix_double_dense.pyx, and then I think this should be called by the generic solve_right method.
For "solve_left", could we just remove the version from matrix_double_dense.pyx altogether, and just use the generic version?
comment:5 in reply to: ↑ 4 Changed 10 years ago by
Replying to jhpalmieri:
Yes, there is some tension in handling inexact matrices with routines designed for exact matrices (which is the root issue here, I believe). NumPy/SciPy
just does not do certain things, such as computing a rank.
I'm hoping to get some guidance from the NumPy
folks here at Sage Days 29 right now and we can chat some more.
comment:6 Changed 10 years ago by
- Status changed from needs_review to needs_work
1) You can use svd to first get a square matrix that you can use to solve the system of linear equations. 2) Why do you copy the code for left and right solve. There is a generic fallback in matrix2.pyx. If you want to replace the docstring for this class, can't you do just that?
comment:7 Changed 10 years ago by
We could use a direct call to SciPy
solve for square nonsingular and analyze/exploit the SVD for nonsingular/rectangular.
comment:8 Changed 10 years ago by
- Status changed from needs_work to needs_review
Finally getting back to this one.
- solve_left and solve_right are confused at the moment - I think this is important to fix.
- The code is different, in part because the error messages are explicit, naming rows and columns when sizes do not match.
- The longer I think about it, the more I think we should never be feeding floating-point matrices into exact routines. We should expose the ScyPy/NumPy/LAPACK routines as much as possible and not pretend the exact routines are going to give reasonable answers as a "fallback."
- I'm happy to try to get solutions to systems with non-square coefficient matrices with floating-point entries, but despite a lot of thought and reading, I'm not convinced of the right approach. In any event, I'd like to fix the problem at hand, and do the non-square somewhere else.
- I am going to set this to "needs review" but I'm happy to entertain more discussion - I'll just need convincing.
Rob
comment:9 Changed 10 years ago by
What I was talking about when I mentioned the fallback is a piece of generic code for solve_left that calls solve_right. Thus you only need to implement one version. And you should, because code duplication makes mistakes more likely and maintenance harder. I will change this when reviewing this, and then will wait for your approval.
comment:10 Changed 10 years ago by
- Reviewers set to Martin Raum
- Status changed from needs_review to positive_review
I didn't touch the implementation of solve_right. I changed my mind for the following simple reason. While the code is quite short and easy to check the documentation is very different (left replaced by right). So it seems quite impractical to duplicate it.
We should think about modifying the documentation so that it gets included in the reference manual. Also there is a pseudoregression in the code. Sage doesn't pretend it could solve equations with double dense matrices as right hand side. Scipy covers also this case, so we should open a new ticket for this.
comment:11 Changed 10 years ago by
- Milestone changed from sage-4.7.1 to sage-4.7.2
comment:12 Changed 10 years ago by
- Description modified (diff)
- Priority changed from major to blocker
comment:13 Changed 10 years ago by
- Merged in set to sage-4.7.2.alpha3
- Resolution set to fixed
- Status changed from positive_review to closed
comment:14 Changed 10 years ago by
- Dependencies set to #11848
Changed 10 years ago by
Reviewer patch. Apply on top of main patch, which causes a signed zero on Itanium 2.
Changed 10 years ago by
Reviewer patch. Apply on top of main patch, which causes doctests to fail on a couple of systems due to noisy zero terms.
comment:15 Changed 10 years ago by
- Description modified (diff)
- Reviewers changed from Martin Raum to Martin Raum, Leif Leonhardy
Attached two reviewer patches fixing doctest errors on some systems.
Third patch to fix doctest errors in sage/matrix/matrix_double_dense.pyx
coming soon.
The patch to the preparser examples requires #11848.
Changed 10 years ago by
Reviewer patch. Apply on top of main patch, which causes another doctest to fail on a couple of systems due to a noisy zero in a vector result.
comment:16 Changed 10 years ago by
- Description modified (diff)
Third reviewer patch (fixing a doctest error in matrix_double_dense.pyx
) is up.
comment:17 Changed 10 years ago by
- Reviewers changed from Martin Raum, Leif Leonhardy to Martin Raum, Leif Leonhardy, Rob Beezer
The three reviewer patches apply, build and pass their tests once #11848 is applied to a recent alpha3 prerelease that already has the main patch. Currently passes all tests in sage/matrix
but will run full tests overnight.
This has a positive review from me based on the reviewer contributions, but is marked as closed right now anyway. I'll only retreat if something unexpected happens overnight.
Leif - thank-you very much for getting these noisy tests cleaned up.
Rob
Changed 10 years ago by
Reviewer patch. Slightly adjust threshold for noisy zero terms in polynomials, needed on MacOS X 10.6 x86_64 with GCC 4.2.1. Apply on top of other patches (i.e., the other reviewer patch to polynomial_element.pyx
).
comment:18 Changed 10 years ago by
- Description modified (diff)
I've attached yet another patch slightly increasing one epsilon, since one doctest still failed on MacOS X 10.6 (due to almost zero terms in a polynomial).
comment:19 Changed 10 years ago by
Yes, I meant to bring this up in the flurry of noisy doctest fixes. I have taken to not cutting it too fine on the zeros. For example, if I get a "zero" on the order of 10^-16
, then I have taken to clipping zeros at 10^-14
. There seems to be too much variability across systems, so if you cut it too close, then some system is eventually going to complain, it seems.
The latest patch looks good and applies and builds find, passing all tests in sage/ring/polynomial
on the alpha3-prerelease I picked up about 10 days ago. So ready to go, again.
comment:20 Changed 10 years ago by
Follow-up ticket further increasing the threshold for noisy zero terms in polynomials (from 2.5e-15, needed for bsd.math
, to 2.7e-15) for MacOS X 10.4 PPC once more: #11901
This might be related to #4932. I'll put this on my list for the bug day coming up.