[SciPy-Dev] Backend support for fftpack - GSoC 2019

[Ralf Gommers]

On Thu, May 9, 2019 at 2:25 PM Peter Bell <PeterBell10 at live.co.uk> wrote:

> >> Should backends be allowed to add functionality or should they be
> completely
>
> >> interchangeable. E.g. unlike fftpack, FFTW has support for long doubles;
>
> >> should this be exposed to SciPy users?
>
> >
>
> > If a user writes backend-specific code, then they may just as well
> bypass the
>
> > backend right, and call the other library directly?
>

After some thought/discussion I'm going to disagree with myself here.


>
> I agree, in fact I quote you in the proposal saying something similar.
>
>
>
> > Although for long double things will work out of the box perhaps,
> because the
>
> > dtype of an input array doesn't come back in the API. So the array can be
>
> > passed straight to the backend.
>
>
>
> This was why I picked it as an example. It doesn't add any new parameters
> or
>
> change the function signature at all, but it does change the result of a
>
> particular call to the API. So it's a very minimal change yet it still
> breaks
>
> the interchangeability of backends.
>

Peter and I had a very interesting discussion in comments on his Google
doc. It relates to the design question about what a backend is allowed to
do (and even higher level, what is a backend for). I'm now relatively sure
I can summarize it well, so I'll give it a try.

"completely interchangeable", as Peter uses above, implies that both API
and semantics of the backend are matching 100% with what SciPy does.
Peter's example from the doc (where _dct_func is the backend-provided
function) may help to illustrate that:

   def dct(x, type=2, n=None, axis=-1, norm=None, overwrite_x=False):

       # From current implementation, also converts float16 to float32

       x = _asfarray(x)
       if x.dtype not in (np.float32, np.float64,

                          np.complex64, np.complex128):
           raise ValueError("type {} is not supported".format(x.dtype))
        return _dct_func(x, type, n, axis, norm, overwrite_x)

So that will do for example array coercion and upcasting the same way as
the SciPy implementation. As a consequence, what a backend offers will then
be different performance - speed and/or accuracy. Let's call this *model 1*.

What my mental model was when I started commenting on the above example is:
a backend has to have a compatible API, but replaces the implementation of
a function completely. Hence can offer different behavior. The downside is
that a user can't blindly switch backends, reading the docs may be
required. The upside is that it allows for more use cases, like:
1. Simple differences, like providing a longdouble implementation
2. Another input type, e.g. cupy arrays that live on a GPU and executing
the fft on the GPU. With the "completely interchangeable model" this may
still work,  but there will be a conversion to-from ndarray on every
function call.
3. Another example of another input type that otherwise won't work at all:
using things that don't fit in memory, like a large Dask array.
Let's call this *model 2*.

I would either way not require _completely_ identical behavior, even if
we'd go with model 1. For example the overwrite_x keyword is an
implementation detail, I think it's fine for any backend to raise an error
if that's True.

Model 2, which I have a preference for, is very similar to the numpy
__array_function__ and __array_ufunc__ protocols (in spirit, not
implementation). I've tried to illustrate that in slides 10-11 of [1].

Both models 1 and 2 allow faster implementations, like mkl-fft and pyfftw,
to cleanly integrate with SciPy (instead of monkeypatching like pyfftw does
now). Model 2 in addition allows writing code that is generic for (for
example) GPU and distributed implementations, at the cost of a little less
uniformity in behavior.

I hope this made sense. Are there other models of a backend than the above
two? Thoughts on what to prefer?

Cheers,
Ralf

[1]
https://www.slideshare.net/RalfGommers/the-evolution-of-array-computing-in-python
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20190517/17661679/attachment-0001.html>