New to julia: Convolution and juliadiff

Specific Domains Numerics
1

Hi,

I am sorry if this is the wrong place to ask this. But I am considering learning Julia for a stochastic analysis project I am working on. I have a numerical scheme to approximately simulate the stochastic integral Y(t) = \int_0^tg(t-s)dW_s where g(t-s) = (t-s)^\alpha,\quad \alpha \in (-\frac{1}{2},0). W is a standard brownian motion. The scheme approximates the integral using the riemann sum approximation \tilde{Y}(\frac{i}{n})= \sum_{k =1}^i g(\frac{i}{n})W_{i-k} on an n sized grid. This can be interpreted as a discrete convolution \tilde{Y}(\frac{i}{n}) = (g*W)_i. So I can run simulations for Y by drawing normal variates and computing the discrete convolution. But what I am effectively seeking is to compute the derivative of the sample mean \mathbb{E}[Y(\frac{i}{n})] = \frac{1}{N}\sum_{j=0}^N \tilde{Y_j}(\frac{i}{n}) with respect to the parameter \alpha. So my idea was to use julia forward diff to do this. But I am not sure how this works together with the convolution routine in julia. Anyone who can give me any guidance on whether this can work or not?

2

If it’s all Julia code it should work. However, this uses the FFT, which I believe will call out to FFTW, a C-code. I think this will block autodifferentiation. If that’s the case then you’ll need to find a native Julia FFT package to them implement the convolution yourself (but conv is written in Julia, so it should be clear how to do that). But before doing that, give it a try: autodifferentiation “just works” on many Julia functions via generic fallbacks already provided in Base.

Also, take a look at Bridge.jl for inspiration on stochastic analysis tools, and StochasticDiffEq.jl for inspiration if it’s for SDEs.

3

@ChrisRackauckas yeah, your instinct is right: I tried this myself and found that conv doesn’t work for Vector{ForwardDiff.Dual}.

@vgdev all of Chris’s advice is good, but if you don’t need the benefits of FFT-powered convolution, you could also just write the basic convolution routine yourself in Julia, and forward diff will then magically work with your convolution function.

4

Maybe the convolution is just a red herring. If my algebra is correct, then \tilde{Y}(\frac{i}{n}) = \tilde{Y}(\frac{i-1}{n}) + g(\frac{i}{n})W_{i} for each sample of W. So, sampling the random walk N times and collecting the sums of Y(\cdot) at each point \frac{i}{n} will be efficient. The gain of using the convolution is (with FFT trickery) getting it to n\log(n) operations over the naive n^2 for computing the sums \tilde{Y}(\cdot) independently for each i. The cumulative sum approach has just order n steps.

5

Thanks @ChrisRackauckas and @rdeits ! Will there be any slowdown in the computing time of the convolution routine when using forwarddiff? I tried something similar using the templatising approach with FADBAD++ in C++ with my own convolution function. As the convolutions involves a lot of (although simple) operations, this resulted in the convolution routine being significantly slower with the FADBAD double type, than with the standard double.

@rdeits the size of my two 1D vectors will probably be in the range of 500. According to this I will prob not get any advantage of the FFT routine. Also can you please share the example you worked on, so I can see how to go about this?

@Dan , it should indeed be done by convolution. I probably simplified the example too much and skipped some notation. But if interested you can have a look at the hybrid scheme here (pp 15-16)

6

Ahh you are right. My algebra was indeed mistaken. But something was weird in the question, since g(\frac{i}{n}) is constant in the sum and can be taken out. The convolution would need it to be g(\frac{k}{n}) ??

7

ForwardDiff.jl has been really well written with a focus on performance, but computing the derivative will always involve somewhat more work than just computing the value of the function. It’s really hard to say how much the slowdown will be, but if your function is scalar → scalar, then I would guess that ForwardDiff.derivative() should be about 2X slower than just calling the function normally (and when taking gradients of functions with vector inputs, forward-mode autodiff will get linearly slower in the number of inputs to your function).

I would encourage you to just try it out and use the BenchmarkTools.jl package to see the performance impact.

As for the example I worked on, all I did was try a very simple function which included conv():

julia> f(α) = sum(conv(α * [5., 6, 7, 8], [1., 2, 3, 4]))
f (generic function with 1 method)

