Быстрое вычисление ранга матрицы над GF (2) - PullRequest
Новая
       23

Быстрое вычисление ранга матрицы над GF (2)

0 голосов
/ 02 июля 2019

Для моего текущего проекта мне нужно иметь возможность рассчитать ранг 64 * 64 матриц с записями из GF (2). Мне было интересно, есть ли у кого-нибудь хорошее решение.

Я использовал для этого pyfinite , но это довольно медленно, так как это чистая реализация на python. Я также пытался цитонизировать код, который использовал, но у меня были проблемы из-за использования 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

Спасибо за любые идеи.

1 Ответ

1 голос
/ 02 июля 2019

Один из способов эффективно представить матрицу над 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

Давайте выполним некоторые тайминги для случайных матриц 64 на 64 по GF2. 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)

Результат, напечатанный на моем компьютере (Intel Core i7 с частотой 2,7 ГГц, macOS 10.14.5, Python 3.7), равен 0.0001984686384, так что для вычисления одного ранга это значение меньше 200 мкс.

200 мкс вполне респектабельно для вычисления чистых рангов Python, но в случае, если это не достаточно быстро, мы можем последовать вашему предложению использовать Cython. Вот функция Cython, которая принимает массив 1d NumPy dtype np.uint64, снова думая о каждом элементе массива как о строке вашей матрицы 64 на 64 над GF2, и возвращает ранг этой матрицы. 

# 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, теперь представленных в виде массивов NumPy типа dtype np.uint64 и формы (64,), я получаю среднее время вычисления ранга в 7,56 мкс, что в 25 раз быстрее, чем на чистом Python версия.

...