在 GF(2) 上快速计算矩阵秩
Fast computation of matrix rank over GF(2)
对于我当前的项目，我需要能够使用 GF(2) 中的条目计算 64*64 矩阵的秩。不知道有没有人有好的解决办法
我一直在为此使用 pyfinite，但它相当慢，因为它是一个纯粹的 python 实现。我也尝试过对我一直在使用的代码进行 cythonise，但由于依赖 pyfinite 而遇到了问题。
我的下一个想法是用 cython 编写我自己的 class，但这对我需要的东西来说似乎有点过分了。
我需要以下功能
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) 上有效表示矩阵的一种方法是将行存储为整数，将每个整数解释为 bit-string。例如，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 np.uint64 的 1d NumPy 数组，再次将数组的每个元素视为 GF2 上的 64×64 矩阵的一行，以及 returns 的等级矩阵.
# 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 矩阵的等效时间，现在表示为 dtype np.uint64 和形状 (64,) 的 NumPy 数组，我得到平均 rank-computation 时间7.56µs，比纯 Python 版本快 25 倍以上。
我编写了一个 Python 程序包 galois，它在 Galois 域上扩展了 NumPy 数组。 Galois 域矩阵上的线性代数是预期用例之一。它是用 Python 编写的，但使用 Numba 进行 JIT 编译以提高速度。它非常快，而且大多数线性代数例程也已编译。 （一个例外，截至 08/11/2021，行缩减例程尚未经过 JIT 编译，但可以添加。）
这是一个使用 galois 库来完成您所描述内容的示例。
创建一个GF(2)数组class并创建一个显式数组和一个随机数组。
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.
上提出问题
对于我当前的项目，我需要能够使用 GF(2) 中的条目计算 64*64 矩阵的秩。不知道有没有人有好的解决办法
我一直在为此使用 pyfinite，但它相当慢，因为它是一个纯粹的 python 实现。我也尝试过对我一直在使用的代码进行 cythonise，但由于依赖 pyfinite 而遇到了问题。
我的下一个想法是用 cython 编写我自己的 class，但这对我需要的东西来说似乎有点过分了。
我需要以下功能
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) 上有效表示矩阵的一种方法是将行存储为整数，将每个整数解释为 bit-string。例如，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 np.uint64 的 1d NumPy 数组，再次将数组的每个元素视为 GF2 上的 64×64 矩阵的一行，以及 returns 的等级矩阵.
# 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 矩阵的等效时间，现在表示为 dtype np.uint64 和形状 (64,) 的 NumPy 数组，我得到平均 rank-computation 时间7.56µs，比纯 Python 版本快 25 倍以上。
我编写了一个 Python 程序包 galois，它在 Galois 域上扩展了 NumPy 数组。 Galois 域矩阵上的线性代数是预期用例之一。它是用 Python 编写的，但使用 Numba 进行 JIT 编译以提高速度。它非常快，而且大多数线性代数例程也已编译。 （一个例外，截至 08/11/2021，行缩减例程尚未经过 JIT 编译，但可以添加。）
这是一个使用 galois 库来完成您所描述内容的示例。
创建一个GF(2)数组class并创建一个显式数组和一个随机数组。
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.
上提出问题