On Nov 12, 6:41 am, Jason Grout <[EMAIL PROTECTED]> wrote:
> John H Palmieri wrote:
> > an error in matrix_double_dense, which I couldn't find on trac (should
> > I create a ticket?):
>
> > sage -t devel/sage/sage/matrix/matrix_double_dense.pyx
> > **********************************************************************
> > File "/Applications/sage/devel/sage/sage/matrix/
> > matrix_double_dense.pyx", line 444:
> > sage: ~A
> > Expected:
> > Traceback (most recent call last):
> > ...
> > LinAlgError: singular matrix
> > Got:
> > [-4.50359962737e+15 9.00719925474e+15 -4.50359962737e+15]
> > [ 9.00719925474e+15 -1.80143985095e+16 9.00719925474e+15]
> > [-4.50359962737e+15 9.00719925474e+15 -4.50359962737e+15]
>
> Hmm. This code was just moved to a numpy backend. Scipy correctly saw
> this matrix as singular on my system. However, it looks like you are
> seeing a spurious inverse, most likely due to numerical issues.
>
> Could you post the output of the following commands:
>
> sage: A = matrix(RDF,3,range(1,10));A
> sage: A.determinant()
> sage: ~A
> sage: b=A.numpy(); b
> sage: import scipy
> sage: import scipy.linalg
> sage: scipy.linalg.det(b)
> sage: scipy.linalg.inv(b)
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]])
That last pair of results looks suspicious, doesn't it?
> > **********************************************************************
> > File "/Applications/sage/devel/sage/sage/matrix/
> > matrix_double_dense.pyx", line 909:
> > sage: max((U*S*V.transpose()-m).list())<1e-15 # check
> > Expected:
> > True
> > Got:
> > False
>
> Again, could you post the output of the following commands:
>
> sage: m = matrix(RDF,3,2,range(6)); m
> sage: U,S,V = m.SVD()
> sage: U*S*V.transpose() # random low order bits
sage: m = matrix(RDF,3,2,range(6)); m
[0.0 1.0]
[2.0 3.0]
[4.0 5.0]
sage: U,S,V = m.SVD()
sage: U*S*V.transpose() # random low order bits
[0.0 1.0]
[2.0 3.0]
[4.0 5.0]
max((U*S*V.transpose()-m).list())
1.7763568394e-15
John
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://www.sagemath.org
-~----------~----~----~----~------~----~------~--~---
