Kristoffer, your example still uses a *non-const* global. Using this benchmark <https://gist.github.com/tlycken/f0dc12b7c211f73f6cc1>, I get the following results:
Non-const global 0.017629 seconds (999.84 k allocations: 15.256 MB) Const global 0.000002 seconds (6 allocations: 192 bytes) Closure 0.000008 seconds (19 allocations: 928 bytes) So yes, a closure is much better than using a non-const global, but using a const global still beats it. However, both are fast enough that it probably won’t be the bottleneck in your application anymore :) Although I can by no means read and understand it, it is also instructive to look at the dumps from @code_llvm to get a picture of the difference in code complexity: Const global: julia> @code_llvm f_const_glob(10^5) define i64 @julia_f_const_glob_2617(i64) { top: %1 = icmp slt i64 %0, 1 br i1 %1, label %L3, label %L.preheader L.preheader: ; preds = %top %2 = icmp sgt i64 %0, 0 %.op = mul i64 %0, 3 %3 = select i1 %2, i64 %.op, i64 0 br label %L3 L3: ; preds = %L.preheader, %top %s.1 = phi i64 [ 0, %top ], [ %3, %L.preheader ] ret i64 %s.1 } Closure: julia> @code_llvm f_clos(10^5, 3) define %jl_value_t* @julia_f_clos_2616(i64, i64) { top: %2 = alloca [5 x %jl_value_t*], align 8 %.sub = getelementptr inbounds [5 x %jl_value_t*]* %2, i64 0, i64 0 %3 = getelementptr [5 x %jl_value_t*]* %2, i64 0, i64 2 %4 = getelementptr [5 x %jl_value_t*]* %2, i64 0, i64 3 store %jl_value_t* inttoptr (i64 6 to %jl_value_t*), %jl_value_t** %.sub, align 8 %5 = getelementptr [5 x %jl_value_t*]* %2, i64 0, i64 1 %6 = load %jl_value_t*** @jl_pgcstack, align 8 %.c = bitcast %jl_value_t** %6 to %jl_value_t* store %jl_value_t* %.c, %jl_value_t** %5, align 8 store %jl_value_t** %.sub, %jl_value_t*** @jl_pgcstack, align 8 store %jl_value_t* null, %jl_value_t** %3, align 8 store %jl_value_t* null, %jl_value_t** %4, align 8 %7 = getelementptr [5 x %jl_value_t*]* %2, i64 0, i64 4 store %jl_value_t* null, %jl_value_t** %7, align 8 %8 = load %jl_value_t** inttoptr (i64 2147534600 to %jl_value_t**), align 8 store %jl_value_t* %8, %jl_value_t** %4, align 8 %9 = load %jl_value_t** inttoptr (i64 2164502072 to %jl_value_t**), align 8 store %jl_value_t* %9, %jl_value_t** %7, align 8 %10 = call %jl_value_t* @jl_f_instantiate_type(%jl_value_t* null, %jl_value_t** %4, i32 2) store %jl_value_t* %10, %jl_value_t** %4, align 8 %11 = call %jl_value_t* @jl_f_svec(%jl_value_t* null, %jl_value_t** null, i32 0) store %jl_value_t* %11, %jl_value_t** %7, align 8 %12 = call %jl_value_t* @jl_f_svec(%jl_value_t* null, %jl_value_t** %4, i32 2) store %jl_value_t* %12, %jl_value_t** %4, align 8 %13 = call %jl_value_t* @jl_box_int64(i64 signext %1) store %jl_value_t* %13, %jl_value_t** %7, align 8 %14 = call %jl_value_t* (i64, ...)* @jl_svec(i64 1, %jl_value_t* %13) store %jl_value_t* %14, %jl_value_t** %7, align 8 %15 = call %jl_value_t* @jl_new_closure(i8* null, %jl_value_t* %14, %jl_value_t* inttoptr (i64 2221717744 to %jl_value _t*)) store %jl_value_t* %15, %jl_value_t** %7, align 8 %16 = load %jl_value_t** @jl_false, align 8 %17 = call %jl_value_t* @jl_method_def(%jl_value_t* inttoptr (i64 507171808 to %jl_value_t*), %jl_value_t** %3, %jl_va lue_t* null, %jl_value_t* null, %jl_value_t* %12, %jl_value_t* %15, %jl_value_t* %16, %jl_value_t* inttoptr (i64 2179532 144 to %jl_value_t*), i32 0) %18 = load %jl_value_t** %3, align 8 %19 = bitcast %jl_value_t* %18 to %jl_value_t* (%jl_value_t*, %jl_value_t**, i32)** %20 = load %jl_value_t* (%jl_value_t*, %jl_value_t**, i32)** %19, align 8 %21 = call %jl_value_t* @jl_box_int64(i64 signext %0) store %jl_value_t* %21, %jl_value_t** %4, align 8 %22 = call %jl_value_t* %20(%jl_value_t* %18, %jl_value_t** %4, i32 1) %23 = load %jl_value_t** %5, align 8 %24 = getelementptr inbounds %jl_value_t* %23, i64 0, i32 0 store %jl_value_t** %24, %jl_value_t*** @jl_pgcstack, align 8 ret %jl_value_t* %22 } On Tuesday, September 29, 2015 at 9:18:29 AM UTC+2, Kristoffer Carlsson wrote: In my experience a closure is much faster than a function that accesses > global variables even if it is passed to another function as an argument. > > Two different examples. Here we pass a function that accesses (non const) > globals to another function: > > a = 3 > > g(f::Function, b, N) = f(b, N) > > function f_glob(b, N) > s = 0 > for i = 1:N > s += a > end > return s > end > > @time g(f_glob, 3, 10^6) > > 3000000 > > 0.020655 seconds (999.84 k allocations: 15.256 MB, 18.40% gc time) > > > Here we instead use an outer function that generates the desired closure > and pass the closure as the function argument: > > function f_clos(b, N, a) > function f_inner(b, N) > s = 0 > for i = 1:N > s += a > end > return s > end > > g(f_inner) > end > > @time f_clos(3, 10^6, 2) > > 2000000 > > 0.000021 seconds (19 allocations: 928 bytes) > > > > > On Tuesday, September 29, 2015 at 8:47:06 AM UTC+2, Tomas Lycken wrote: >> >> The performance penalty from global variables comes from the fact that >> non-const global variables are type unstable, which means that the compiler >> has to generate very defensive, and thus slow, code. The same thing is, for >> now, true also for anonymous functions and functions passed as parameters >> to other functions, so switching one paradigm for the other will >> unfortunately not help much for performance. I haven't looked closely at >> your code, but the [FastAnonymous.jl]( >> https://github.com/timholy/FastAnonymous.jl) package can probably help >> with making the anonymous function version faster. >> >> If the data you want to pass along doesn't change, using *const* globals >> is fine (at least from a performance perspective, you might have other >> reasons to avoid them too...), i.e. you can do the following and it will >> also be faster: >> >> ``` >> const data = 3 >> >> function model(x) >> # use data somehow >> end >> >> # use the model function in your optimization >> ``` >> >> If you need the data to change, you can still use a `const` global >> variable, but set it to some mutable (and type-stable!) object, for example >> >> ``` >> type MyData >> data::Int >> end >> >> const data = MyData(3) >> >> function model(x) >> data.data += 1 >> # ... >> end >> >> # use the model function in your optimization >> # it can now modify the global data on each invocation >> ``` >> >> This doesn't necessarily provide the cleanest interface, but it should be >> more performant than your current solution. >> >> // T >> >> >> On Monday, September 28, 2015 at 6:53:22 PM UTC+2, Christopher Fisher >> wrote: >>> >>> Thanks for your willingness to investigate the matter more closely. I >>> cannot post the exact code I am using (besides its rather long). However, I >>> posted a toy example that follows the same basic operations. Essentially, >>> my code involves an outer function (SubLoop) that loops through a data set >>> with multiple subjects. The model is fit to each subject's data. The other >>> function (LogLike) computes the log likelihood and is called by optimize. >>> The first set of code corresponds to the closure method and the second set >>> of code corresponds to the global variable method. In both cases, the code >>> executed in about .85 seconds over several runs on my computer and has >>> about 1.9% gc time. I'm also aware that my code is probably not optimized >>> in other regards. So I would be receptive to any other advice you might >>> have. >>> >>> >>> >>> >>> using Distributions,Optim >>> >>> function SubLoop1(data1) >>> >>> function LogLike1(parms) >>> >>> L = pdf(Normal(parms[1],exp(parms[2])),SubData) >>> >>> LL = -sum(log(L)) >>> >>> end >>> >>> #Number of Subjects >>> >>> Nsub = size(unique(data1[:,1],1),1) >>> >>> #Initialize per subject Data >>> >>> SubData = [] >>> >>> for i = 1:Nsub >>> >>> idx = data1[:,1] .== i >>> >>> SubData = data1[idx,2] >>> >>> parms0 = [1.0;1.0] >>> >>> optimize(LogLike1,parms0,method=:nelder_mead) >>> >>> end >>> >>> end >>> >>> >>> >>> N = 10^5 >>> >>> #Column 1 subject index, column 2 value >>> >>> Data = zeros(N*2,2) >>> >>> for sub = 1:2 >>> >>> Data[(N*(sub-1)+1):(N*sub),:] = [sub*ones(N) rand(Normal(10,2),N)] >>> >>> end >>> >>> @time SubLoop1(Data) >>> >>> >>> >>> >>> >>> >>> >>> >>> >>> >>> >>> Using Distributions, Optim >>> >>> function SubLoop2(data1) >>> >>> global SubData >>> >>> #Number of subjects >>> >>> Nsub = size(unique(data1[:,1],1),1) >>> >>> #Initialize per subject data >>> >>> SubData = [] >>> >>> for i = 1:Nsub >>> >>> idx = data1[:,1] .== i >>> >>> SubData = data1[idx,2] >>> >>> parms0 = [1.0;1.0] >>> >>> optimize(LogLike2,parms0,method=:nelder_mead) >>> >>> end >>> >>> end >>> >>> >>> >>> function LogLike2(parms) >>> >>> L = pdf(Normal(parms[1],exp(parms[2])),SubData) >>> >>> LL = -sum(log(L)) >>> >>> end >>> >>> >>> >>> N = 10^5 >>> >>> #Column 1 subject index, column 2 value >>> >>> Data = zeros(N*2,2) >>> >>> for sub = 1:2 >>> >>> Data[(N*(sub-1)+1):(N*sub),:] = [sub*ones(N) rand(Normal(10,2),N)] >>> >>> end >>> >>> @time SubLoop2(Data) >>> >>> >>> >>> On Monday, September 28, 2015 at 11:24:13 AM UTC-4, Kristoffer Carlsson >>> wrote: >>>> >>>> From only that comment alone it is hard to give any further advice. >>>> >>>> What overhead are you seeing? >>>> >>>> Posting runnable code is the best way to get help. >>>> >>>>