[SciPy-Dev] PR for overlap-add convolution

[Todd]

I have submitted a PR [1] for a new function, "oaconvolve", which is an
implementation of the overlap-add algorithm for convolution.  This
algorithm is much faster than conventional FFT-based convolution when both
inputs are large and one is much bigger than the other.  I would be
interested in getting your feedback.  Any improvements on the algorithm
would be greatly appreciated, and I am open to any suggestions on any
aspect of the PR.

My plans were previously discussed on the mailing list [2].  The original
request for this feature dates back to 2009 [3].

More details below for those interested:

The function supports all arguments fftconvolve supports, including working
on an arbitrary subset of axes of multidimensional arrays.  It will also
automatically fall back on using fftconvolve in situations where using the
overlap-add method doesn't apply (such as when the two inputs are the same
shape).

This implementation works in a vectorized manner.  The input arrays are
reshaped so the frequency-domain convolution of all overlapping blocks can
be calculated in a single step rather than using a loop like most
implementations I have seen.  This means that it should work efficiently on
any FFT backend that supports vectorization, unlike loops that would be
slowed down by FFT backends with long plan creation times.  Currently only
pocketfft is supported by scipy, but there are plans to support additional
FFT backends [4].  This approach should just work for any of them.
pocketfft has very low plan creation time so this is less of an issue, but
it also supports vectorizing multidimensional FFTs so this approach is
still highly efficient for pocketfft.

Some potential issues:

We previously discussed whether this should be its own function or part of
fftconvolve.  The performance and memory characteristics of this function
are significantly different than those of fftconvolve and would really be
used in different situations, so the consensus seemed to be that this
should be its own function.  I would ultimately like "convolve" to be able
to use oaconvolve when it is the best choice, but that is an issue for
later.

There is a "method" argument for choosing whether to automatically fall
back on fftconvolve or raise an error.  Being able to prevent falling back
on fftconvolve is essentially for testing.  We need to be able to make sure
that oaconvolve is actually being tested and not the fftconvolve fallback.
But this may not be the best way to do that.

I have abstracted much of the shared functionality between oaconvolve and
fftconvolve to several hidden functions, which should make it easier to add
shared features to both.  But if this is considered too complicated I can
make the functions independent.

The tests are slow.  There are many potential corner cases with particular
combinations of sizes so I have a test marked slow to check a large variety
of sizes. This means it is disabled in some situations, but it does need to
be run at least sometimes.  This is only for one dimension though.
Checking more this way would likely be prohibitively slow.

There is a detailed docstring, but I didn't add anything to the tutorial.
Should I?

Are there any additional references anyone things I should add?

Is the example okay?  Is there another one someone thinks would be useful?

I included a bugfix that this code depends on.  Should I split it into
another pull request or is it okay having them together?  The bugfix has to
be accepted before the new function can be accepted since oaconvolve will
not work without it.


[1] https://github.com/scipy/scipy/pull/10869
[2] https://mail.python.org/pipermail/scipy-dev/2019-June/023555.html
[3] https://github.com/scipy/scipy/issues/1364
[4] https://mail.python.org/pipermail/scipy-dev/2019-May/023506.html
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20190926/6822ce9b/attachment.html>