如何得到一个Numpy矩阵的整数特征向量

Question

我有一个Numpy矩阵，例如，numpy.matrix([[-1, 2],[1, -2]], dtype='int')。我想得到它的整数值特征向量，如果有的话;例如，上面矩阵的numpy.array([[-1], [1]])。Numpy返回的是浮点数的特征向量，缩放成单位长度。

可以在Sage中这样做，在Sage中可以指定矩阵的字段(即数据类型)，并且对矩阵进行的操作将尊重指定的字段。

知道如何在Python中很好地做到这一点吗?

Answer 1

我个人对以下解决方案感到满意:我在Python中调用sage，让sage计算我想要的。sage以数学为导向，在涉及实数以外的领域的计算中相当通用。

下面是我的脚本compute_intarrs.py，它需要安装sage。注意，它有点慢。

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 *nn')
f.write('A = matrix(ZZ, %s)nn' % 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

结果:

In [1]: %run compute_intarrs.py

[array([1, 1])]

Answer 2

你可以用dtype = object和fractions.Fraction类做一些很酷的事情，例如

>>> 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)

不幸的是，np.linalg模块在做任何事情之前都会将所有内容转换为float，因此您不能期望直接获得整数或有理数的解决方案。但是，在计算之后，您总是可以执行以下操作:

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)

它会很慢，并且对舍入错误非常敏感:

>>> 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)

您可以使用Fraction的limit_denominator方法来缓解其中的一些问题，但可能不会那么健壮。