【问题标题】:Can I use my own Python class with numpy or some other matrix library?我可以将我自己的 Python 类与 numpy 或其他一些矩阵库一起使用吗?
【发布时间】:2011-07-31 11:51:25
【问题描述】:

我希望能够使用 Python 类作为元素进行矩阵运算——在本例中,是一个简单的 Galois field 实现。它实现了必要的__add____mul____sub__等。

起初,我认为numpy arrays 应该可以做到这一点,使用dtype 参数,但从the dtype documentation 看来,dtype 似乎不能是任意Python 类。例如,我有一个类 Galois 进行模 2 运算:

>>> from galois import Galois
>>> Galois(1) + Galois(0)
Galois(1)
>>> Galois(1) + Galois(1)
Galois(0)

我可以尝试在 numpy 中使用它:

>>> import numpy as np
>>> a = np.identity(4, Galois)
>>> a
array([[1, 0, 0, 0],
       [0, 1, 0, 0],
       [0, 0, 1, 0],
       [0, 0, 0, 1]], dtype=object)

但是如果我对矩阵进行操作,元素不会遵循我的类的方法:

>>> b = np.identity(4, Galois)
>>> a+b
array([[2, 0, 0, 0],
       [0, 2, 0, 0],
       [0, 0, 2, 0],
       [0, 0, 0, 2]], dtype=object)

有什么办法可以用 numpy 来完成这项工作?

是否有任何其他 Python 矩阵库可以对任意类数类进行矩阵运算(包括求逆)?

更新

感谢到目前为止的回答。但我仍然无法像我希望的那样真正使用它。加法和乘法看起来不错,但不是矩阵求逆。例如,让我们尝试从forward S-box affine transform matrix 中获取AES inverse S-box affine transform matrix

class Galois(object):
    MODULO = 2

    def __init__(self, val):
        self.val = int(val) % self.MODULO

    def __add__(self, val):
        return self.__class__((self.val + int(val)) % self.MODULO)
    def __sub__(self, val):
        return self.__class__((self.val - int(val)) % self.MODULO)
    def __mul__(self, val):
        return self.__class__((self.val * int(val)) % self.MODULO)
    def __int__(self):
        return self.val
    def __repr__(self):
        return "%s(%d)" % (self.__class__.__name__, self.val)
    def __float__(self):
        return float(self.val)

if __name__ == "__main__":
    import numpy as np

    Gv = np.vectorize(Galois)

    a = Gv(np.identity(8)) + Gv(np.eye(8,8,-1)) + Gv(np.eye(8,8,-2)) + Gv(np.eye(8,8,-3)) + Gv(np.eye(8,8,-4)) + Gv(np.eye(8,8,4)) + Gv(np.eye(8,8,5)) + Gv(np.eye(8,8,6)) + Gv(np.eye(8,8,7))
    print np.matrix(a)
    print np.matrix(a).I

结果:

[[Galois(1) Galois(0) Galois(0) Galois(0) Galois(1) Galois(1) Galois(1)
  Galois(1)]
 [Galois(1) Galois(1) Galois(0) Galois(0) Galois(0) Galois(1) Galois(1)
  Galois(1)]
 [Galois(1) Galois(1) Galois(1) Galois(0) Galois(0) Galois(0) Galois(1)
  Galois(1)]
 [Galois(1) Galois(1) Galois(1) Galois(1) Galois(0) Galois(0) Galois(0)
  Galois(1)]
 [Galois(1) Galois(1) Galois(1) Galois(1) Galois(1) Galois(0) Galois(0)
  Galois(0)]
 [Galois(0) Galois(1) Galois(1) Galois(1) Galois(1) Galois(1) Galois(0)
  Galois(0)]
 [Galois(0) Galois(0) Galois(1) Galois(1) Galois(1) Galois(1) Galois(1)
  Galois(0)]
 [Galois(0) Galois(0) Galois(0) Galois(1) Galois(1) Galois(1) Galois(1)
  Galois(1)]]
[[ 0.4  0.4 -0.6  0.4  0.4 -0.6  0.4 -0.6]
 [-0.6  0.4  0.4 -0.6  0.4  0.4 -0.6  0.4]
 [ 0.4 -0.6  0.4  0.4 -0.6  0.4  0.4 -0.6]
 [-0.6  0.4 -0.6  0.4  0.4 -0.6  0.4  0.4]
 [ 0.4 -0.6  0.4 -0.6  0.4  0.4 -0.6  0.4]
 [ 0.4  0.4 -0.6  0.4 -0.6  0.4  0.4 -0.6]
 [-0.6  0.4  0.4 -0.6  0.4 -0.6  0.4  0.4]
 [ 0.4 -0.6  0.4  0.4 -0.6  0.4 -0.6  0.4]]

