1

对于我当前的项目,我需要能够使用来自 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

感谢您的任何想法。

4

2 回答 2

4

在 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从与其最低有效位对应的列中删除所有条目。如果该行为零,则它没有最低有效位并且对排名没有贡献,因此我们只需将其丢弃并继续。

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×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 函数,它采用 dtype 的一维 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的 NumPy 数组(64,),我得到的平均秩计算时间为 7.56µs,比纯 Python 版本快 25 倍以上。

于 2019-07-02T19:34:35.147 回答
3

我编写了一个 Python 包galois,它在 Galois 字段上扩展了 NumPy 数组。伽罗瓦域矩阵上的线性代数是预期的用例之一。它是用 Python 编写的,但使用 Numba 进行 JIT 编译以提高速度。它非常快,并且大多数线性代数例程也被编译。(一个例外,截至 2021 年 8 月 11 日,行缩减例程尚未经过 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)

NumPy 数组有一个添加的方法,称为row_reduce()对矩阵执行高斯消元。您还可以在 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 上打开一个问题。

于 2021-08-11T16:15:21.860 回答