Opened 11 years ago

Closed 10 years ago

Last modified 10 years ago

#7852 closed defect (fixed)

solve_left for RDF matrices is WRONG — at Version 16

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 Reviewers: Martin Raum, Leif Leonhardy
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
            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


  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

to the Sage library.

Change History (20)

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.


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.

Note: See TracTickets for help on using tickets.