Velocity Reviews > arbitrary precision linear algebra

# arbitrary precision linear algebra

Ben123
Guest
Posts: n/a

 03-02-2011
Hello. I have a written Python program which currently uses numpy to
perform linear algebra operations. Specifically, I do matrix*matrix,
matrix*vector, numpy.linalg.inv(matrix), and linalg.eig(matrix)
operations. Now I am interested in allowing arbitrary precision. I
have tried gmpy, bigfloat, mpmath, and decimal but I have been unable
to easily implement any with my current program. I suspect I have to
change some commands but I am unsure what.

My question is which of the arbitrary precision implementations will
most easily handle linear algebra? I don't care about speed, just ease
of use. Online tutorials for arbitrary precision linear algebra
operations would be useful.

For example, it looks like mpmath can handle matrix operations
http://fredrik-j.blogspot.com/search?q=matrix
but I was unable to find a clear tutorial. The tutorials for most of
the arbitrary precision implementations demonstrate simple scalar
examples.

Ben123
Guest
Posts: n/a

 03-02-2011
On Mar 2, 8:42*am, Ben123 <(E-Mail Removed)> wrote:
> Hello. I have a written Python program which currently uses numpy to
> perform linear algebra operations. Specifically, I do matrix*matrix,
> matrix*vector, numpy.linalg.inv(matrix), and linalg.eig(matrix)
> operations. Now I am interested in allowing arbitrary precision. I
> have tried gmpy, bigfloat, mpmath, and decimal but I have been unable
> to easily implement any with my current program. I suspect I have to
> change some commands but I am unsure what.
>
> My question is which of the arbitrary precision implementations will
> most easily handle linear algebra? I don't care about speed, just ease
> of use. Online tutorials for arbitrary precision linear algebra
> operations would be useful.
>
> For example, it looks like mpmath can handle matrix operationshttp://fredrik-j.blogspot.com/search?q=matrix
> but I was unable to find a clear tutorial. The tutorials for most of
> the arbitrary precision implementations demonstrate simple scalar
> examples.
>

Hello again. I forgot to mention I'm using
Python 2.6.4
mpmath 0.17
bigfloat 0.2.1
gmp 5.01
gmpy2 2.0.0a1
mpfr 3.0.0
all on Ubuntu x64

Arthur Mc Coy
Guest
Posts: n/a

 03-02-2011
What do you mean by "arbitrary precision" ? Each method of calculating
of something has its own precision...

Ben123
Guest
Posts: n/a

 03-02-2011
On Mar 2, 9:04*am, Arthur Mc Coy <(E-Mail Removed)> wrote:
> What do you mean by "arbitrary precision" ? Each method of calculating
> of something has its own precision...

If you are unfamiliar with arbitrary precision, I'm referring to
http://en.wikipedia.org/wiki/Arbitra...ion_arithmetic

Suppose I find the eigenvalues of a matrix and the eigenvalues range
from 1 to 0.0001. This can be handled by numpy in Python because the
smallest eigenvalue is larger than then numerical precision of 1E-19.
However, if the range of eigenvalues is 1 to 1E-40, then I will need
to increase the precision of all calculations leading up to finding
the eigenvalues.

I am working with complex valued matrices and I expect to get real
eigenvalues back (based on the physics of the system). The precision
of numpy is apparent from the imaginary component of the eigenvalues I
find, currently 1E-19 or 1E-20. I need better precision for small
eigenvalues.

In case you are curious, the complex-valued matrices are 20x20.

Thanks

Arthur Mc Coy
Guest
Posts: n/a

 03-02-2011
On Mar 2, 5:26*pm, Ben123 <(E-Mail Removed)> wrote:
> On Mar 2, 9:04*am, Arthur Mc Coy <(E-Mail Removed)> wrote:
>
> > What do you mean by "arbitrary precision" ? Each method of calculating
> > of something has its own precision...

>
> If you are unfamiliar with arbitrary precision, I'm referring tohttp://en..wikipedia.org/wiki/Arbitrary-precision_arithmetic
>
> Suppose I find the eigenvalues of a matrix and the eigenvalues range
> from 1 to 0.0001. This can be handled by numpy in Python because the
> smallest eigenvalue is larger than then numerical precision of 1E-19.
> However, if the range of eigenvalues is 1 to 1E-40, then I will need
> to increase the precision of all calculations leading up to finding
> the eigenvalues.
>
> I am working with complex valued matrices and I expect to get real
> eigenvalues back (based on the physics of the system). The precision
> of numpy is apparent from the imaginary component of the eigenvalues I
> find, currently 1E-19 or 1E-20. I need better precision for small
> eigenvalues.
>
> In case you are curious, the complex-valued matrices are 20x20.
>
> Thanks

You probably have to change the method of finding eigenvalues.
Which one do you use? Power or algebraic ?
Do you use Gaussian method to simplify matrices ?

Languages can't support infinitely large or small numbers, so try to
multiply the inner variables by 10^n to increase their values if this
will not involve on the method. For example, I did this when was
calculating geometric means of computer benchmarks.
In such way you will be storing the number of zeros as n.

