快速计算矩阵等级超过GF(2)



对于我当前的项目,我需要能够使用GF(2(的条目计算64*64个矩阵的排名。我想知道是否有人有一个很好的解决方案。

我一直在为此使用Pyfinite,但是由于它是纯粹的Python实现,因此相当慢。我还尝试过cythonise我一直使用的代码,但由于依靠Pyfinite而遇到了问题。

我的下一个想法是在Cython写自己的课,但这似乎有点过分杀伤。

我需要以下功能

matrix = GF2Matrix(size=64) # creating a 64*64 matrix
matrix.setRow(i, [1,0,1....,1]) # set row using list
matrix += matrix2 # addition of matrices
rank(matrix) # then computing the rank

感谢任何想法。

有效地表示gf(2(上的矩阵的一种方法是将行存储为整数,将每个整数解释为小弦。因此,例如,4 x-4矩阵

[0 1 1 0]
[1 0 1 1]
[0 0 1 0]
[1 0 0 1]

(具有等级3(可以表示为整数的列表[6, 13, 4, 9]。在这里,我认为第一列是对应整数的最小一点,而最后一列则与最重要的位相对应,但是反向惯例也可以使用。

使用此表示,可以使用Python的位整数操作有效地执行行操作:添加的^&用于乘法。然后,您可以使用标准高斯淘汰方法来计算排名。

这是一些合理效率的代码。给定代表矩阵的非负整数的集合rows如上所述,我们反复删除列表中的最后一行,然后使用该行消除从列中的所有1条目,与其最小的位相对应。如果该行为零,则它没有最少的显着位,并且不会有助于等级,因此我们只是将其丢弃并继续。

def gf2_rank(rows):
    """
    Find rank of a matrix over GF2.
    The rows of the matrix are given as nonnegative integers, thought
    of as bit-strings.
    This function modifies the input list. Use gf2_rank(rows.copy())
    instead of gf2_rank(rows) to avoid modifying rows.
    """
    rank = 0
    while rows:
        pivot_row = rows.pop()
        if pivot_row:
            rank += 1
            lsb = pivot_row & -pivot_row
            for index, row in enumerate(rows):
                if row & lsb:
                    rows[index] = row ^ pivot_row
    return rank

让我们在GF2上运行一些随机64乘64矩阵的时间。random_matrices是创建随机64 x 64矩阵集合的函数:

import random
def random_matrix():
    return [random.getrandbits(64) for row in range(64)]
def random_matrices(count):
    return [random_matrix() for _ in range(count)]

这是计时代码:

import timeit
count = 1000
number = 10
timer = timeit.Timer(
    setup="ms = random_matrices({})".format(count),
    stmt="[gf2_rank(m.copy()) for m in ms]",
    globals=globals())
print(min(timer.repeat(number=number)) / count / number)

在我的机器上打印的结果(2.7 GHz Intel Core i7,MacOS 10.14.5,Python 3.7(是0.0001984686384,因此,单个排名计算的200µs触摸是一个触摸。

200µs对于纯Python等级计算是非常可观的,但是如果这不够快,我们可以遵循您的建议使用Cython。这是一个Cython函数,它采用了1D numpy数组dtype np.uint64,再次将数组的每个元素视为gf2上64 x-64矩阵的行,并返回该矩阵的等级。

# cython: language_level=3, boundscheck=False
from libc.stdint cimport uint64_t, int64_t
def gf2_rank(uint64_t[:] rows):
    """
    Find rank of a matrix over GF2.
    The matrix can have no more than 64 columns, and is represented
    as a 1d NumPy array of dtype `np.uint64`. As before, each integer
    in the array is thought of as a bit-string to give a row of the
    matrix over GF2.
    This function modifies the input array.
    """
    cdef size_t i, j, nrows, rank
    cdef uint64_t pivot_row, row, lsb
    nrows = rows.shape[0]
    rank = 0
    for i in range(nrows):
        pivot_row = rows[i]
        if pivot_row:
            rank += 1
            lsb = pivot_row & -pivot_row
            for j in range(i + 1, nrows):
                row = rows[j]
                if row & lsb:
                    rows[j] = row ^ pivot_row
    return rank

运行64 x 64矩阵的等效时间,现在以DTYPE np.uint64和Shape (64,)的numpy阵列表示,我得到的平均排名时间为7.56µs,比纯Python版本快25倍。

我写了一个python package galois,该软件包在galois字段上扩展了numpy阵列。Galois场矩阵上的线性代数是预期用例之一。它用python编写,但使用numba编制了JIT的速度。它非常快,大多数线性代数例程也已编译。(一个例外,截至08/11/2021,尚未编译JIT行减少例程,但可以添加。(

这是使用galois库来执行您描述的示例。

创建GF(2)数组类并创建一个显式数组和一个随机数组。

In [1]: import numpy as np                                                                                                                                                                     
In [2]: import galois                                                                                                                                                                          
In [3]: GF = galois.GF(2)                                                                                                                                                                      
In [4]: A = GF([[0, 0, 1, 0], [0, 1, 1, 1], [1, 0, 1, 0], [1, 0, 1, 0]]); A                                                                                                                                                                
Out[4]: 
GF([[0, 0, 1, 0],
    [0, 1, 1, 1],
    [1, 0, 1, 0],
    [1, 0, 1, 0]], order=2)
In [5]: B = GF.Random((4,4)); B                                                                                                                                                                
Out[5]: 
GF([[1, 1, 1, 0],
    [1, 1, 1, 0],
    [1, 1, 0, 0],
    [0, 0, 1, 0]], order=2)

您可以更新整个行(按照您的要求(。

In [6]: B[0,:] = [1,0,0,0]; B                                                                                                                                                                  
Out[6]: 
GF([[1, 0, 0, 0],
    [1, 1, 1, 0],
    [1, 1, 0, 0],
    [0, 0, 1, 0]], order=2)

矩阵算术与普通的二元运算符合作。这是矩阵添加和矩阵乘法。

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

在称为row_reduce()的Numpy阵列中增加了一种方法,该方法在矩阵上执行高斯消除。您也可以在Galois字段数组上调用标准Numpy线性代数函数并获得正确的结果。

In [9]: A.row_reduce()                                                                                                                                                                         
Out[9]: 
GF([[1, 0, 0, 0],
    [0, 1, 0, 1],
    [0, 0, 1, 0],
    [0, 0, 0, 0]], order=2)
In [10]: np.linalg.matrix_rank(A)                                                                                                                                                              
Out[10]: 3

希望这会有所帮助!如果需要其他功能,请在GitHub上打开一个问题。

最新更新