need eigenvectors of complex matrix to arbitrary precision

[Ben123]

Hello. I've written a sage program which produces a complex matrix. I
want to find the eigenvalues and associated eigenvectors. I also want
to use arbitrary precision. I don't care about speed. I've read old
posts to this group on this topic, but am unsure how to proceed.
Currently I'm using the following method and using sage 4.6.1

precision_digits=30
nop=5 # rank of matrix
MS_nop_comp=MatrixSpace(ComplexField(precision_digits),nop,nop)
tmat=MS_nop_comp(0) # zero-ize the values
ttdag=MS_nop_comp(0)

# I realize there are more efficient methods of getting a random
matrix, but this is explicit
for a in range(nop):
for b in range(nop):
tmat[a,b]=random()+I*random()

ttdag=tmat*tmat.conjugate().transpose() # get a Hermitian matrix
print 'ttdag is'
print ttdag
print 'eigenvalues of ttdag are '
print ttdag.eigenvalues() # eigenvalues of Hermitian matrix should be
real. Imaginary component is due to finite precision.
# I can get better precision here by increasing precision_digits

#print ttdag.eigenmatrix_right()
# IndexError: list index out of range

print ttdag.eigenvectors_right()
# this is not returning the eigenvectors, even when precision is
increased to 500

How can I find the eigenvectors of a complex Hermitian matrix with
arbitrary precision?

Thanks,

[Ben123]

Update: after reading #10346 on
http://www.sagemath.org/mirror/src/changelogs/sage-4.6.2.txt
I upgraded to 4.6.2 and am still having the same problem (no
eigenvectors specified, even with 500 digits of precision).

sage: version()
'Sage Version 4.6.2, Release Date: 2011-02-25'
sage: !uname -a
Linux ben-desktop 2.6.31-14-generic #48-Ubuntu SMP Fri Oct 16 14:05:01
UTC 2009 x86_64 GNU/Linux

[Jason Grout]

One option: you might look at using the alglib library; at one time, the
author was writing a Sage interface, but that work has stalled for
several months. Alglib appears to have a python interface, though.

Alglib: http://www.alglib.net/

Rob Beezer's been doing some work on the numerical linear algebra in
Sage, so he also might have something to add...

Thanks,

Jason

[achrzesz]

sage: gp_console()
? \p 100
realprecision = 115 significant digits (100 digits displayed)
? a=matrix(3,3,k,m,random(1.0))+I*matrix(3,3,k,m,random(1.0));
? m=a*conj(a)~;
? mateigen(m)

[Ben123]

Hello. Thank you for showing me the equivalent process in PARI/GP. I
think this implies I need to transfer the complex Hermitian matrix
into gp_console() to find eigenvectors, then transfer them back to
Sage. I'll update this post when I figure out how to do that.

In the mean time, if someone has a simple way of doing that please let
me know. My updated Sage code lacks the last step.

precision_digits=30
nop=5 # rank of matrix
MS_nop_comp=MatrixSpace(ComplexField(precision_digits),nop,nop)
tmat=MS_nop_comp(0) # zero-ize the values
ttdag=MS_nop_comp(0)

for a in range(nop):
for b in range(nop):
tmat[a,b]=random()+I*random()
ttdag=tmat*tmat.conjugate().transpose() # get a Hermitian matrix

# Find eigenvectors of ttdag in PARI/GP, pass results back to sage

Thanks for the help,

[achrzesz]

# NO WARRANTY

precision_digits=30
nop=5 # rank of matrix
MS_nop_comp=MatrixSpace(ComplexField(precision_digits),nop,nop)
tmat=MS_nop_comp(0) # zero-ize the values
ttdag=MS_nop_comp(0)
for a in range(nop):
for b in range(nop):
tmat[a,b]=random()+I*random()
ttdag=tmat*tmat.conjugate().transpose()

m=gp(ttdag)
ev=gp.mateigen(m).sage()
print ev[:-1]

[achrzesz]

Since GP/PARI default precision is 38 (on my system)
(and I dont know how to set \p 100 from sage)
it would be probably better to replace
precision_digits=30
by
precision_digits=38

[achrzesz]

#NOT CHECKED
precision_digits=100

nop=5 # rank of matrix
MS_nop_comp=MatrixSpace(ComplexField(precision_digits),nop,nop)
tmat=MS_nop_comp(0) # zero-ize the values
ttdag=MS_nop_comp(0)
for a in range(nop):
for b in range(nop):
tmat[a,b]=random()+I*random()
ttdag=tmat*tmat.conjugate().transpose()
m=gp(ttdag)

gp('default(realprecision,100)')
ev=gp.mateigen(m).sage()
print ev[:-1]

[achrzesz]

Compare such version

precision_digits=100
nop=5 # rank of matrix
MS_nop_comp=MatrixSpace(ComplexField(precision_digits),nop,nop)

tmat=MS_nop_comp.random_element()

ttdag=tmat*tmat.conjugate().transpose()
m=gp(ttdag)
gp('default(realprecision,100)')
ev=gp.mateigen(m).sage()
print ev[:-1]

[Ben123]

Thanks for the implementation. It would have taken me a long time to
figure that out. It appears the "ev" (eigenvectors) are returned, with
the last one being normalized to have entries of unity.