Python: уменьшенная форма эшелона строк (mod p) очень большой матрицы

Я хочу найти уменьшенную форму эшелона строки (в поле F_q) большой матрицы. Я попробовал следующий код. Хотя я использовал библиотеку gmpy2 для ускорения, программе все еще не хватало памяти. потому что моя входная матрица очень большая (100 x 2 ^ 15), и p тоже очень большое (| p | = 256 бит). Может ли кто-нибудь предложить, как уменьшить сложность этого алгоритма.

Спасибо

def invmodp(a, p):
    return gmpy2.invert(a,p)

def division_mod(a, b, p): #a/b mod p
    invert = invmodp(b, p)
    return (a * invert) %p

def row_echelon_form(M, p):
   lead = 0
   rowCount = len(M)
   columnCount = len(M[0])
   for r in range(rowCount):
       if lead >= columnCount:
           return
       i = r
       while M[i][lead] == 0:
           i += 1
           if i == rowCount:
               i = r
               lead += 1
               if columnCount == lead:
                   return
    M[i],M[r] = M[r],M[i]
    lv = M[r][lead]
    M[r] = [ division_mod(mrx, lv, p) for mrx in M[r]]
    for i in range(rowCount):
        if i != r:
            lv = M[i][lead]
            M[i] = [ (iv - lv*rv)%p for rv,iv in zip(M[r],M[i])]
    lead += 1
return M

python gaussian gmpy
person santa    schedule 09.07.2015    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Мне удалось сэкономить несколько секунд времени работы, используя gmpy2.divm вместо division_mod. Я не смог сделать каких-либо других значительных улучшений. Следующая программа создает случайную матрицу 100 x 2 ^ 15 и вычисляет форму эшелона строк примерно за 3 минуты и потребляет 425 МБ памяти.

import gmpy2

bits = 256
r = 100
c = 2**15

p = gmpy2.next_prime(2**bits - 1234)
seed = gmpy2.random_state(42)

M = []
for i in range(r):
    M.append([gmpy2.mpz_urandomb(seed, bits) for j in range(c)])

def row_echelon_form(M, p):
    lead = 0
    rowCount = len(M)
    columnCount = len(M[0])
    for r in range(rowCount):
        if lead >= columnCount:
            return
        i = r
        while M[i][lead] == 0:
            i += 1
            if i == rowCount:
                i = r
                lead += 1
                if columnCount == lead:
                    return

        M[i],M[r] = M[r],M[i]
        lv = M[r][lead]
        M[r] = [ gmpy2.divm(mrx, lv, p) for mrx in M[r]]
        for i in range(rowCount):
            if i != r:
                lv = M[i][lead]
                M[i] = [ (iv - lv*rv) % p for rv,iv in zip(M[r],M[i])]
        lead += 1
    return M

N = row_echelon_form(M, p)

Если использование памяти превысит 500 МБ, в вашей версии gmpy2 может быть утечка памяти. Или я неправильно интерпретировал ваши требования и матрица значительно больше.

Отказ от ответственности: я поддерживаю gmpy2.

person casevh    schedule 10.07.2015