Полиномиальное умножение БПФ в Python с использованием встроенного Numpy.fft

Question

Я хочу быстро умножить два полинома в python. Поскольку мои многочлены довольно большие (> 100000) элементов, и я должен умножить их много. Ниже вы найдете мой подход,

from numpy.random import seed, randint
from numpy import polymul, pad 
from numpy.fft import fft, ifft
from timeit import default_timer as timer

length=100

def test_mul(arr_a,arr_b):  #inbuilt python multiplication

    c=polymul(arr_a,arr_b)

    return c    

def sb_mul(arr_a,arr_b):    #my schoolbook multiplication

    c=[0]*(len(arr_a) + len(arr_b) - 1 )
    for i in range( len(arr_a) ):
        for j in range( len(arr_b) ):
            k=i+j
            c[k]=c[k]+arr_a[i]*arr_b[j]
    return c    


def fft_test(arr_a,arr_b):  #fft based polynomial multuplication

    arr_a1=pad(arr_a,(0,length),'constant')
    arr_b1=pad(arr_b,(0,length),'constant')
    a_f=fft(arr_a1)
    b_f=fft(arr_b1)

    c_f=[0]*(2*length)

    for i in range( len(a_f) ):
        c_f[i]=a_f[i]*b_f[i]

    return c_f


if __name__ == '__main__':

    seed(int(timer()))
    random=1
    if(random==1):
        x=randint(1,1000,length)
        y=randint(1,1000,length)
    else:
        x=[1]*length
        y=[1]*length

    start=timer()
    res=test_mul(x,y)
    end=timer()
    print("time for built in pol_mul", end-start)

    start=timer()
    res1=sb_mul(x,y)
    end=timer()
    print("time for schoolbook mult", end-start)

    res2=fft_test(x,y)

    print(res2)

    #########check############
    if( len(res)!=len(res1) ):
        print("ERROR");

    for i in range( len(res) ):
        if( res[i]!=res1[i] ):
            print("ERROR at pos ",i,"res[i]:",res[i],"res1[i]:",res1[i])

Теперь, вот мой подход в деталях, 1. Сначала я попробовал себя в наивной реализации Schoolbook со сложностью O (n ^ 2). Но, как вы можете ожидать, это оказалось очень медленно.

2. Во-вторых, я узнал polymul в библиотеке Numpy. Эта функция намного быстрее, чем предыдущая. Но я понял, что это тоже сложность O (n ^ 2). Вы можете видеть, что при увеличении длины k время увеличивается в k ^ 2 раза.

3. Мой третий подход заключается в том, чтобы попробовать умножение на основе БПФ с использованием встроенных функций БПФ. Я следовал хорошо известному подходу, также описанному здесь , но я не смог заставить его работать.

Теперь мои вопросы,

1. Где я ошибаюсь в подходе, основанном на БПФ? Подскажите, пожалуйста, как мне это исправить?

2. Является ли моё наблюдение, что функция polymul имеет сложность O (n ^ 2)?

Пожалуйста, дайте мне знать, если у вас есть какие-либо вопросы. Заранее спасибо.

Answer 1

1. Где я ошибаюсь в подходе, основанном на БПФ?Подскажите, пожалуйста, как мне это исправить?

Основная проблема заключается в том, что в подходе, основанном на FFT, вы должны выполнять обратное преобразование после умножения, но этот шаг отсутствуетиз вашего кода.С этим пропущенным шагом ваш код должен выглядеть следующим образом:

def fft_test(arr_a,arr_b):  #fft based polynomial multiplication

    arr_a1=pad(arr_a,(0,length),'constant')
    arr_b1=pad(arr_b,(0,length),'constant')
    a_f=fft(arr_a1)
    b_f=fft(arr_b1)

    c_f=[0]*(2*length)

    for i in range( len(a_f) ):
        c_f[i]=a_f[i]*b_f[i]

    return ifft(c_f)

Обратите внимание, что также может быть несколько возможностей для улучшений:

• Заполнение нулями может быть обработано напрямую путем передачитребуемая длина БПФ в качестве второго аргумента (например, a_f = fft(arr_a, length))
• Умножение коэффициента в вашем цикле for может напрямую обрабатываться numpy.multiply.
• Если коэффициенты полинома являются действительнымизатем вы можете использовать numpy.fft.rfft и numpy.fft.irfft (вместо numpy.fft.fft и numpy.fft.ifft) для некоторого дополнительного повышения производительности.

Таким образом, реализация для вещественных входных данных может выглядеть следующим образом:

from numpy.fft import rfft, rifft
from numpy import multiply
def fftrealpolymul(arr_a,arr_b):  #fft based real-valued polynomial multiplication

    L = len(arr_a1) + len(arr_b2) - 1
    a_f=rfft(arr_a,L)
    b_f=rfft(arr_b,L)

    return irfft(multiply(a_f,b_f))

Является ли мое наблюдение о том, что функция polymul имеет сложность O (n ² ), правильную?

Это также является производительностью, которую я наблюдаю, исоответствует доступному коду в моей простой установке (версия 1.15.4, и в более поздней версии 1.16.1 изменений в этой части не наблюдается).