首页 > 解决方案 > 带有对象的 Numpy 矩阵求逆

问题描述

我正在使用 gf256 库进行 galois 场数学运算,并将它放在一个 numpy 矩阵中。尽管使用它调用np.linalg.inv()时,它会引发错误。

以上是总结,下面是详细内容:

import numpy as np
from gf256 import GF256 as gf
npgf = np.vectorize(gf)

arr = np.identity(4, np.uint8) * 10
gfarr = npgf(arr)

毕竟,gfarr看起来像这样

array([[GF256(0b00001010), GF256(0b00000000), GF256(0b00000000),
        GF256(0b00000000)],
       [GF256(0b00000000), GF256(0b00001010), GF256(0b00000000),
        GF256(0b00000000)],
       [GF256(0b00000000), GF256(0b00000000), GF256(0b00001010),
        GF256(0b00000000)],
       [GF256(0b00000000), GF256(0b00000000), GF256(0b00000000),
        GF256(0b00001010)]], dtype=object)

np.linalg.inv(gfarr)抛出此错误

Traceback (most recent call last):
  File "<pyshell#152>", line 1, in <module>
    np.linalg.inv(gfarr)
  File "[python3.6]\lib\site-packages\numpy\linalg\linalg.py", line 528, in inv
    ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
TypeError: No loop matching the specified signature and casting
was found for ufunc inv

矩阵绝对是可逆的,GF256 类支持所有常用运算符。是否可以使用 numpy 进行这项工作?

标签: pythonnumpygalois-fieldfinite-field

解决方案


np.linalg.inv将使用浮点数调用矩阵求逆的 BLAS/LAPACK 实现,但您需要在矩阵求逆过程中使用伽罗瓦域算术。为此,NumPy 数组需要拦截或覆盖对np.linalg.invin的调用__array_function__()。的矩阵求逆A可以使用[A | I]伽罗瓦域上的高斯消元法完成,得到[I | A^-1].

我有一个类似的用例,所以我创建了一个名为galois的 Python 包,它在 Galois 字段上扩展了 NumPy 数组。它使用 Numba 将 NumPy ufunc 替换为 JIT 编译的 ufunc。这意味着数组算术与普通 NumPy 算术一样快,或者几乎一样快。请参阅此性能比较

它还支持线性代数并覆盖相关np.linalg函数。因此,您正在寻找的矩阵求逆开箱即用。这是使用您的矩阵的示例。

In [1]: import numpy as np                                                                                                                                                                     

In [2]: import galois                                                                                                                                                                          

In [3]: GF = galois.GF(2**8)                                                                                                                                                                   

In [4]: print(GF.properties)                                                                                                                                                                   
GF(2^8):
  characteristic: 2
  degree: 8
  order: 256
  irreducible_poly: x^8 + x^4 + x^3 + x^2 + 1
  is_primitive_poly: True
  primitive_element: x

In [5]: A = GF.Identity(4) * GF(10); A                                                                                                                                                         
Out[5]: 
GF([[10,  0,  0,  0],
    [ 0, 10,  0,  0],
    [ 0,  0, 10,  0],
    [ 0,  0,  0, 10]], order=2^8)

In [6]: A_inv = np.linalg.inv(A); A_inv                                                                                                                                                        
Out[6]: 
GF([[221,   0,   0,   0],
    [  0, 221,   0,   0],
    [  0,   0, 221,   0],
    [  0,   0,   0, 221]], order=2^8)

In [7]: A @ A_inv                                                                                                                                                                              
Out[7]: 
GF([[1, 0, 0, 0],
    [0, 1, 0, 0],
    [0, 0, 1, 0],
    [0, 0, 0, 1]], order=2^8)

推荐阅读