Proposal: faster cwt() with FFT based convolution

[alsauve]

Hi there,

The cwt() code currently use a discrete convolution algorithm which is ~O( len(data) * scale).
While this is usually fine for common values of data and scale < 100, it becomes an issue quickly for EEG like signals with 100k samples or sound.

To overcome this limit I implemented a selective algorithm which switch to
~O( len(data) * log(len(data)) with a FFT based convolution when its complexity is lower.

I propose to merge this version with the pywt tree if you guys find it appropriate.

The source code has been implemented as a drop-in replacement to the function cwt() and is available on the scaleogram project as fastcwt() in this file.

The output has been unit tested with all continuous wavelet to be the same as the current pywt.cwt() output at floating point precision level, which is expected from FFT based convolution when circular convolution is avoided.

Example::

    %time (coef1, freq1) = fastcwt(np.arange(140000), np.arange(2,200), 'cmorl1-1')
    => CPU times: user 12.6 s, sys: 2.2 s, total: 14.8 s
    => Wall time: 14.9 s

    %time (coef1, freq1) = pywt.cwt(np.arange(140000), np.arange(2,200), 'cmorl1-1')
    => CPU times: user 1min 50s, sys: 401 ms, total: 1min 51s
    => Wall time: 1min 51s

This file is part of a a visualization tool designed for wavelet beginners and is wrapped around PyWavelets.

[grlee77]

Hi @alsauve

I think this would be a welcome addition. There is a clear performance benefit and it doesn't look like any new dependencies would need to be added. Out of curiosity is your benchmark above using a NumPy configured to use Intel's mkl_fft (i.e. as in the default NumPy builds provided by Anaconda?). If so then the FFT is probably multi-threaded already, but if not then it should be even faster with a multi-threaded FFT.

The output has been unit tested with all continuous wavelet to be the same as the current pywt.cwt() output at floating point precision level.

Great. Since the output is compatible with the existing function, it is probably easiest from a user standpoint to just keep the current cwt function name in PyWavelets rather than add a new function. Just configure the existing cwt to call the FFT variant internally when warranted. I'm not yet sure about the best API choice, but one could potentially allow the user to override the default choice to explicitly choose the FFT or convolution-based approach.

Also, because the FFT has axis support (which is not present in np.convolve), it should also be pretty easy to extend this to add axis support for batch application of 1D CWT as requested in #445.

One other thought: If I recall correctly, NumPy's FFTs always use double precision (although this is no longer true when NumPy is using mkl_fft internally)? For backwards compatiblity with existing user code we should ensure that single precision inputs still give single precision outputs even if the FFT itself ends up using double precision.

[alsauve]

Hi @grlee77 ,

I'm not yet sure about the best API choice, but one could potentially allow the user to override the default choice to explicitly choose the FFT or convolution-based approach.
This is an interresting question as I was not aware about mkl_fft, I checked my installation is fairly standard from pip3 on ubuntu 18.04

Actually it is implemented this way in scaleogram (at least the bases are there), this is a good practice to have choice to revert to working default.

Also, because the FFT has axis support (which is not present in np.convolve), it should also be pretty easy to extend this to add axis support for batch application of 1D CWT as requested in #445.

Seems feasible but never used or seen any 2D CWT.

One other thought: If I recall correctly, NumPy's FFTs always use double precision (although this is no longer true when NumPy is using mkl_fft internally)? For backwards compatiblity with existing user code we should ensure that single precision inputs still give single precision outputs even if the FFT itself ends up using double precision.

Did not knew about mkl-fft for numpy
...trying....
Well it crashes on my Core2 Duo with either python 2.7 and 3.6 ... postponing this one ;)

For double precision: attention is payed to the dtype value in fastcwt() but results indicates that both versions have issues:

  # pywt
  In [6]: sig = np.ones(100000, dtype=np.float128) # <=== 128 here
  In [7]: c1,f1 = pywt.cwt(sig, np.arange(200,201), 'mexh')
  In [8]: c1.dtype
  Out[8]: dtype('float64') # <==== 64 here :-/
  
  In [9]: import scaleogram as scg
  In [10]: c2,f2 = scg.fastcwt(sig, np.arange(200,201), 'mexh')
  In [11]: c2.dtype
  Out[11]: dtype('float64') # <==== 64 here too :-/

Thanks for pointing that, I found the missing bit in my function and corrected it for homogeneous input/output, but the case of complex wavelet is more tricky to implement,
one have to promote the input real dtype on complex type to the same precision into complex.
Don't know how to do that in a simple way.

[grlee77]

Are you interested in starting a PR with what you have? You do not need to worry about batch (>1D) support in the initial PR and any new tests would most likely go in the file pywt/tests/test_cwt_wavelets.py. I can help fix up the dtype conversions, etc. to match PyWavelets conventions.

In regards to precision, perhaps I was thinking of how all the DWT routines handle dtypes. It looks like the output array for cwt is always initialized in double-precision. That inconsistency with DWT behavior is a separate issue that we may want to address in a future release. (Those discrete transforms preserve precision for single or double precision and any other dtypes get converted to an appropriate floating point type)

[alsauve]

Sure, I'll be glad to provide a PR for 1D convolutions and add relevant unit tests.
Would this kind of design be appropriate to you?

cwt_legacy();
cwt_fft_based();
_CALL_CWT = cwt_fft_based
set_cwt_fun(callable);
cwt(*args, **kw):
   return  _CALL_CWT(*args, **kw)

For precision, currently the dtype output matching is fixed for real values, and just while writing this I found an elegant way of promoting a real array to complex using the same precision : multiplication by 1j

In [38]: for dtype_test in [np.float16, np.float32, np.float64, np.float128]: 
    ...:     print( dtype_test, '=>', (np.zeros(1, dtype=dtype_test)*1j).dtype)
<class 'numpy.float16'> => complex64
<class 'numpy.float32'> => complex64
<class 'numpy.float64'> => complex128
<class 'numpy.float128'> => complex256

The dtype np.complex32 matching np.float16 does not appear to exist within numpy 1.13

[grlee77]

Okay, great to hear.

Let's just keep the dtype handling consistent with the current CWT for now (i.e. just using float64). Preserving single precision can be added later in a separate PR as that may need a deprecation cycle to avoid unexpected changes for users.

For API, I am not a big fan of having _CALL_CWT or set_cwt_fun methods. I would just append an additional keyword argument with a name such as method or implementation to the existing cwt function where the values of the keyword could be one of {'auto', 'convolution', 'fft'}. Then a simple if/else based on the value can be used to select the internal implementation. If auto is specified some heuristic such as you proposed above could be used to auto-select the faster approach for the given data size and scales.

Please proceed to make a PR and we can refine further based on the implementation there.

[alsauve]

Ok, fine with the roadmap with a preservation of function output type for user side stability and a keyword argument for processing mode selection. Thinking again about it, the keword argument is homogeneous with many others library, so this is just good practice for the end user experience.
Let do that.

[grlee77]

closed by #490