FFT design considerations

Hi!

Having an FFT in Stan would be really cool!

One of the very cool applications of a FFT transform (and its inverse) is the possibility to calculate the convolution as a simple product of two functions in the transformed space. This is for example handy if one wants to get the difference distribution of two quantities; one simply FFT transforms the function, multiplies them and then takes the inverse of that product.

Long story short: If this is possible with what is planned then by all means it should be taken in consideration for the design. That would be quite useful I think.

Best,
Sebastian

1 Like

It’s very inconvenient to do this without complex numbers. I think the options here are

  1. Implement a complex type and vectorized complex arithmetic.
  2. Establish a convention for treating 2-by-n matrices as complex vectors, providing complex arithmetic.
  3. Not expose any complex arithmetic, but give users a convolution function.

I think 1 has the character of tilting at windmills for pretty thin benefit, it’s a truly huge project. 2 seems like a design disaster (picturing nightmares if our handling of matrices evolves or eigen’s does). We can’t prevent users from doing this in an ad hoc way but I don’t want to bless it. 3 is a bit unsatisfying but has the benefit of letting us add the pieces one at a time. If after a few months we have a zoo of signal processing functions, we may want to re-evaluate this.

Thel

Complex type isn’t enough of a priority now. And I
think the derivatives are pretty tricky to do in general.

Is the second option what we were talking about before?
If so, I don’t see a problem with that.

I don’t understand the third option. Is the convolution
why people want these operations? Or is it more things
which is why you think libraries may grow?

  • Bob

Forget convolution for now – there are more subtle issues involved due to the discretization assumed in the DFT. Right now let’s just start with

real[] fft_amps(real[] x);
real[] fft_phases(real[] x);
real[,] fft_amps_phases(real[] x);
2 Likes

With regards to 2 I’m just saying that I don’t want to give people a bunch of functions like:

real[2,] pointwise_multiply_as_if_complex_vector(real[2,-] z, real[2,-] w)

as a brittle workaround for the lack of a complex type. People may slice the thing we give them and do all kinds of manipulation, I’m just kind of wary of creating a pseudo-type. So I’m saying that we shouldn’t really worry about convolution for now, but later if there’s a need we can specific functions for convolution.

Michael’s signatures are exactly right, though I’m tempted to add one more: real[,] fft_real_imag(real[] x);

I’d suggest starting with just one function, getting
feedback and merging, then going on to the rest.

  • Bob

That’s absolutely the plan. Just figured Sebastian deserved a response on convolutions.

Hi!

Ok; let’s keep things simple. The Fourier transform is super handy for a convolution, but I am not an expert on the discrete FFT nor is this an easy request as this entails to transform and back-transform. Still, if we think that convolutions are feasible at some point, then we should probably design the internal code with this application in mind? But as Bob says, getting started with the simple function is a great start.

Sebastian

The actual FFT code is sitting inside of Eigen, so it’s
out of our hands unless we want to rewrite it.

  • Bob

Oh, did not know that. This makes it much easier…just expose what they did sounds as the way to go then.

That’s the issue—how to do that. Their API is based
on complex numbers, which we don’t have in Stan.

  • Bob

Hi Bob,

I am working on using convolution operators to model spatial population dynamics. Right now my go-to solution is to use nested loops in the function block to implement the convolution. Ideally, the simulations can be made much faster by solving convolution by FFT.

I saw that the Stan team just released support for complex numbers but I could not find a FFT (iFFT) function. Is the stan team working on providing FFT related functions to solve convolutions? thank you in advance

1 Like