See "method lookup hosting" and "function-valued argument inlining" in https://github.com/JuliaLang/julia/issues/3440. One good workaround is to pass a type instead of a function. If you define:
type X; end evaluate(::X) = ... then you can pass X() to a function and evaluate(t) and will be as fast as if you defined a function in the global scope. In this particular case, you might consider just using the kde function from StatsBase, which uses an FFT-based algorithm that should be much faster. Simon On Wednesday, March 5, 2014 8:57:09 PM UTC-5, Antonio Linero wrote: > > Below is a (highly non-optimal) implementation of kernel density > estimation. In principle, I'd like to pass in different kernel functions > (e.g. use the same code with a Gaussian kernel or triweight kernel). > However, I get a huge performance hit if I try to pass the kernel function > in. The following code calculates, for the Gaussian kernel, the KDE on a > grid of size 1000 given 1000 training data points. The first implementation > gets the kernel function from the global environment while the second takes > the kernel function passed in as an argument. > > Why does this give different performance and is there any way to repair > it? The version that does not pass the kernel function is (1) at least > three times faster and (2) allocates far more memory. > > using Rmath > > > function kern(x::Float64, y::Float64, h::Float64) > return Rmath.dnorm(x, y, h) > end > > > function eval_kern(x::Array{Float64, 1}, X::Array{Float64, 1}, h::Float64) > N = length(x) > K = length(X) > out = zeros(N) > for i in 1:N > for k in 1:K > out[i] += kern(x[i], X[k], h) > end > end > return out / N > end > > > function eval_kern(x::Array{Float64, 1}, X::Array{Float64, 1}, h::Float64, > Kern::Function) > N = length(x) > K = length(X) > out = zeros(N) > for i in 1:N > for k in 1:K > out[i] += Kern(x[i], X[k], h) > end > end > return out / N > end > > N = 1000 > X = randn(N) > h = std(X) * N ^ (-.2) * 1.06 > x_gridd = linspace(-4.0, 4.0, 1000); > > @time y_gridd = eval_kern(x_gridd, X, h); > @time y_gridd = eval_kern(x_gridd, X, h, kern); > > This gives the following results for me: > > julia> @time y_gridd = eval_kern(x_gridd, X, h); > elapsed time: 0.062120766 seconds (422240 bytes allocated) > > julia> @time y_gridd = eval_kern(x_gridd, X, h, kern); > elapsed time: 0.249457403 seconds (104151016 bytes allocated) > > > > I've experimented with a few things. Using Rmath.dnorm directly rather > than kern cuts down on the memory allocation for the first implementation > but doesn't seem to change the run time.. If I use > out[i] += Kern(x[i]::Float64, X[k]::Float64, h::Float64)::Float64 > instead I get within a factor of 2 in terms of speed and a 50% decreasing > in memory allocation. >