Yes, interesting what are you calculating.

Ben123
Guest
Posts: n/a

 03-02-2011
On Mar 2, 10:21*am, Arthur Mc Coy <(E-Mail Removed)> wrote:
> On Mar 2, 5:26*pm, Ben123 <(E-Mail Removed)> wrote:
>
>
>
>
>
>
>
>
>
> > On Mar 2, 9:04*am, Arthur Mc Coy <(E-Mail Removed)> wrote:

>
> > > What do you mean by "arbitrary precision" ? Each method of calculating
> > > of something has its own precision...

>
> > If you are unfamiliar with arbitrary precision, I'm referring tohttp://en.wikipedia.org/wiki/Arbitrary-precision_arithmetic

>
> > Suppose I find the eigenvalues of a matrix and the eigenvalues range
> > from 1 to 0.0001. This can be handled by numpy in Python because the
> > smallest eigenvalue is larger than then numerical precision of 1E-19.
> > However, if the range of eigenvalues is 1 to 1E-40, then I will need
> > to increase the precision of all calculations leading up to finding
> > the eigenvalues.

>
> > I am working with complex valued matrices and I expect to get real
> > eigenvalues back (based on the physics of the system). The precision
> > of numpy is apparent from the imaginary component of the eigenvalues I
> > find, currently 1E-19 or 1E-20. I need better precision for small
> > eigenvalues.

>
> > In case you are curious, the complex-valued matrices are 20x20.

>
> > Thanks

>
> You probably have to change the method of finding eigenvalues.
> Which one do you use? Power or algebraic ?

I'm not sure what you mean by this. As I mentioned, in Python I am
using linalg.eig() from numpy on complex matrices. I have not
investigated how this is implemented.

> Do you use Gaussian method to simplify matrices ?

No

>
> Languages can't support infinitely large or small numbers, so try to
> multiply the inner variables by 10^n to increase their values if this
> will not involve on the method. For example, I did this when was
> calculating geometric means of computer benchmarks.

Currently I have values between 1 and 1E-300 (not infinitely small). I
don't see how scaling by powers of 10 will increase precision.

> In such way you will be storing the number of zeros as n.

Are you saying python cares whether I express a number as 0.001 or
scaled by 10^5 to read 100? If this is the case, I'm still stuck. I
need the full range of eigenvalues from 1 to 1E-300, so the entire
range could be scaled by 1E300 but I would still need better precision
than 1E19

>
> Yes, interesting what are you calculating.

Arthur Mc Coy
Guest
Posts: n/a

 03-02-2011
> Are you saying python cares whether I express a number as 0.001 or
> scaled by 10^5 to read 100? If this is the case, I'm still stuck. I
> need the full range of eigenvalues from 1 to 1E-300, so the entire
> range could be scaled by 1E300 but I would still need better precision
> than 1E19

If python package can't compute 1E-300, then you can't use it.
handled by specific python package.

Then you need to think of the way to do that. Mathematically. If no
Programming is limiting, that is true.

I have heard about grid computing. You may find scientists working on
But I think they do not an answer as well.

Also google uses its matrix to rank web pages. It computes the maximum
eigenvalue from the matrix which contain near zero entries too. Maybe
you can find how do they store those values.

Sorry, can't help anymore. I also have computing problems which I
can't yet solve

Robin Becker
Guest
Posts: n/a

 03-02-2011
On 02/03/2011 16:39, Ben123 wrote:
............
>> Languages can't support infinitely large or small numbers, so try to
>> multiply the inner variables by 10^n to increase their values if this
>> will not involve on the method. For example, I did this when was
>> calculating geometric means of computer benchmarks.

>
> Currently I have values between 1 and 1E-300 (not infinitely small). I
> don't see how scaling by powers of 10 will increase precision.
>
>> In such way you will be storing the number of zeros as n.

>
> Are you saying python cares whether I express a number as 0.001 or
> scaled by 10^5 to read 100? If this is the case, I'm still stuck. I
> need the full range of eigenvalues from 1 to 1E-300, so the entire
> range could be scaled by 1E300 but I would still need better precision
> than 1E19
>

........

If you enter a number as 1e-19 then python will treat as a float; by default I
think that float is IEEE double precision so you're getting a 48 bit mantissa
(others may have better details). That means you've already lost any idea of
arbitrary precision.

When you say you have numbers like 1E-300 are those actually numerically zero or
have you some valid inputs that vary over a huge range. It should be possible to
compute determinants/inverses etc to arbitrary precision as those are known to
be polynomial functions of the input elements. However, eigenvalues are not.

eg

[0 2]
[1 0]

has eigenvalues +/- sqrt(2) so even though you can represent the matrix in
finite precision the eigenvalues require infinite precision.

Eigenvalues are roots of a polynomial in the elements and root solving may
require an infinite number of steps so it will be difficult with arbitrary
matrices to keep arbitrary precision.
--
Robin Becker

Ben123
Guest
Posts: n/a

 03-02-2011