julia> ForwardDiff.derivative(f, 1.0)
ERROR: MethodError: no method matching conv(::Array{ForwardDiff.Dual{1,Float64},1}, ::Array{Float64,1})

which showed that conv() wasn’t implemented for the special Dual type that ForwardDiff uses. On the other hand, a naive implementation of conv() should be pretty easy to write, and should “just work” with Dual inputs.

For example, let’s write a function (not a convolution, but something simpler and representative):

function foo(x, y)
   x' * y
end

And maybe the function you’re interested in takes data vectors x and y and a parameter a:

function bar(a, x, y)
  foo(a .* x, y)
end

Then, for some given x and y, we can easily take the derivative of bar w.r.t. a, evaluated at a = 1.0:

julia> x = rand(3)
julia> y = rand(3)
julia> ForwardDiff.derivative(a -> bar(a, x, y), 1.0)
1.018944053791661

(note the use of an anonymous function to create a new function of one argument (a) that has the current values of x and y “baked in”).

We can check performance easily:

julia> @benchmark bar($a, $x, $y)
BenchmarkTools.Trial:
  memory estimate:  112 bytes
  allocs estimate:  1
  --------------
  minimum time:     50.699 ns (0.00% GC)
  median time:      52.859 ns (0.00% GC)
  mean time:        58.809 ns (7.41% GC)
  maximum time:     1.238 ÎĽs (86.53% GC)
  --------------
  samples:          10000
  evals/sample:     987

julia> @benchmark ForwardDiff.derivative(a -> bar(a, $x, $y), $a)
BenchmarkTools.Trial:
  memory estimate:  272 bytes
  allocs estimate:  6
  --------------
  minimum time:     136.237 ns (0.00% GC)
  median time:      141.569 ns (0.00% GC)
  mean time:        167.207 ns (12.90% GC)
  maximum time:     3.976 ÎĽs (89.75% GC)
  --------------
  samples:          10000
  evals/sample:     873

so the derivative is ~2.5x slower, which is pretty close to my wild guess from earlier. Fortunately, that slowdown factor should not depend much on the number of computations you do, only on the number of inputs to your function.

1 Like
8

Even for shorter vectors you can get a speedup, since one of the vectors you are convolving with is constant. Computing the FFT of this vector once, and just scalar multiplying with an FFT of a random walk (basically half the work) and doing the reverse FFT.

In any case, there may be a way to sample the FFT of a random walk directly (sounds like something mathematicians would like to calculate) and then there is even more to gain.

1 Like
9

@dan yes, sorry, it depends on k too. There is an optimal discretisation b_k that can be made, so it should be somthing like g(\frac{b_k}{n})

10

@rdeits thanks a lot for the useful examples! So, an expected slowdown of roughly 2x then.
Given that I probably need to have the possibility of doing N =1,000,000 simulations I have in my current C++ implementation carried out the path generations in parallel, to make this go as fast as possible. I.e. if I have 4 cores, each core does ~250,000 convolutions each. Is there any way I can carry out the convolutions in parallel when also using forwardDiff in Julia?

@Dan Do you mean “for even longer vectors” ? Thanks for pointing out that one of the vectors are indeed constant, and its FFT can thus be precomputed.

11

Hi guys. So I have been working on this. But I am struggling to get it to run. (keep in mind that I need to make this run as quickly as possible). The function carrying out the numerical scheme below, gives me an error when trying to get the gradient at the line where i fill the kernel vector using the b function. How should I solve this?

@time ForwardDiff.derivative(a ->gen_volterra_for(a,s,n,rand_1,rand_2), -0.43)

MethodError: Cannot convert an object of type ForwardDiff.Dual{1,Float64} to an object of type Int64
This may have arisen from a call to the constructor Int64(…),
since type constructors fall back to convert methods.
in derivative at ForwardDiff\src\derivative.jl:5
in gen_volterra_for at volterra:67strong text

function gen_volterra_for(a,s,n,rand_1,rand_2)
  G = zeros(eltype(s),s+1); #Kernel function
  Z_1 = zeros(eltype(a),s); #first normal
  Z_2 = zeros(eltype(a),s+1); # 2nd normal
   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+1] = get_sigma_Z2(a,n)*(get_rho(a)*rand_1[i]+sqrt(1-get_rho(a)^2.0)*rand_2[i]); #set correct variance of random variate
  end
  return (direct_conv(Z_1,G)[1:s+1]+Z_2); #Convolute and return
