11

Suppose the convolution of a general number of discrete probability density functions needs to be calculated. For the example below there are four distributions which take on values 0,1,2 with the specified probabilities:

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]])

The convolution can be found like this:

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

The probabilities of seeing 0,1,...,8 are then given by

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

This part is the bottleneck in my code and it seems there must be something available to vectorize this operation. Does anyone have a suggestion for making it faster?

Alternatively, a solution where you could use

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)

and get the pairwise convolutions

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

would also help tremendously.

Mark Dickinson
  • 29,088
  • 9
  • 83
  • 120
Forzaa
  • 1,465
  • 4
  • 15
  • 27
  • According to the numpy docs, the arguments to `np.convolve` can only be 1-dimensional. So I guess, there's not much to vectorize here. But maybe its worth to use a different convolution like scipy's fft based one? http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.fftconvolve.html – SmCaterpillar Mar 06 '15 at 14:48
  • @SmCaterpillar I played around with that a bit but my knowledge about convolutions is too limited to understand what's going on there. The version here I understand, but I haven't a clue how to specify the weights for the fft version. – Forzaa Mar 06 '15 at 14:52
  • What do you mean by weight? I tried both and both convolutions give the same result for your question. However, the fft one was much slower (due to overhead, your toy problem is too small, maybe when the pdfs themselves contain more values, you actually get a speed increase). – SmCaterpillar Mar 06 '15 at 14:58
  • @SmCaterpillar I suppose you are again using the for loop for the scipy version and convolute one by one. I'd like to avoid the for loop and apply the operation on all rows of pdfs immediately. – Forzaa Mar 06 '15 at 15:04
  • I was looking at this version of convolve for the record http://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.filters.convolve.html – Forzaa Mar 06 '15 at 15:13
  • You can parallelize it with multiprocessing. Also you have many 0 values. Skip convolution for them. – User Mar 06 '15 at 19:58
  • Can you comment on the typical sizes involved? What's the typical shape of the `pdfs` array? – Mark Dickinson Mar 24 '15 at 09:49
  • @Mark There could be 10 of these pdfs with 100 entries each. – Forzaa Mar 24 '15 at 10:01

1 Answers1

20

You can compute the convolution of all your PDFs efficiently using fast fourier transforms (FFTs): the key fact is that the FFT of the convolution is the product of the FFTs of the individual probability density functions. So transform each PDF, multiply the transformed PDFs together, and then perform the inverse transform. You'll need to pad each input PDF with zeros to the appropriate length to avoid effects from wraparound.

This should be reasonably efficient: if you have m PDFs, each containing n entries, then the time to compute the convolution using this method should grow as (m^2)n log(mn). The time is dominated by the FFTs, and we're effectively computing m + 1 independent FFTs (m forward transforms and one inverse transform), each of an array of length no greater than mn. But as always, if you want real timings you should profile.

Here's some code:

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

Applying this to your example, here's what I get:

>>> 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])

That's the basic idea. If you want to tweak this, you might also look at numpy.fft.rfft (and its inverse, numpy.fft.irfft), which take advantage of the fact that the input is real to produce more compact transformed arrays. You might also be able to gain some speed by padding the rows array with zeros so that the total number of columns is optimal for performing an FFT. The definition of "optimal" here would depend on the FFT implementation, but powers of two would be good targets, for example. Finally, there are some obvious simplifications that can be made when creating rows if all the input arrays have the same length. But I'll leave these potential enhancements to you.

Mark Dickinson
  • 29,088
  • 9
  • 83
  • 120
  • Why not use ``scipy.signal.fftconvolve()`` (http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.fftconvolve.html)? – Dietrich Mar 24 '15 at 20:39
  • 2
    @Dietrich: Because (unless I'm missing something) that only convolves two arrays at a time, and using it repeatedly would involve a lot of unnecessary transforming and untransforming. – Mark Dickinson Mar 24 '15 at 20:42
  • @MarkDickinson Could you elaborate how we can match the output (the density probabilities) to the actual outcomes? That is how do we calculate the outcomes to which these probabilities belong? What is the purpose of `result_length`? Why are we adding some zeros to each array, which we never fill anyway since we fill the rows array only till `:len(array)`?. – user2974951 Aug 18 '20 at 09:01
  • 1
    @user2974951 It depends a bit on what you're doing. Typically, you're convolving the PDFs of random variables `X_1`, `X_2`, ..., `X_n` in order to get the PDF of `X_1 + X_2 + ... + X_n`. In that case, each `X_i` needs to be discrete with the possible values evenly spaced using some spacing `s`, say, and that spacing needs to match for all of the `X_i`. And then the outcomes correspond to evenly-spaced values with spacing `s`, starting at the sum of the minima of each `X_i` and going to the sum of the maxima of each `X_i`. – Mark Dickinson Aug 21 '20 at 11:28
  • @user2974951 `result_length` is, as the name suggests, the length of the resulting convolution array. We need to pad each input array to that length so that the FFTs are compatible. – Mark Dickinson Aug 21 '20 at 11:30