Sage: Ticket #4502: [with patch; positive review] numerical noise in matrix_double_dense on intel mac os X 10.5: inverting a singular matrix
https://trac.sagemath.org/ticket/4502
<p>
(This has only been reported on intel macs running 10.4 or 10.5.)
</p>
<p>
From <a class="ext-link" href="http://groups.google.com/group/sage-devel/browse_frm/thread/fae59a9a9a53b8c0"><span class="icon"></span>sage-devel</a>:
</p>
<pre class="wiki">sage: A = matrix(RDF,3,range(1,10));A
[1.0 2.0 3.0]
[4.0 5.0 6.0]
[7.0 8.0 9.0]
sage: A.determinant()
0.0
sage: -A
[-1.0 -2.0 -3.0]
[-4.0 -5.0 -6.0]
[-7.0 -8.0 -9.0]
sage: b = A.numpy(); b
array([[ 1., 2., 3.],
[ 4., 5., 6.],
[ 7., 8., 9.]])
sage: import scipy
sage: import scipy.linalg
sage: scipy.linalg.det(b)
0.0
sage: scipy.linalg.inv(b)
array([[ -4.50359963e+15, 9.00719925e+15, -4.50359963e+15],
[ 9.00719925e+15, -1.80143985e+16, 9.00719925e+15],
[ -4.50359963e+15, 9.00719925e+15, -4.50359963e+15]])
</pre><p>
This leads to a doctest failure for <code>matrix_double_dense.py</code>.
</p>
en-usSagehttps://trac.sagemath.org/chrome/site/logo_sagemath_trac.png
https://trac.sagemath.org/ticket/4502
Trac 1.1.6jasonFri, 14 Nov 2008 05:50:30 GMT
https://trac.sagemath.org/ticket/4502#comment:1
https://trac.sagemath.org/ticket/4502#comment:1
<p>
I think the easiest way to avoid this (until we update to the latest scipy and test this again) is to check the scipy.linalg.det explicitly before attempting doing an inverse. If that determinant is zero (or close enough to it), then return the same error as the doctest returns (i.e., some scipy error).
</p>
TicketGeorgSWeberSun, 16 Nov 2008 09:49:51 GMTattachment set
https://trac.sagemath.org/ticket/4502
https://trac.sagemath.org/ticket/4502
<ul>
<li><strong>attachment</strong>
set to <em>4502-crude_hack.patch</em>
</li>
</ul>
<p>
patch against Sage 3.2.rc1
</p>
TicketGeorgSWeberSun, 16 Nov 2008 09:54:01 GMTsummary changed
https://trac.sagemath.org/ticket/4502#comment:2
https://trac.sagemath.org/ticket/4502#comment:2
<ul>
<li><strong>summary</strong>
changed from <em>numerical noise in matrix_double_dense on intel mac os X 10.5: inverting a singular matrix</em> to <em>[with patch; needs review] numerical noise in matrix_double_dense on intel mac os X 10.5: inverting a singular matrix</em>
</li>
</ul>
<p>
Well, this makes the doctest failure go away on my Intel Mac. (I followed the proposal from Jason.)
</p>
<p>
But please someone have a look if it is worth the price paid.
</p>
TicketjhpalmieriSun, 16 Nov 2008 16:57:59 GMT
https://trac.sagemath.org/ticket/4502#comment:3
https://trac.sagemath.org/ticket/4502#comment:3
<p>
It fixes the doctest failure for me, too. One or two questions:
</p>
<p>
Could the error message
</p>
<pre class="wiki">raise ValueError, "self must be square"
</pre><p>
be changed to something like "matrix must be square"? Actually, why is this here? I get an error message about non-square matrices (I think produced by scipy) before applying the patch.
</p>
TicketmhamptonTue, 18 Nov 2008 16:46:58 GMTsummary changed
https://trac.sagemath.org/ticket/4502#comment:4
https://trac.sagemath.org/ticket/4502#comment:4
<ul>
<li><strong>summary</strong>
changed from <em>[with patch; needs review] numerical noise in matrix_double_dense on intel mac os X 10.5: inverting a singular matrix</em> to <em>[with patch; with review; needs work] numerical noise in matrix_double_dense on intel mac os X 10.5: inverting a singular matrix</em>
</li>
</ul>
<p>
I don't think this fix is worth the performance penalty, which for me looks like about a 25% hit. The determinant is not a great measure of singularity. I think it would be better to just change the doctest.
</p>
TicketjasonTue, 18 Nov 2008 16:58:53 GMTattachment set
https://trac.sagemath.org/ticket/4502
https://trac.sagemath.org/ticket/4502
<ul>
<li><strong>attachment</strong>
set to <em>trac_4502.patch</em>
</li>
</ul>
TicketjasonTue, 18 Nov 2008 16:59:41 GMTsummary changed
https://trac.sagemath.org/ticket/4502#comment:5
https://trac.sagemath.org/ticket/4502#comment:5
<ul>
<li><strong>summary</strong>
changed from <em>[with patch; with review; needs work] numerical noise in matrix_double_dense on intel mac os X 10.5: inverting a singular matrix</em> to <em>[with patch; needs review] numerical noise in matrix_double_dense on intel mac os X 10.5: inverting a singular matrix</em>
</li>
</ul>
<p>
doctest changed in trac_4502.patch. Apply instead of previous patch.
</p>
TicketmabshoffTue, 18 Nov 2008 17:01:47 GMTsummary changed
https://trac.sagemath.org/ticket/4502#comment:6
https://trac.sagemath.org/ticket/4502#comment:6
<ul>
<li><strong>summary</strong>
changed from <em>[with patch; needs review] numerical noise in matrix_double_dense on intel mac os X 10.5: inverting a singular matrix</em> to <em>[with patch; positive review] numerical noise in matrix_double_dense on intel mac os X 10.5: inverting a singular matrix</em>
</li>
</ul>
<p>
Positive review on trac_4502.patch
</p>
<p>
Cheers,
</p>
<p>
Michael
</p>
TicketmabshoffTue, 18 Nov 2008 17:44:41 GMTsummary changed
https://trac.sagemath.org/ticket/4502#comment:7
https://trac.sagemath.org/ticket/4502#comment:7
<ul>
<li><strong>summary</strong>
changed from <em>[with patch; positive review] numerical noise in matrix_double_dense on intel mac os X 10.5: inverting a singular matrix</em> to <em>[with patch; needs work] numerical noise in matrix_double_dense on intel mac os X 10.5: inverting a singular matrix</em>
</li>
</ul>
<p>
Sorry, I spoke too soon. This needs work since it does not apply as is.
</p>
<p>
Cheers,
</p>
<p>
Michael
</p>
TicketGeorgSWeberTue, 18 Nov 2008 21:23:33 GMTattachment set
https://trac.sagemath.org/ticket/4502
https://trac.sagemath.org/ticket/4502
<ul>
<li><strong>attachment</strong>
set to <em>4502-next_try.patch</em>
</li>
</ul>
<p>
patch against Sage 3.2.rc1; apply only this one
</p>
TicketGeorgSWeberTue, 18 Nov 2008 21:26:40 GMTsummary changed
https://trac.sagemath.org/ticket/4502#comment:8
https://trac.sagemath.org/ticket/4502#comment:8
<ul>
<li><strong>summary</strong>
changed from <em>[with patch; needs work] numerical noise in matrix_double_dense on intel mac os X 10.5: inverting a singular matrix</em> to <em>[with patch; needs review] numerical noise in matrix_double_dense on intel mac os X 10.5: inverting a singular matrix</em>
</li>
</ul>
<p>
Just changes the failing doctest (and only that single one) to be #random, and should apply cleanly against Sage 3.2.rc1. Now on my box, the doctest for this file passes.
</p>
TicketmabshoffTue, 18 Nov 2008 21:28:54 GMTsummary changed
https://trac.sagemath.org/ticket/4502#comment:9
https://trac.sagemath.org/ticket/4502#comment:9
<ul>
<li><strong>summary</strong>
changed from <em>[with patch; needs review] numerical noise in matrix_double_dense on intel mac os X 10.5: inverting a singular matrix</em> to <em>[with patch; needs work] numerical noise in matrix_double_dense on intel mac os X 10.5: inverting a singular matrix</em>
</li>
</ul>
<p>
<code>#random</code> does not fix it since the exception still gets thrown. I tried it, so that is how I found out. In IRC we concluded to nuke the doctests for now and deal with that once <a class="missing wiki">SciPy?</a> 0.7 is out.
</p>
<p>
Cheers,
</p>
<p>
Michael
</p>
TicketGeorgSWeberTue, 18 Nov 2008 21:43:00 GMT
https://trac.sagemath.org/ticket/4502#comment:10
https://trac.sagemath.org/ticket/4502#comment:10
<p>
Hi John,
</p>
<blockquote class="citation">
<p>
Could the error message
</p>
</blockquote>
<pre class="wiki">raise ValueError, "self must be square"
</pre><blockquote class="citation">
<p>
be changed to something like "matrix must be square"? Actually, why is this here? I get an error message about non-square matrices (I think produced by scipy) before applying the patch.
</p>
</blockquote>
<p>
This test and error message were copied verbatim from lines 759/760 of that same file (function "def determinant(self)"). In my first patch, I introduced the additional calculation of the determinant in the function for inverting the matrix, and I wanted to make sure that all necessary sanity checks were in place.
</p>
<p>
Would it be possible for you to provide not only a concrete wording, but a new trac ticket with patch, if you like some other wording better for that error message? I'll review it then. (Maybe just removing this test from the .pyx file and relying on scipy producing the error message for us also for determinants on non-square matrices is an option here.)
</p>
<p>
Cheers,<br />
gsw
</p>
TicketGeorgSWeberTue, 18 Nov 2008 21:48:49 GMT
https://trac.sagemath.org/ticket/4502#comment:11
https://trac.sagemath.org/ticket/4502#comment:11
<blockquote class="citation">
<p>
#random does not fix it since the exception still gets thrown. I tried it, so that is how I found out.
</p>
</blockquote>
<p>
Strange. So does this mean the doctest passing as follows:
</p>
<pre class="wiki">susanne-webers-computer:~/Public/sage/sage-3.2.rc1/devel/sage georgweber$ sage -t sage/matrix/matrix_double_dense.pyx
sage -t devel/sage-test2/sage/matrix/matrix_double_dense.pyx
[4.2 s]
----------------------------------------------------------------------
All tests passed!
Total time for all tests: 4.2 seconds
</pre><p>
is not enough here? I'm a bit puzzled, which exception exactly gets "thrown nevertheless" and when? I'm fine with nuking the doctest and postponing the issue till scipy 0.7 is there and included, I ask just out of curiosity.
</p>
TicketmabshoffTue, 18 Nov 2008 21:52:42 GMT
https://trac.sagemath.org/ticket/4502#comment:12
https://trac.sagemath.org/ticket/4502#comment:12
<p>
With
</p>
<pre class="wiki">diff -r 50c6651a1286 sage/matrix/matrix_double_dense.pyx
--- a/sage/matrix/matrix_double_dense.pyx Tue Nov 18 08:58:27 2008 -0800
+++ b/sage/matrix/matrix_double_dense.pyx Tue Nov 18 13:49:26 2008 -0800
@@ -441,7 +441,7 @@
sage: A.determinant() < 10e-12
True
- sage: ~A
+ sage: ~A # random
Traceback (most recent call last):
...
LinAlgError: singular matrix
</pre><p>
I get
</p>
<pre class="wiki">sage -t devel/sage/sage/matrix/matrix_double_dense.pyx
**********************************************************************
File "/scratch/mabshoff/release-cycle/sage-3.2.rc2/devel/sage/sage/matrix/matrix_double_dense.pyx", line 444:
sage: print "ignore this"; ~A # random
Exception raised:
Traceback (most recent call last):
File "/scratch/mabshoff/release-cycle/sage-3.2.rc2/local/bin/ncadoctest.py", line 1231, in run_one_test
self.run_one_example(test, example, filename, compileflags)
File "/scratch/mabshoff/release-cycle/sage-3.2.rc2/local/bin/sagedoctest.py", line 38, in run_one_example
OrigDocTestRunner.run_one_example(self, test, example, filename, compileflags)
File "/scratch/mabshoff/release-cycle/sage-3.2.rc2/local/bin/ncadoctest.py", line 1172, in run_one_example
compileflags, 1) in test.globs
File "<doctest __main__.example_15[10]>", line 1, in <module>
print "ignore this"; ~A # random###line 444:
sage: print "ignore this"; ~A # random
File "matrix_double_dense.pyx", line 460, in sage.matrix.matrix_double_dense.Matrix_double_dense.__invert__ (sage/matrix/matrix_double_dense.c:3878)
File "/scratch/mabshoff/release-cycle/sage-3.2.rc2/local/lib/python2.5/site-packages/scipy/linalg/basic.py", line 306, in inv
if info>0: raise LinAlgError, "singular matrix"
LinAlgError: singular matrix
**********************************************************************
1 items had failures:
1 of 11 in __main__.example_15
***Test Failed*** 1 failures.
For whitespace errors, see the file /scratch/mabshoff/release-cycle/sage-3.2.rc2/tmp/.doctest_matrix_double_dense.py
[2.4 s]
exit code: 1024
</pre>
TicketGeorgSWeberTue, 18 Nov 2008 21:57:36 GMT
https://trac.sagemath.org/ticket/4502#comment:13
https://trac.sagemath.org/ticket/4502#comment:13
<p>
Ahhh, sorry for the noise --- the exception does not fly on my box because my box believes the matrix would be invertible and provides an inverse matrix --- but on other computers the matrix is not invertible and there, the exception comes up.
</p>
<p>
Would that not be an issue to enhance the "#random" pragma then?
</p>
<p>
More precisely: Allowing on certain boxes the calculation to finish (but throw away the result), and on other boxes having an exception (but again throwing the outcome)?
</p>
<p>
Just point out to me the respective parts of the IRC conversation if all that info is in there ... and if you find the time and have patience to do so ... ;-)
</p>
TicketmabshoffTue, 18 Nov 2008 22:00:55 GMT
https://trac.sagemath.org/ticket/4502#comment:14
https://trac.sagemath.org/ticket/4502#comment:14
<p>
The issue is the following:
</p>
<pre class="wiki">[10:34am] wstein|afk: mabshoff -- all that #random does is the following.
[10:35am] wstein|afk: sage: foo # random
[10:35am] wstein|afk: gets replaced by
[10:35am] wstein|afk: sage: _ = foo
[10:35am] wstein|afk: That's it.
[10:35am] mabshoff: mhh
[10:35am] wstein|afk: So it won't work at all with exceptions.
[10:35am] wstein|afk: Better might be:
[10:35am] wstein|afk: sage: foo # random
[10:35am] wstein|afk: gets replaced by
[10:35am] wstein|afk: sage: foo
[10:35am] wstein|afk: ...
[10:35am] wstein|afk: If that worked.
[10:35am] wstein|afk: I didn't know about ... when I wrote # random.
</pre><p>
The goal now is to get 3.2 out and the easiest way to do this is to nuke the affected doctests and then deal with <code>#random</code> properly in 3.2.1.
</p>
<p>
Cheers,
</p>
<p>
Michael
</p>
TicketGeorgSWeberTue, 18 Nov 2008 22:17:05 GMT
https://trac.sagemath.org/ticket/4502#comment:15
https://trac.sagemath.org/ticket/4502#comment:15
<p>
Thanks a lot!!!
</p>
TicketGeorgSWeberTue, 18 Nov 2008 22:17:47 GMTattachment set
https://trac.sagemath.org/ticket/4502
https://trac.sagemath.org/ticket/4502
<ul>
<li><strong>attachment</strong>
set to <em>4502-next_try.2.patch</em>
</li>
</ul>
<p>
apply only this one (doctest just nuked)
</p>
TicketGeorgSWeberTue, 18 Nov 2008 22:19:15 GMTsummary changed
https://trac.sagemath.org/ticket/4502#comment:16
https://trac.sagemath.org/ticket/4502#comment:16
<ul>
<li><strong>summary</strong>
changed from <em>[with patch; needs work] numerical noise in matrix_double_dense on intel mac os X 10.5: inverting a singular matrix</em> to <em>[with new patch; needs review] numerical noise in matrix_double_dense on intel mac os X 10.5: inverting a singular matrix</em>
</li>
</ul>
<p>
Is it that what you mean by "nuke"ing?
</p>
TicketmabshoffTue, 18 Nov 2008 22:26:01 GMT
https://trac.sagemath.org/ticket/4502#comment:17
https://trac.sagemath.org/ticket/4502#comment:17
<p>
Yes, but two things:
</p>
<ul><li>a comment would be nice that the issue needs to be fixed with a reference to the ticket
</li><li>The other test right before, i.e. "sage: A.inverse().det()" also needs to be commented out
</li></ul><p>
Once the above two happen with a new patch I will give this a positive review.
</p>
<p>
Cheers,
</p>
<p>
Michael
</p>
TicketGeorgSWeberWed, 19 Nov 2008 07:33:24 GMT
https://trac.sagemath.org/ticket/4502#comment:18
https://trac.sagemath.org/ticket/4502#comment:18
<p>
OK, no problem. (Except for the time it will take --- another 12 hours or so till I am back home.)
</p>
TicketGeorgSWeberWed, 19 Nov 2008 18:07:04 GMTattachment set
https://trac.sagemath.org/ticket/4502
https://trac.sagemath.org/ticket/4502
<ul>
<li><strong>attachment</strong>
set to <em>4502-next_try.3.patch</em>
</li>
</ul>
<p>
addressing the reviewers wishes (hopefully)
</p>
TicketGeorgSWeberWed, 19 Nov 2008 18:08:08 GMTsummary changed
https://trac.sagemath.org/ticket/4502#comment:19
https://trac.sagemath.org/ticket/4502#comment:19
<ul>
<li><strong>summary</strong>
changed from <em>[with new patch; needs review] numerical noise in matrix_double_dense on intel mac os X 10.5: inverting a singular matrix</em> to <em>[with even newer patch; needs review] numerical noise in matrix_double_dense on intel mac os X 10.5: inverting a singular matrix</em>
</li>
</ul>
TicketmabshoffWed, 19 Nov 2008 19:43:55 GMTsummary changed
https://trac.sagemath.org/ticket/4502#comment:20
https://trac.sagemath.org/ticket/4502#comment:20
<ul>
<li><strong>summary</strong>
changed from <em>[with even newer patch; needs review] numerical noise in matrix_double_dense on intel mac os X 10.5: inverting a singular matrix</em> to <em>[with patch; positive review] numerical noise in matrix_double_dense on intel mac os X 10.5: inverting a singular matrix</em>
</li>
</ul>
<p>
Positive review for 4502-next_try.3.patch. Thanks.
</p>
<p>
Cheers,
</p>
<p>
Michael
</p>
TicketmabshoffWed, 19 Nov 2008 19:44:24 GMTstatus changed; resolution set
https://trac.sagemath.org/ticket/4502#comment:21
https://trac.sagemath.org/ticket/4502#comment:21
<ul>
<li><strong>status</strong>
changed from <em>new</em> to <em>closed</em>
</li>
<li><strong>resolution</strong>
set to <em>fixed</em>
</li>
</ul>
<p>
Merged in Sage 3.2.rc2
</p>
Ticket