How do I get the integer eigenvectors of a Numpy matrix?
I have a Numpy matrix for example numpy.matrix([[-1, 2],[1, -2]], dtype='int')
. I want to get integer -valued eigenvectors, if any; for example, numpy.array([[-1], [1]])
for the above matrix. What Numpy returns is a float eigenvector, scaled to one length.
This can be done in Sage, where you can specify the field (i.e., data type) of the matrix, and the operations performed on the matrix will match the specified field.
Any idea how to do this in Python? Thank you very much in advance.
source to share
I am personally happy with the following solution: I called sage
in Python and sage
figured out what I want. sage
being mathematically oriented, it is quite versatile in calculations involving fields other than real ones.
Below is my script compute_intarrs.py
and installation is required sage
. Remember, this is a little slow.
import subprocess
import re
import numpy as np
# construct a numpy matrix
mat = np.matrix([[1,-1],[-1,1]])
# convert the matrix into a string recognizable by sage
matstr = re.sub('\s|[a-z]|\(|\)', '', mat.__repr__())
# write a (sage) python script "mat.py";
# for more info of the sage commands:
# www.sagemath.org/doc/faq/faq-usage.html#how-do-i-import-sage-into-a-python-script
# www.sagemath.org/doc/tutorial/tour_linalg.html
f = open('mat.py', 'w')
f.write('from sage.all import *\n\n')
f.write('A = matrix(ZZ, %s)\n\n' % matstr)
f.write('print A.kernel()') # this returns the left nullspace vectors
f.close()
# call sage and run mat.py
p = subprocess.Popen(['sage', '-python', 'mat.py'], stdout=subprocess.PIPE)
# process the output from sage
arrstrs = p.communicate()[0].split('\n')[2:-1]
arrs = [np.array(eval(re.sub('(?<=\d)\s*(?=\d|-)', ',', arrstr)))
for arrstr in arrstrs]
print arrs
Result:
In [1]: %run compute_intarrs.py
[array([1, 1])]
source to share
You can do some cool things with dtype = object
and fractions.Fraction
like
>>> A = np.array([fractions.Fraction(1, j) for j in xrange(1, 13)]).reshape(3, 4)
>>> A
array([[1, 1/2, 1/3, 1/4],
[1/5, 1/6, 1/7, 1/8],
[1/9, 1/10, 1/11, 1/12]], dtype=object)
>>> B = np.array([fractions.Fraction(1, j) for j in xrange(1, 13)]).reshape(4, 3)
>>> B
array([[1, 1/2, 1/3],
[1/4, 1/5, 1/6],
[1/7, 1/8, 1/9],
[1/10, 1/11, 1/12]], dtype=object)
>>> np.dot(A, B)
array([[503/420, 877/1320, 205/432],
[3229/11760, 751/4620, 1217/10080],
[1091/6930, 1871/19800, 1681/23760]], dtype=object)
Unfortunately, the module np.linalg
converts everything to float
before doing anything, so you cannot expect to get solutions directly as whole or rational. But you can always do the following after your calculations:
def scale_to_int(x) :
fracs = [fractions.Fraction(j) for j in x.ravel()]
denominators = [j.denominator for j in fracs]
lcm = reduce(lambda a, b: max(a, b) / fractions.gcd(a, b) * min(a, b),
denominators)
fracs = map(lambda x : lcm * x, fracs)
gcd = reduce(lambda a, b: fractions.gcd(a, b), fracs)
fracs = map(lambda x: x / gcd, fracs)
return np.array(fracs).reshape(x.shape)
It will be slow and very sensitive to rounding errors:
>>> scale_to_int(np.linspace(0, 1, 5)) # [0, 0.25, 0.5, 0.75, 1]
array([0, 1, 2, 3, 4], dtype=object)
>>> scale_to_int(np.linspace(0, 1, 4)) # [0, 0.33333333, 0.66666667, 1]
array([0, 6004799503160661, 12009599006321322, 18014398509481984], dtype=object)
You can mitigate some of these using the method limit_denominator
Fraction
, but probably not all will be reliable.
source to share