Getting started with directConv functions

Table of Contents

1 Comparison with Julia native conv() function

With Julia native conv function

α=rand(5);
β=rand(10);
γ_ftt=conv(α,β)
14-element Array{Float64,1}:
 0.216624 
 0.416602 
 0.646538 
 0.660492 
 1.61936  
 1.15249  
 1.22127  
 1.72914  
 1.44925  
 1.29137  
 0.94214  
 0.65695  
 0.262058 
 0.0200272

conv function returns an array of size \(|\alpha|+|\beta|-1\), here

14

With DirectConvolution.jl one must first define a filter with offset=0, then use the directConv function:

α_filter=LinearFilter(α,0)
γ_direct=directConv(α_filter,β)
DirectConvolution.LinearFilter_Default{Float64,5}([0.0570966, 0.159702, 0.932971, 0.746808, 0.675408], 0)
10-element Array{Float64,1}:
 0.00823509
 0.0778808 
 0.322487  
 1.11396   
 1.45592   
 1.46361   
 1.37137   
 1.36605   
 1.82005   
 1.46266

A difference with the native Julia conv function is that the returned signal as the same size as the initial \(\beta\) signal.

Let's use first check that the results are identical (for their common parts):

norm(γ_ftt[1:endof(γ_direct)]-γ_direct)
4.718447854656915e-16

2 Performance

The following is not a complete benchmark, but gives some ideas. Please also note that for the moment directConv is a pure Julia, straightforward implementation.

With the given α, β and γ we get:

r_ftt=@benchmark(conv(α,β))
BenchmarkTools.Trial: 
  memory estimate:  9.17 KiB
  allocs estimate:  132
  --------------
  minimum time:     36.946 μs (0.00% GC)
  median time:      43.541 μs (0.00% GC)
  mean time:        51.735 μs (8.59% GC)
  maximum time:     21.249 ms (63.96% GC)
  --------------
  samples:          10000
  evals/sample:     1
r_direct=@benchmark(directConv(α_filter,β))
BenchmarkTools.Trial: 
  memory estimate:  432 bytes
  allocs estimate:  5
  --------------
  minimum time:     551.717 ns (0.00% GC)
  median time:      573.893 ns (0.00% GC)
  mean time:        694.201 ns (13.76% GC)
  maximum time:     34.393 μs (96.89% GC)
  --------------
  samples:          10000
  evals/sample:     187

For these short filters we have a

median(r_ftt.times)/median(r_direct.times)
75.86953726308728

speed up.

For a more realistic size for \(\beta\), we still get this speedup:

β=rand(1000);
r_ftt=@benchmark(conv(α,β));
r_direct=@benchmark(directConv(α_filter,β));
median(r_ftt.times)/median(r_direct.times)
24.075015278060704

Typical size for \(\alpha\) is around \(5\) for applications like wavelets transform. If we increase this size the direct approach is less and less competitive compared to the FFT approach. However on my computer it is still advantageous to use it for an \(\alpha\) size of \(40\).

α=rand(40);
α_filter=LinearFilter(α,0);
r_ftt=@benchmark(conv(α,β));
r_direct=@benchmark(directConv(α_filter,β));
median(r_ftt.times)/median(r_direct.times)
2.5229993396169776

3 Versatility

3.1 Parallel computation

The directConv functions can be used in various ways. For instance we can compute a convolution by splitting destination

α=LinearFilter(rand(10),0);
β=rand(100);

γ_ref = directConv(α,β);        # compute all in one call

γ = Vector{Float64}(length(β));
directConv!(α,-1,β,γ,1:25);     # independant computations 1:25
directConv!(α,-1,β,γ,25:75);    #                         25:75      
directConv!(α,-1,β,γ,75:100);   #                         75:100 

norm(γ_ref-γ) # check that results are identical
0.0

This can be used for parallel computations. (no done yet)

3.2 Boundary extensions

Boundary extensions are rigorously implemented and are free of all limitations. For instance you can even use a filter α bigger than the signal β, with mirror symmetry on left and periodic extension of the right:

α=LinearFilter([0,0,0,1,0,0,0],-6); # with offset=-6, mimic a filter of support 6:12
β=[1,2,3,4];
γ = zeros(Int,20);
directConv!(α,-1,β,γ,1:20,MirrorBE,PeriodicBE); 
plot(1:20, γ, markershape = :hexagon)

be.png

With zero padding we can also check commutativity \(\alpha \circledast \beta = \beta \circledast \alpha\):

α=rand(4);
β=rand(6);
tαβ=zeros(6);
directConv!(LinearFilter(α,0),-1,β,tαβ,1:6);
tβα=zeros(6);
directConv!(LinearFilter(β,0),-1,α,tβα,1:6);
norm(tαβ-tβα)
0.0

3.3 Adjoint

From \[ \gamma[k]=\sum\limits_{i\in\Omega^\alpha}\alpha[i]\beta[k+\lambda i],\text{ with }\lambda\in\mathbb{Z}^* \] with free \(\lambda\) we can compute convolution (\(\lambda=-1\)) cross-correlation (\(\lambda=+1\)) and more generally adjoint operator: \[ \forall \lambda\in\mathbb{Z}^*,\ \langle \alpha\bigodot\limits_{\lambda}v ,w \rangle_E = \langle v , \alpha\bigodot\limits_{-\lambda} w \rangle_F \]

This relation only holds true for zero padding and periodic boundary extension:

α=LinearFilter(rand(5));
v=rand(8);
w=rand(8);

dot(directConv(α,-3,v),w) - dot(directConv(α,+3,w),v)
dot(directConv(α,-3,v,PeriodicBE,PeriodicBE),w) - dot(directConv(α,+3,w,PeriodicBE,PeriodicBE),v)
0.0
0.0

Author: picaud

Created: 2018-06-25 Mon 18:58

Validate