Ask Your Question
0

Set the precision of imported methods

asked 11 years ago

gundamlh gravatar image

updated 11 years ago

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!

Preview: (hide)

1 Answer

Sort by » oldest newest most voted
2

answered 11 years ago

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.

Preview: (hide)
link

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 ( 11 years ago )

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 ( 11 years ago )

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: 11 years ago

Seen: 1,466 times

Last updated: Nov 22 '13