Performance Tips

This section briefly summarises a few common techniques to ensure good performance when using Turing. We refer to julialang.org for general techniques to ensure good performance of Julia programs.

Use multivariate distributions

It is generally preferable to use multivariate distributions if possible.

The following example:

@model gmodel(x) = begin
m ~ Normal()
for i = 1:length(x)
x[i] ~ Normal(m, 0.2)
end
end


can be directly expressed more efficiently using a simple transformation:

@model gmodel(x) = begin
m ~ Normal()
x ~ MvNormal(fill(m, length(x)), 0.2)
end


Turing currently provides support for two different automatic differentiation (AD) backends. Generally, try to use :forwarddiff for models with few parameters and :reversediff, :tracker or :zygote for models with large parameter vectors or linear algebra operations. See Automatic Differentiation for details.

Special care for :tracker and :zygote

In case of :tracker and :zygote, it is necessary to avoid loops for now. This is mainly due to the reverse-mode AD backends Tracker and Zygote which are inefficient for such cases. ReverseDiff does better but vectorized operations will still perform better.

Avoiding loops can be done using filldist(dist, N) and arraydist(dists). filldist(dist, N) creates a multivariate distribution that is composed of N identical and independent copies of the univariate distribution dist if dist is univariate, or it creates a matrix-variate distribution composed of N identical and idependent copies of the multivariate distribution dist if dist is multivariate. filldist(dist, N, M) can also be used to create a matrix-variate distribution from a univariate distribution dist. arraydist(dists) is similar to filldist but it takes an array of distributions dists as input. Writing a custom distribution with a custom adjoint is another option to avoid loops.

For efficient gradient-based inference, e.g. using HMC, NUTS or ADVI, it is important to ensure the model is type-stable. We refer to julialang.org for a general discussion on type-stability.

The following example:

@model tmodel(x, y) = begin
p,n = size(x)
params = Vector{Real}(undef, n)
for i = 1:n
params[i] ~ truncated(Normal(), 0, Inf)
end

a = x * params
y ~ MvNormal(a, 1.0)
end


can be transformed into the following type-stable representation:

@model tmodel(x, y, ::Type{T}=Vector{Float64}) where {T} = begin
p,n = size(x)
params = T(undef, n)
for i = 1:n
params[i] ~ truncated(Normal(), 0, Inf)
end

a = x * params
y ~ MvNormal(a, 1.0)
end


Note that you can use @code_warntype to find type instabilities in your model definition.

For example consider the following simple program

@model tmodel(x) = begin
p = Vector{Real}(undef, 1);
p[1] ~ Normal()
p = p .+ 1
x ~ Normal(p[1])
end


we can use

using Random

model = tmodel(1.0)
varinfo = Turing.VarInfo(model)
spl = Turing.SampleFromPrior()

@code_warntype model.f(Random.GLOBAL_RNG, model, varinfo, spl, Turing.DefaultContext())


to inspect the type instabilities in the model.

Reuse Computations in Gibbs Sampling

Often when performing Gibbs sampling, one can save computational time by caching the output of expensive functions. The cached values can then be reused in future Gibbs sub-iterations which do not change the inputs to this expensive function. For example in the following model:

@model demo(x) = begin
a ~ Gamma()
b ~ Normal()
c = function1(a)
d = function2(b)
x .~ Normal(c, d)
end
alg = Gibbs(MH(:a), MH(:b))
sample(demo(zeros(10)), alg, 1000)


when only updating a in a Gibbs sub-iteration, keeping b the same, the value of d doesn’t change. And when only updating b, the value of c doesn’t change. However, if function1 and function2 are expensive and are both run in every Gibbs sub-iteration, a lot of time would be spent computing values that we already computed before. Such a problem can be overcome using Memoization.jl. Memoizing a function lets us store and reuse the output of the function for every input it is called with. This has a slight time overhead but for expensive functions, the savings will be far greater.

To use Memoization.jl, simply define memoized versions of function1 and function2 as such:

using Memoization

@memoize memoized_function1(args...) = function1(args...)
@memoize memoized_function2(args...) = function2(args...)


Then define the Turing model using the new functions as such:

@model demo(x) = begin
a ~ Gamma()
b ~ Normal()
c = memoized_function1(a)
d = memoized_function2(b)
x .~ Normal(c, d)
end