Does there exists 1D FFT/IFFT implementation in pure Julia (i.e. no FFTW)

Specific Domains Numerics
1

Hi,
I need a pure julia FFT implementation as I am planning to use it in correspondence with the ForwardDiff package. Does this exist?

2

You could adapt the BigFloat fft implementation in FastTransforms.jl: https://github.com/MikaelSlevinsky/FastTransforms.jl/blob/master/src/fftBigFloat.jl

1 Like
3

I implemented an experimental one in: https://github.com/JuliaLang/julia/pull/6193

1 Like
4

Hmmm, it seems I will be running into problems with forwarddiffs dual numbers and the resulting complex numbers from the FFT.

5

Complex numbers are indeed not supported by ForwardDiff. If they were, a better way to differentiate through FFT would be to teach ForwardDiff the derivatives of FFTs: https://github.com/JuliaDiff/ForwardDiff.jl/pull/165 and refs

2 Likes
6

Yes, since it is a linear transformation, the derivative of an FFT with respect to its inputs is simply another FFT.

7

For a vector x ou can always describe its Fourier Transform as y = D x where D is the DFT Matrix.

You can even generalize it for cases which are “Zero Padded”.
Since the Derivative is just {D}^{T} all you need is a function which generates those matrices.

Is there such built in function for that in Julia?

8

Okey, so what I am really after is a way to find the derivative of the Volterra process \int_0^tg_\alpha(t-s) dW_s wrt. \alpha \in (-0.5,0) where g_\alpha(t-s) = (t-s)^\alpha and W a brownian motion. I am approximating this integral by looking at at the Riemann sum approximation \tilde{Y}(\frac{i}{n})= \sum_{k =1}^i g(\frac{b_k}{n})W_{i-k} on an n sized grid, i = 0,1,2,\dots. This can be interpreted as a discrete convolution \tilde{Y}(\frac{i}{n}) = (g*W)_i, which can be computed efficiently by FFT (Ihave gotten this to work with ForwardDiff using a direct naive convolution, but its too slow). Furthermore, if I need to do several runs of this numerical scheme to find an estimate of \mathbb{E}[\tilde{Y}(\frac{i}{n})] \sim \frac{1}{N}\sum_j^N\tilde{Y}^{(j)}(\frac{i}{n}). I have to perform this convolution each time I simulate a path of this process. So to speed up things I would like to compute the FFT of g(b_k/n) once (since its deterministic) and multiply it with the FFT of a new set of random variates W, and then take the IFFT to arrive at the convolution.

The code to my current implementation (with some slight changes to the actual approximation scheme):

function gen_volterra(a,s,n,G,rand_1,rand_2)
  Z_1 = zeros(eltype(a),s); #first normal
  Z_2 = zeros(eltype(a),s); # 2nd normal
  #G = zeros(eltype(a),s+1); #Kernel function
  @fastmath @inbounds @simd for i = 1:s
    #G[i+1] = b(a,i+1)/(n^a); #fill kernel vector <<----this is where the error occurs
    Z_1[i] = get_sigma_Z1(n)*rand_1[i] #set correct variance of random variate
    Z_2[i] = get_sigma_Z2(a,n)*(get_rho(a)*rand_1[i]+sqrt(1-get_rho(a)*get_rho(a))*rand_2[i]); #set correct variance of random variate
  end
  return [0;(direct_conv(Z_1,G)[1:s]+Z_2)]; #Convolute
  #return G
end

function direct_conv{T}(a::Array{T}, b::Array{T})
    m = length(a)
    n = length(b)
    c = zeros(T,m+n-1)
    @fastmath @inbounds @simd for j=1:m
        @fastmath @inbounds @simd for k=1:n
            c[j+k-1] += a[j]*b[k]
        end
    end
    return c
end

and addtional helper functions:

function b(a, k)
  return (k^(a+1.0)-(k-1.0)^(a+1.0))/(a+1.0);
end

function get_sigma_Z1(n)
  return 1.0/sqrt(n);
end

function get_sigma_Z2(a,n)
  return 1.0/sqrt(((2.0*a+1.0)*n^(2.0*a+1.0)))
end

function get_rho(a)
  return sqrt(2.0*a+1.0)/(a+1.0)
end

function gen_G(a,n,s)
  i = range(1,s)
  return [0;map(i->b(a,i+1)/(n^a),i)]
end

Example of generating 1 path:

julia> a= -0.43
julia> n= 500
julia> s= 500
julia> G = gen_G(a,n,s)
julia> rand_1 = randn(s)
julia> rand_2 = randn(s)
julia> path = gen_volterra(a,s,n,G,rand_1,rand_2)

But the operation below will then be quite slow due to the number of floating point operations in the direct_conv method.

m = ForwardDiff.derivative(a ->gen_volterra(a,s,n,G,rand_1,rand_2),-0.43)
1 Like
9

You can just FFT the columns of the identity matrix. But usually it is preferable to work with the DFT matrix implicitly, since multiplying by it (or its transpose — it is symmetric) is computed efficiently by the fft function (and the inverse or conjugate by ifft).

10

Why do you need AD for this? Just do the derivative wrt alpha, which is again a convolution that you can do with an FFT.

11

Because this process is used as an input in a SDE for another process, which I didn’t include here for clairity. It is actually wrt to this other process I want the gradient, so the stuff outlined above is a subproblem as this is where the performance bottleneck is. But a solution might be to use this