Ask Your Question
0

Set the precision of imported methods

asked 2013-11-21 11:01:49 +0100

gundamlh gravatar image

updated 2013-11-21 11:16:33 +0100

How to set the precision of "numpy" methods?

e.g., for calculating the singular values of a matrix using numpy methods

sage: R = RealField(100)
sage: R
Real Field with 100 bits of precision
sage: A = matrix(R ,2,2, [1.746, 0.940, 1.246, 1.898])
sage: A
[ 1.7460000000000000000000000000 0.94000000000000000000000000000]
[ 1.2460000000000000000000000000  1.8980000000000000000000000000]
sage: A.parent()
Full MatrixSpace of 2 by 2 dense matrices over Real Field 
with 100 bits of precision
sage: A = np.array(A) # the precision (100 bits) does not preserve in "numpy"
sage: U,sig,V = numpy.linalg.svd(A)
sage: sig
array([  1.08490731e+06,   1.97535694e+00])
# the precision (100 bits) does not preserve in "numpy", hey only 8 digits here

NumPy Data types

sage: np.array([1, 2, 3], dtype='f')

array([ 1., 2., 3.], dtype=float32)

That's a choice. However, I'd like some more convenient method like "numpy.set_digits(100)"...


Thanks in advance!

edit retag flag offensive close merge delete

1 Answer

Sort by ยป oldest newest most voted
2

answered 2013-11-22 16:48:21 +0100

tmonteil gravatar image

You can preserve the precision of your matrix while transferring it to numpy as follows:

sage: R = RealField(100)
sage: A = matrix(R ,2,2, [1.746, 0.940, 1.246, 1.898])
sage: import numpy as np
sage: a = np.array(A, dtype=object)
sage: a
array([[1.7460000000000000000000000000, 0.94000000000000000000000000000],
       [1.2460000000000000000000000000, 1.8980000000000000000000000000]], dtype=object)
sage: a[0]
array([1.7460000000000000000000000000, 0.94000000000000000000000000000], dtype=object)
sage: a[0][0]
1.7460000000000000000000000000
sage: a[0][0].parent()
Real Field with 100 bits of precision

Unfortunately, this will not solve your second problem:

sage: from scipy import linalg
sage: b = linalg.svd(a) ; b
(array([[-0.65092234, -0.75914432],
       [-0.75914432,  0.65092234]]),
 array([ 2.92405178,  0.73277362]),
 array([[-0.71216394, -0.70201319],
       [-0.70201319,  0.71216394]]))
sage: b[0][0][0]
-0.65092234401411786
sage: type(b[0][0][0])
<type 'numpy.float64'>

As you can see, the svd operation took your entries back to double precision floating point numbers. If you type:

sage: linalg.svd??

You will see that scipy uses lapack to solve the singular value decomposition, and the lapack library only work with simple or double precision floating-point numbers, not mpfr (arbitrary precision) numbers. For the same reason, Sage offers a SVD method for matrices over RDF and not over RR.

What may be hard to understand is why there is no generic method for doing that in any precision. Actually, it is not so easy since there are a lot of numerical stability issues when doing linear algebra over floating-point numbers, this is why Sage and scipy currently relies on a specialized library.

edit flag offensive delete link more

Comments

Thanks! > it is not so easy since there are a lot of numerical stability issues when doing linear algebra over floating-point numbers< I will start reading the book "matrix computation". Have you read any books by G. W. Stewart?

gundamlh gravatar imagegundamlh ( 2013-11-23 05:30:49 +0100 )edit

No i didn't, i will have a look, thanks for the reference. In general, naive implementations of the algorithms you learned in math may be sensituve to numerical stability (e.g. the unstability of the Gauss elimination depends on the conditionning of the matrix).

tmonteil gravatar imagetmonteil ( 2013-11-23 05:57:27 +0100 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2013-11-21 11:01:49 +0100

Seen: 1,417 times

Last updated: Nov 22 '13