Opened 12 years ago

Closed 10 years ago

# solve_left for RDF matrices is WRONG

Reported by: Owned by: was was blocker sage-4.7.2 linear algebra jason sage-4.7.2.alpha3 Rob Beezer, Leif Leonhardy Martin Raum, Leif Leonhardy, Rob Beezer N/A #11848

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

to the Sage library.

### comment:1 Changed 12 years ago by jason

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

### comment:2 Changed 11 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 11 years ago by rbeezer

• Authors set to Rob Beezer

### comment:4 follow-up: ↓ 5 Changed 11 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 11 years ago by rbeezer

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