arbitrary precision linear algebra

[geremy condra]

On Wed, Mar 2, 2011 at 10:47 AM, Ben123 <ben.is.located at gmail.com> wrote:
> On Mar 2, 12:22 pm, geremy condra <debat... at gmail.com> wrote:
>> On Wed, Mar 2, 2011 at 10:21 AM, geremy condra <debat... at gmail.com> wrote:
>> > On Wed, Mar 2, 2011 at 6:42 AM, Ben123 <ben.is.loca... at gmail.com> 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.
>>
>> >> Thanks in advance
>>
>> > 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
>>
>> Apologies, forgot the links:
>>
>> http://www.sagemath.org/doc/constructions/linear_algebra.htmlhttp://www.sagemath.org/doc/reference/sage/rings/complex_field.html
>>
>> Geremy Condra
>
> I'm not sufficiently familiar with Sage, but from
> http://www.sagemath.org/doc/constructions/linear_algebra.html
>
> "currently Sage does not implement multiprecision numerical
> eigenvalues/eigenvectors"
>
> I'll ask on the Sage forums about this. In the mean time, I'm still
> trying to get arbitrary precision linear algebra in Python

I'd suggest you read that slightly more carefully. It's speaking
specifically of RR and CC, ie, double-width reals and complex values.
Using RealField and ComplexField- the arbitrary precision real and
complex fields- seems to be working. Using the earlier example:

sage: M1 = Matrix(RealField(1000), [[0, 2], [1, 0]])
sage: M2 = Matrix(RR, [[0, 2], [1, 0]])
sage: M1.eigenvalues()
[1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157273501384623091229702492483605585073721264412149709993583141322266592750559275579995050115278206057147010955997160597027453459686201472851741864088919860955232923048430871432145083976260362799525140799,
-1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157273501384623091229702492483605585073721264412149709993583141322266592750559275579995050115278206057147010955997160597027453459686201472851741864088919860955232923048430871432145083976260362799525140799]
sage: M2.eigenvalues()
[1.41421356237310, -1.41421356237310]

Converting the first of the latter values to an element of
RealField(1000) yields much what I would expect from higher precision
arithmetic:

 R = RealField(1000)
sage: x = M1.eigenvalues()[0]
sage: y = R(M2.eigenvalues()[0])
sage: x
1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157273501384623091229702492483605585073721264412149709993583141322266592750559275579995050115278206057147010955997160597027453459686201472851741864088919860955232923048430871432145083976260362799525140799
sage: y
1.41421356237309514547462185873882845044136047363281250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

So, while I don't know for a fact that it's using the precision you
need, it certainly does seem to be using high precision arithmetic
here. Furthermore, repeating it for various precisions seems to
increase the difference, as would be expected from better
approximations, and the number of digits in the result is consistent
with the interpretation that it is using the specified precision.

All of this to say that it seems to be doing what you want it to do.

Geremy Condra