Opened 11 years ago

Closed 10 years ago

Last modified 10 years ago

#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:

Status badges

Description (last modified by leif)

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

  1. trac_7852-solve-linear-systems-CDF.patch
  2. trac_7852-fix_noise_errors_in_preparser_examples.reviewer.patch
  3. trac_7852-fix_noise_errors_in_polys.reviewer.patch
  4. trac_7852-fix_noisy_zero_error_in_matrix_double_dense.reviewer.patch
  5. trac_7852-adjust_noisy_zero_term_threshold_for_polys.reviewer.patch

to the Sage library.

Attachments (5)

trac_7852-solve-linear-systems-CDF.patch (14.7 KB) - added by rbeezer 10 years ago.
trac_7852-fix_noise_errors_in_preparser_examples.reviewer.patch (1.1 KB) - added by leif 10 years ago.
Reviewer patch. Apply on top of main patch, which causes a signed zero on Itanium 2.
trac_7852-fix_noise_errors_in_polys.reviewer.patch (1.6 KB) - added by leif 10 years ago.
Reviewer patch. Apply on top of main patch, which causes doctests to fail on a couple of systems due to noisy zero terms.
trac_7852-fix_noisy_zero_error_in_matrix_double_dense.reviewer.patch (846 bytes) - added by leif 10 years ago.
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.
trac_7852-adjust_noisy_zero_term_threshold_for_polys.reviewer.patch (978 bytes) - added by leif 10 years ago.
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).

Download all attachments as: .zip

Change History (25)

comment:1 Changed 11 years ago by jason

This might be related to #4932. I'll put this on my list for the bug day coming up.

Changed 10 years ago by rbeezer

comment:2 Changed 10 years ago by rbeezer

  • 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 rbeezer

  • Authors set to Rob Beezer

comment:4 follow-up: Changed 10 years ago by jhpalmieri

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 rbeezer

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 mraum

  • 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 rbeezer

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 rbeezer

  • Status changed from needs_work to needs_review

Finally getting back to this one.

  1. solve_left and solve_right are confused at the moment - I think this is important to fix.
  1. The code is different, in part because the error messages are explicit, naming rows and columns when sizes do not match.
  1. 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."
  1. 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.
  1. 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 mraum

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 mraum

  • 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 jdemeyer

  • Milestone changed from sage-4.7.1 to sage-4.7.2

comment:12 Changed 10 years ago by leif

  • Description modified (diff)
  • Priority changed from major to blocker

comment:13 Changed 10 years ago by leif

  • 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 leif

  • Dependencies set to #11848

Changed 10 years ago by leif

Reviewer patch. Apply on top of main patch, which causes a signed zero on Itanium 2.

Changed 10 years ago by leif

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 leif

  • 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 leif

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 leif

  • Description modified (diff)

Third reviewer patch (fixing a doctest error in matrix_double_dense.pyx) is up.

comment:17 Changed 10 years ago by rbeezer

  • Authors changed from Rob Beezer to Rob Beezer, Leif Leonhardy
  • 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 leif

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 leif

  • 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 rbeezer

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 leif

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

Note: See TracTickets for help on using tickets.