Более быстрая свертка функций плотности вероятности в Python

Question

Предположим, что необходимо вычислить свертку общего числа функций плотности дискретной вероятности. Для примера ниже есть четыре распределения, которые принимают значения 0,1,2 с заданными вероятностями:

import numpy as np 
pdfs = np.array([[0.6,0.3,0.1],[0.5,0.4,0.1],[0.3,0.7,0.0],[1.0,0.0,0.0]]) 

свертка можно найти так:

pdf = pdfs[0]   
for i in range(1,pdfs.shape[0]): 
    pdf = np.convolve(pdfs[i], pdf) 

Вероятности видеть 0, 1, ..., 8 затем дается

array([ 0.09 , 0.327, 0.342, 0.182, 0.052, 0.007, 0. , 0. , 0. ]) 

Эта часть является узким местом в моем коде, и кажется, что должно быть что-то доступно для векторизации этой операции. Есть ли у кого-нибудь предложение сделать это быстрее?

В качестве альтернативы, решение, в котором вы могли бы использовать

pdf1 = np.array([[0.6,0.3,0.1],[0.5,0.4,0.1]]) 
pdf2 = np.array([[0.3,0.7,0.0],[1.0,0.0,0.0]]) 
convolve(pd1,pd2) 

и получить попарных сверток

array([[ 0.18, 0.51, 0.24, 0.07, 0. ], 
     [ 0.5, 0.4, 0.1, 0. , 0. ]]) 

также поможет чрезвычайно.

Answer 1

Вы можете эффективно вычислить свертку всех ваших PDF-файлов с использованием быстрых преобразований Фурье (FFT): ключевым фактом является то, что FFT of the convolution является продуктом БПФ отдельных функций плотности вероятности. Таким образом, преобразуйте каждый PDF, соедините преобразованные PDF-файлы вместе, а затем выполните обратное преобразование. Вам нужно будет поместить каждый входной PDF с нулями на соответствующую длину, чтобы избежать эффектов от wraparound.

Это должно быть достаточно эффективным: если у вас есть m PDF-файлов, каждый из которых содержит n записи, то время, чтобы вычислить свертку, используя этот метод должен расти, как (m^2)n log(mn). Время доминирует БПФ, и мы эффективно вычисляем независимые БПФ m и одно обратное преобразование, каждый из массива длиной не более mn. Но, как всегда, если вам нужны реальные тайминги, вы должны профиль.

Вот код:

import numpy.fft 

def convolve_many(arrays): 
    """ 
    Convolve a list of 1d float arrays together, using FFTs. 
    The arrays need not have the same length, but each array should 
    have length at least 1. 

    """ 
    result_length = 1 + sum((len(array) - 1) for array in arrays) 

    # Copy each array into a 2d array of the appropriate shape. 
    rows = numpy.zeros((len(arrays), result_length)) 
    for i, array in enumerate(arrays): 
     rows[i, :len(array)] = array 

    # Transform, take the product, and do the inverse transform 
    # to get the convolution. 
    fft_of_rows = numpy.fft.fft(rows) 
    fft_of_convolution = fft_of_rows.prod(axis=0) 
    convolution = numpy.fft.ifft(fft_of_convolution) 

    # Assuming real inputs, the imaginary part of the output can 
    # be ignored. 
    return convolution.real 

Применяя это к вашему примеру, вот что я получаю:

>>> convolve_many([[0.6, 0.3, 0.1], [0.5, 0.4, 0.1], [0.3, 0.7], [1.0]]) 
array([ 0.09 , 0.327, 0.342, 0.182, 0.052, 0.007]) 

Это основная идея. Если вы хотите изменить это, вы также можете посмотреть на numpy.fft.rfft (и его обратный, numpy.fft.irfft), которые используют тот факт, что вход вещественный для создания более компактных преобразованных массивов. Вы также можете получить некоторую скорость, заполнив массив rows нулями, чтобы общее количество столбцов было оптимальным для выполнения БПФ. Определение «оптимального» здесь будет зависеть от реализации БПФ, но, например, силы двух будут хорошими целями. Наконец, есть некоторые очевидные упрощения, которые могут быть сделаны при создании rows, если все входные массивы имеют одинаковую длину. Но я оставлю эти потенциальные улучшения для вас.