On Mar 2, 11:28*am, Robin Becker <(E-Mail Removed)> wrote:
> On 02/03/2011 16:39, Ben123 wrote:
> ...........
>
>
>
>
>
>
>
> >> Languages can't support infinitely large or small numbers, so try to
> >> multiply the inner variables by 10^n to increase their values if this
> >> will not involve on the method. For example, I did this when was
> >> calculating geometric means of computer benchmarks.

>
> > Currently I have values between 1 and 1E-300 (not infinitely small). I
> > don't see how scaling by powers of 10 will increase precision.

>
> >> In such way you will be storing the number of zeros as n.

>
> > Are you saying python cares whether I express a number as 0.001 or
> > scaled by 10^5 to read 100? If this is the case, I'm still stuck. I
> > need the full range of eigenvalues from 1 to 1E-300, so the entire
> > range could be scaled by 1E300 but I would still need better precision
> > than 1E19

>
> .......
>
> If you enter a number as 1e-19 then python will treat as a float; by default I
> think that float is IEEE double precision so you're getting a 48 bit mantissa
> (others may have better details). That means you've already lost any ideaof
> arbitrary precision.
>
> When you say you have numbers like 1E-300 are those actually numerically zero or
> have you some valid inputs that vary over a huge range. It should be possible to
> compute determinants/inverses etc to arbitrary precision as those are known to
> be polynomial functions of the input elements. However, eigenvalues are not.

I meant that I expect the eigenvalues to be in a range from 1 to
1E-300. The values in the matrix are from sine and cosine values and
have range [-10,10] for the real and imaginary parts. Thus, I could
specify the matrix elements to arbitrary precision prior to
diagonalization. However, I'm looking for the resulting eigenvalues to
have precision to 1E-300

>
> eg
>
> [0 2]
> [1 0]
>
> has eigenvalues +/- sqrt(2) so even though you can represent the matrix in
> finite precision the eigenvalues require infinite precision.

Here's an example of a matrix I am diagonalizing:
from numpy import *
exmpl=array([[ 9.39582700e-04 +1.21760986e-21j, 3.32958159e-04
-6.03359588e-04j, \
-3.71551527e-04 -1.77143102e-04j, 2.87113322e-04
-9.84374562e-04j], \
[ 3.32958159e-04 +6.03359588e-04j, 5.05441331e-04
+6.80604208e-21j, \
-1.79123354e-05 -3.01368276e-04j, 7.33866807e-04
-1.64459142e-04j], \
[ -3.71551527e-04 +1.77143102e-04j, -1.79123354e-05
+3.01368276e-04j, \
1.80324968e-04 -2.63622461e-21j, 7.20508934e-05
+4.43394732e-04j], \
[ 2.87113322e-04 +9.84374562e-04j, 7.33866807e-04
+1.64459142e-04j, \
7.20508934e-05 -4.43394732e-04j, 1.11903651e-03
-6.61744490e-21j]])

The eigenvalues I get back using linalg.eig(exmpl)[0] are

2.74438550e-03 +7.67536143e-20j
6.54324852e-12 +2.03438929e-20j
1.39471366e-13 +4.68335525e-20j
1.76591707e-12 +2.20243218e-20j])

Here I would want to have better precision in the eigenvalues, a
smaller imaginary component.

To use your example, I'd like to increase the number of digits shown
in the eigenvalue:
from numpy import *
test=array([[0, 2],[1, 0]])
linalg.eig(test)[0]

but using

import mpmath
from mpmath import *
mp.dps = 50
from numpy import *
test=array([[0, 2],[1, 0]])
linalg.eig(test)[0]

does not help

>
> Eigenvalues are roots of a polynomial in the elements and root solving may
> require an infinite number of steps so it will be difficult with arbitrary
> matrices to keep arbitrary precision.
> --
> Robin Becker

geremy condra
Guest
Posts: n/a

 03-02-2011
On Wed, Mar 2, 2011 at 6:42 AM, Ben123 <(E-Mail Removed)> wrote:
> Hello. I have a written Python program which currently uses numpy to
> perform linear algebra operations. Specifically, I do matrix*matrix,
> matrix*vector, numpy.linalg.inv(matrix), and linalg.eig(matrix)
> operations. Now I am interested in allowing arbitrary precision. I
> have tried gmpy, bigfloat, mpmath, and decimal but I have been unable
> to easily implement any with my current program. I suspect I have to
> change some commands but I am unsure what.
>
> My question is which of the arbitrary precision implementations will
> most easily handle linear algebra? I don't care about speed, just ease
> of use. Online tutorials for arbitrary precision linear algebra
> operations would be useful.
>
> For example, it looks like mpmath can handle matrix operations
> http://fredrik-j.blogspot.com/search?q=matrix
> but I was unable to find a clear tutorial. The tutorials for most of
> the arbitrary precision implementations demonstrate simple scalar
> examples.
>

Have you looked at Sage[0]? I don't know for a fact, but you should be
able to define a matrix over RealField(precision_in_bits) and then
take the eigenvalue of it. I don't know if it will actually produce
the precision you need though.

Geremy Condra