How to get the integer eigenvectors of a Numpy matrix?

Question

I have a Numpy matrix, for example, numpy.matrix([[-1, 2],[1, -2]], dtype='int'). I want to get its integer-valued eigenvectors, if any; for example, numpy.array([[-1], [1]]) for the above matrix. What Numpy returns are eigenvectors in floating numbers, scaled to have unit length.

One can do this in Sage, where one can specify the field (i.e., data type) of the matrix and operations done on the matrix will respect the field one specifies.

Any idea of how to do this nicely in Python? Many thanks in advance.

Answer 1

I am personally content with the following solution: I called sage in Python and let sage compute what I want. sage, being math-oriented, is rather versatile in computations involving fields other than reals.

Below is my script compute_intarrs.py and it requires sage be installed. Be aware it 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])]

Answer 2

You can do some pretty cool things with dtype = object and the fractions.Fraction class, e.g.

>>> 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 np.linalg module converts everything to float before doing anything, so you can't expect to get solutions directly as integers or rationals. But you can always do the following after your computations:

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 round-off 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 could mitigate some of that using the limit_denominator method of Fraction, but probably will not be all that robust.