不是我希望的结果。似乎对于矩阵求逆,numpy 只是将矩阵转换为浮点数,然后用普通实数进行求逆。

【问题讨论】:

    标签: python matrix numpy


    【解决方案1】:

    您可以将object 用作dtype,这将允许任意Python 对象。我不认为有任何方法可以专门化一个 numpy 数组来只接受一个特定类的 Python 对象。

    【讨论】:

    • 谢谢。我用一个例子更新了我的问题。我错过了什么吗?
    • 问题在于你的单位矩阵的创建,而不是它的操作。如果您用Galois 对象强制替换ab 数组的每个元素,那么您会发现a+b 为您提供了一个充满Galois(0) 的数组。
    • 谢谢!这是有道理的,而且似乎有效。虽然我希望 numpy 的 np.identitynp.eye 等会使用其 __init__ 函数创建 dtype 类型的元素。
    • 明确使用dtype=object。如果你给我们一个 dtype 的类,我们只是假设你的意思是dtype=object。抱歉,我们没有为您的班级制作特殊的 dtype。
    【解决方案2】:

    您是否查看过sage,特别是galois_group

    看来您是在重新发明轮子。但如果你坚持这样做,你可以考虑subclass ndarray

    【讨论】:

    • 圣人无疑是优秀的。但是我认为 Python 比 sage 有一定的优势——例如更轻量级,更适合 Windows。但无论如何......你能确认圣人galois_group 在圣人矩阵中表现良好吗?
    • @Craig:我无法确认任何事情。但是现在看到您的更新,我希望sage 将是您可以拥有的最多的pythonic 工具。试一试,自己判断。谢谢
    【解决方案3】:

    我不知道有一种方法可以让矩阵元素根据任意 Python 类运行。但是,可以全局更改某些操作的行为,如以下示例所示,

    import numpy as np
    from numpy import set_numeric_ops
    
    from numpy import poly1d as poly
    from numpy import identity as idt
    
    def gfadd(x,y):
        return np.add(x,y) % 2
    
    set_numeric_ops(add=gfadd)
    
    a = idt(4,np.int)
    print a+a
    

    生产,

    [[0 0 0 0]
     [0 0 0 0]
     [0 0 0 0]
     [0 0 0 0]]
    

    和,

    p = poly([1,0,1])
    print p+p
    

    给予,

    0
    

    但是,您可能希望继承 ndarray

    【讨论】:

      【解决方案4】:

      以下是如何使用另一个数组创建和初始化 Numpy 对象数组:

      import numpy as np
      
      class G:
          def __init__(self, x):
              self.x = x
      
      I = np.identity(5)
      Gv = np.vectorize(G)
      GG = Gv(I)
      
      print GG[0,0].x
      print GG[0,1].x
      

      【讨论】:

      • 然而,object dtypes 很慢并且通常被避免。创建具有 NumPy 数组属性的对象是更好的做法。
      【解决方案5】:

      关于矩阵求逆的更新:使用 NumPy 矩阵求逆对伽罗瓦域上的矩阵求逆是行不通的。 NumPy 将实际反转矩阵的任务委托给LAPACK,这是一个用 Fortran 编写的线性代数库。 LAPACK 当然完全不知道 Python 类和运算符,并且永远无法调用 Galois 类的方法。此外,他们的矩阵求逆算法使用了比较运算符(如<>),这对伽罗瓦域的元素没有意义。

      因此,您可以选择自己实现矩阵求逆或使用其中一种可用的实现。例如,SymPy 对 Galois 域的支持有限。 PARI/GP 支持伽罗瓦域和some Python bindings

      【讨论】:

      • SymPy 也有对 Galois 域的undocumented 支持。但到目前为止,我发现有sympy.polys.domains.FF,所以我可以做GF2 = sympy.polys.domains.FF(2),然后做GF2(1) + GF2(1),这给了我SymmetricModularInteger2(0)。但是,我不知道如何在 SymPy 矩阵中使用它。
      • 似乎 SymPy 还没有被设计用于处理其他类型的矩阵(除非我遗漏了什么),但 thisthis 表明其他人也想要,并且正在考虑关于它。
      • @Craig:感谢您的反馈。似乎没有太多人对有限域上的矩阵算术感兴趣。如果您只需要对密集矩阵进行矩阵求逆,那么我建议您自己实现高斯消除。这不是太多的代码,而且在有限的字段上你甚至不需要旋转。
      猜你喜欢
      • 1970-01-01
      • 2011-01-29
      • 2019-06-18
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2014-09-11
      相关资源
      最近更新 更多