GF（2）上矩阵秩的快速计算

Question

对于我当前的项目，我需要能够计算64*64矩阵的秩与GF（2）的条目。我想知道是否有人能找到一个好的解决办法。在

我一直在使用pyfinite来实现这一点，但由于它是纯python实现，所以速度相当慢。我也尝试过对我一直使用的代码进行cythonis，但是由于依赖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

谢谢你的建议。在

Answer 1

有效表示GF（2）上矩阵的一种方法是将行存储为整数，将每个整数解释为位字符串。例如，4乘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项。如果该行为零，那么它没有最低有效位，并且对排名没有贡献，所以我们简单地丢弃它并继续。在

^{pr2}$

让我们对GF2上的随机64×64矩阵运行一些计时。random_matrices是一个创建随机64×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.7ghz英特尔酷睿i7，macOS 10.14.5，python3.7）上打印的结果是0.0001984686384，因此对于单列计算来说，这是一个200微秒以下的触摸。在

在这种情况下，我们可以很快地使用Python 200来计算。这是一个Cython函数，它接受一个1d NumPy数组np.uint64，同样将数组中的每个元素看作是GF2上64×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×64矩阵（现在表示为dtypenp.uint64和shape(64,)的NumPy数组）运行等效计时，我得到的平均秩计算时间为7.56µs，比纯Python版本快25倍以上。在