end

Various other utility 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 sqrt(1.0/n);
end

function get_sigma_Z2(a,n)
  return sqrt(1.0/((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 direct_conv{T}(a::Array{T}, b::Array{T})
    m = length(a)
    n = length(b)
    c = zeros(T,m+n-1)
     @inbounds @simd for j=1:m
        @inbounds @simd for k=1:n
            c[j+k-1] += a[j]*b[k]
        end
    end
    return c
end
12

Did you mean eltype(s) ? Because s is probably an Int and the RHS of G[...] isn’t

1 Like
13

Obviously @Dan. The setup now runs. To put this in a more specific setting, the code below use this process to generate the path of a stock. And I want the gradient of the sample mean over several stock paths at time T =1.0, w.r.t. to the parameter a, for calibration/market data fitting purposes. However, the output from the derivative is NaN, and I have a hard time figuring out where it goes wrong. Maybe its not possible to get forwardiff to work for this problem?

using ForwardDiff
using Gadfly
using DualNumbers
using BenchmarkTools


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 sqrt(1.0/n);
end

function get_sigma_Z2(a,n)
  return sqrt(1.0/((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 direct_conv{T}(a::Array{T}, b::Array{T})
    m = length(a)
    n = length(b)
    c = zeros(T,m+n-1)
     @inbounds @simd for j=1:m
        @inbounds @simd for k=1:n
            c[j+k-1] += a[j]*b[k]
        end
    end
    return c
end

function gen_volterra_for(a,s,n,rand_1,rand_2)
  Z_1 = zeros(eltype(a),s); #first normal
  Z_2 = zeros(eltype(a),s+1); # 2nd normal
  G = zeros(eltype(a),s+1); #Kernel function
   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+1] = get_sigma_Z2(a,n)*(get_rho(a)*rand_1[i]+sqrt(1-get_rho(a)^2.0)*rand_2[i]); #set correct variance of random variate
  end
  return (direct_conv(Z_1,G)[1:s+1]+Z_2); #Convolute
  #return G
end

function gen_path_for(a,rho,eta,xi,S_0,s,n,dt, rand_nums)
  vol = zeros(eltype(a),s);
  dZ = zeros(eltype(a),s);
  stock = zeros(eltype(a),s+1);
  stock[1] = S_0;
  volterra = gen_volterra_for(a,s,n,rand_nums[1:n],rand_nums[n+1:2*n]);
  S_incr = 0.0;
  t = 0.0;
  for i=1:s
    vol[i] = xi*exp(eta*sqrt(2.0*a+1.0)*volterra[i]-0.5*eta^2*t^(2.0*a+1.0));
    dZ[i] = rho*get_sigma_Z1(n)*rand_nums[i]+ sqrt(1.0+rho^2.0)*rand_nums[2*n+i]*dt;
    S_incr = S_incr +sqrt(vol[i])*dZ[i]-0.5*vol[i]*dt;
    stock[i+1] = S_0*exp(S_incr);
    t += dt;
  end
  return stock;
end

function gen_paths(a,num_paths)
  rho = -0.9;
  eta = 1.9;
  S_0 = 1.0
  xi = 0.235^2;
  n = 312;
  s = round(Int,n*1.0;)
  dt = 1.0/n;
  mean = 0;
  for p = 1:num_paths
    rand_nums = randn(3*n);
    mean +=gen_path_for(a,rho,eta,xi,S_0,s,n,dt,rand_nums)[s];
  end
  mean =mean/num_paths;
  return mean;
end

@time m = ForwardDiff.derivative(a ->gen_paths(a,1000),-0.43)
14

Seems like this is the problematic piece:

In particular this part t^(2.0*a+1.0) when t = 0.0. I recreated the error using:

julia> f(x,t) = t^(2.0*x+1.0)
f (generic function with 1 method)

julia> ForwardDiff.derivative(x->f(x,0.),-0.43)
NaN

This seems like an error with ForwardDiff. To work around this I used the following:

    if t == 0.0
      vol[i] = xi*exp(eta*sqrt(2.0*a+1.0)*volterra[i]);
    else
      vol[i] = xi*exp(eta*sqrt(2.0*a+1.0)*volterra[i]-0.5*eta^2*t^(2.0*a+1.0));
    end

Here is the full code with a couple of other minor changes:

using ForwardDiff
using Gadfly
using DualNumbers
using BenchmarkTools

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 sqrt(1.0/n);
end

function get_sigma_Z2(a,n)
  return sqrt(1.0/((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 direct_conv{T}(a::Array{T}, b::Array{T})
    m = length(a)
    n = length(b)
    c = zeros(T,m+n-1)
    @inbounds @simd for j=1:m
        @inbounds @simd for k=1:n
            c[j+k-1] += a[j]*b[k]
        end
    end
    return c
end

function gen_volterra_for(a,s,n,rand_1,rand_2)
  Z_1 = zeros(eltype(a),s); #first normal
  Z_2 = zeros(eltype(a),s+1); # 2nd normal
  G = zeros(eltype(a),s+1); #Kernel function
  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+1] = get_sigma_Z2(a,n)*(get_rho(a)*rand_1[i]+sqrt(1-get_rho(a)^2.0)*rand_2[i]); #set correct variance of random variate
  end
  return direct_conv(Z_1,G)[1:s+1]+Z_2; #Convolute
  #return G
end

function gen_path_for(a,rho,eta,xi,S_0,s,n,dt, rand_nums)
  vol = zeros(eltype(a),s);
  dZ = zeros(eltype(a),s);
  stock = zeros(eltype(a),s+1);
  stock[1] = S_0;
  volterra = gen_volterra_for(a,s,n,rand_nums[1:n],rand_nums[n+1:2*n]);
  S_incr = zero(a); #Note this change
  t = 0.0;
  for i=1:s
    if t == 0.0
      vol[i] = xi*exp(eta*sqrt(2.0*a+1.0)*volterra[i]);
    else
      vol[i] = xi*exp(eta*sqrt(2.0*a+1.0)*volterra[i]-0.5*eta^2*t^(2.0*a+1.0));
    end
    dZ[i] = rho*get_sigma_Z1(n)*rand_nums[i]+ sqrt(1.0+rho^2.0)*rand_nums[2*n+i]*dt;
    S_incr = S_incr + sqrt(vol[i])*dZ[i]-0.5*vol[i]*dt;
    stock[i+1] = S_0*exp(S_incr);
    t += dt;
  end
  return stock;
end

function gen_paths(a,num_paths)
  rho = -0.9;
  eta = 1.9;
  S_0 = 1.0
  xi = 0.235^2;
  n = 312;
  s = n; #Note this change
  dt = 1.0/n;
  mean = zero(a); #Note this change
  for p = 1:num_paths
    rand_nums = randn(3*n);
    temp = gen_path_for(a,rho,eta,xi,S_0,s,n,dt,rand_nums)[s];
    mean += temp;
  end
  mean = mean / num_paths;
  return mean;
end

@time m = ForwardDiff.derivative(a -> gen_paths(a,1000),-0.43)

Also it might be a good idea to check that ForwardDiff correctly differentiates the rest of your code by possibly comparing it with a finite difference approximation.

15

Out of curiosity, since I know little of Ito calculus: computing the expected derivative using the kernel \alpha g(t-s)^{\alpha-1} gives the wrong result?

16

Well, now we are looking at \frac{\partial}{\partial \alpha}\mathbb{E}_0[S_T] approximated by \frac{\partial}{\partial \alpha}\frac{1}{N}\sum_i^NS_T^{(i)}. But, as usual we have that the stock price is a martingale so, \mathbb{E}_0[S_T]=S_0=1.0. So in my head i should get values of \frac{\partial}{\partial \alpha}\frac{1}{N}\sum_i^NS_T^{(i)} converging to zero as i increase N.

17

Okey, guys this now runs, thanks so much for the help so far! To take things a bit further, I have now implemented a code such that I can compute the gradient of the implied volatility generated by the stochastic model with respect to my parameters (a, \rho, \eta ). The estimate of the implied vol for a given strike K is loosely defined by \hat\sigma_{imp} = BlackScholes^{-1}(\mathbb{E}_0[S_T-K]), which computation is facilitated by the blsimpv(spot,K,r,T,Price,d) function in FinancialToolbox.jl package (Fortunately this supports dual numbers!! ).

So, the code below computes the gradient wrt. to the implied volatility for an at-the-money option (K=S_0) using forwarddiff and compares it to finite difference estimates(Which are the same :raised_hands:). However, I really need some help with speeding this code up, as I am planning to use this for calibration to market implied volatilities (Is forwarddiff compatible with multithreading? )

using ForwardDiff
using Gadfly
using DualNumbers
using BenchmarkTools
using FinancialToolbox

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 sqrt(1.0/n);
end

function get_sigma_Z2(a,n)
  return sqrt(1.0/((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 direct_conv{T}(a::Array{T}, b::Array{T})
    m = length(a)
    n = length(b)
    c = zeros(T,m+n-1)
    @inbounds @simd for j=1:m
        @inbounds @simd for k=1:n
            c[j+k-1] += a[j]*b[k]
        end
    end
    return c
end

function gen_volterra_for(a,s,n,rand_1,rand_2)
  Z_1 = zeros(eltype(a),s); #first normal
  Z_2 = zeros(eltype(a),s+1); # 2nd normal
  G = zeros(eltype(a),s+1); #Kernel function
  @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+1] = get_sigma_Z2(a,n)*(get_rho(a)*rand_1[i]+sqrt(1-get_rho(a)^2.0)*rand_2[i]); #set correct variance of random variate
  end
  return direct_conv(Z_1,G)[1:s+1]+Z_2; #Convolute
  #return G
end

function gen_path_for(a,rho,eta,xi,S_0,s,n,dt, rand_nums)
  vol = zeros(eltype(a),s);
  dZ = zeros(eltype(a),s);
  stock = zeros(eltype(a),s+1);
  stock[1] = S_0;
  volterra = gen_volterra_for(a,s,n,rand_nums[1:n],rand_nums[n+1:2*n]);
  S_incr = zero(a);
  t = zero(a);
  @inbounds @simd for i=1:s
    if t == 0.0
      vol[i] = xi;
    else
      vol[i] = xi*exp(eta*sqrt(2.0*a+1.0)*volterra[i]-0.5*eta^2*t^(2.0*a+1.0));
    end
    dZ[i] = rho*get_sigma_Z1(n)*rand_nums[i]+ sqrt(1.0+rho^2.0)*rand_nums[2*n+i]*dt;
    S_incr = S_incr + sqrt(vol[i])*dZ[i]-0.5*vol[i]*dt;
    stock[i+1] = S_0*exp(S_incr);
    t += dt;
  end
  return stock;
end

function gen_paths(x::Vector,num_paths)
  rng = MersenneTwister(1234);
  #rho = -0.9;
  #eta = 1.9;
  S_0 = 1.0;
  K =S_0;
  xi = 0.235^2;
  T=1.0;
  n = 312;
  s = s = round(Int,n*T;)
  dt = 1.0/n;
  mean = zero(a);
  @inbounds @simd for p = 1:num_paths
    rand_nums = randn(rng,3*n);
    temp = max(gen_path_for(x[1],x[2],x[3],xi,S_0,s,n,dt,rand_nums)[s]-K,0);
    mean += temp;
  end
  mean = mean / num_paths;
  imp_vol = blsimpv(S_0,K,0.0,T, mean)

  return imp_vol;
end

a = -0.40;
rho = -0.9;
eta = 1.9;

x = [a,rho,eta]
@time m = ForwardDiff.gradient(x -> gen_paths(x,10000),x)

##Comparing with finite difference estimates
@time imp_vol = gen_paths(x,10000)
delta =0.000001;
x_delta = [a+delta,rho,eta]
imp_vol_1 = gen_paths(x_delta, 10000)
fd_est_a = (imp_vol_1-imp_vol)/delta_alpha

x_delta = [a,rho+delta,eta]
imp_vol_1 = gen_paths(x_delta, 10000)
fd_est_a = (imp_vol_1-imp_vol)/delta_alpha

x_delta = [a,rho,eta+delta]
imp_vol_1 = gen_paths(x_delta, 10000)
fd_est_a = (imp_vol_1-imp_vol)/delta_alpha
18

I’m down to help, but where can I grab the FinancialToolbox package? It doesn’t seem to be registered.

19

FinancialToolbox here

Also see another topic, where we have investigated why the power operator ^ is slow in Julia: Why are calls to this function slow, compared to Python+Numpy?

1 Like
20

@jrevels did you have any hints for me, on how to speed this up? thx