Turing: problem using MCMCthreads()

Specific Domains Statistics
1

I am getting an error when trying to use MCMCthreads() for sampling. The following code works fine:

chain = sample(MSM(m, S, model),
    MH(:θt => AdvancedMH.RandomWalkProposal(MvNormal(zeros(size(m,1)), tuning*Σp))),
    length; init_params=m, discard_initial=burnin)

However, the following code

chain = sample(MSM(m, S, model),
    MH(:θt => AdvancedMH.RandomWalkProposal(MvNormal(zeros(size(m,1)), tuning*Σp))),
    MCMCThreads(), length, nchains;  init_params=m, discard_initial=burnin)

gives an error “Provided initial value doesn’t match the dimension of the model”.

Before either of these are called, nchains, length, burnin and m are all defined, and are the same for both of the above calls. Also, Julia is started with julia -t 4, and the environment variable JULIA_NUM_THREADS=4 is set.

The code is from the file https://github.com/mcreel/SimulatedNeuralMoments.jl/blob/main/examples/MN/MNexample.jl

Strangely enough, the file https://github.com/mcreel/SimulatedNeuralMoments.jl/blob/main/examples/SV/SVexample.jl

contains identical code:

chain = sample(MSM(m, S, model),
    MH(:θt => AdvancedMH.RandomWalkProposal(MvNormal(zeros(size(m,1)), tuning*Σp))),
    MCMCThreads(), length, nchains; init_params=m, discard_initial=burnin)

# single thread
#=
chain = sample(MSM(m, S, model),
    MH(:θt => AdvancedMH.RandomWalkProposal(MvNormal(zeros(3), tuning*Σp))),
    length; init_params = m, discard_initial=burnin)
=#

and both of these run without problems. If anyone sees the problem, that would be much appreciated. If not, I guess I should file an issue with Turing.

2

I’m not quite sure what is wrong. Can you produce the behavior with a simpler example, such as estimating mu and sigma of a normal distribution? If so, it might be worth filing an issue.

3

The odd thing is that the two files, MNexample.jl and SVexample.jl, are essentially identical from Turing’s point of view - they just implement two different statistical models. I will try to boil down the file that has the problem to a minimum working example.

1 Like
4

Here is a MWE, or better said, M non-working E.

With the 5 or 4 dimensional problems, the error is

 
        nested task error: Provided initial value doesn't match the dimension of the model

With the 3 dimensional problem, the error is

ERROR: LoadError: ArgumentError: not enough initial parameters (expected 4, received 3

In all cases, julia was started with julia -t 4 --proj where the project has the needed packages. If the block using MCMCThreads() is commented out, the single thread version works in all cases.

using Turing, AdvancedMH, LinearAlgebra


#lb = [0.0, 0.0, 0.0, 0.0, 0.05]
#ub = [3.0, 3.0, 1.0, 3.0, 0.95] 
#m = [1.0, 1.0, 0.2, 1.8, 0.4] 

lb = [0.0, 0.0, 0.0, 0.0]
ub = [3.0, 3.0, 1.0, 3.0] 
m = [1.0, 1.0, 0.2, 1.8] 

#lb = [0.0, 0.0, 0.0]
#ub = [3.0, 3.0, 1.0] 
#m = [1.0, 1.0, 0.2] 

k = size(m,1)
V = 1.0* Matrix(I, k, k)

macro Prior()
    return :( arraydist([Uniform(lb[i], ub[i]) for i = 1:size(lb,1)]) )
end

transf = bijector(@Prior) # transforms draws from prior to draws from  ℛⁿ 
transformed_prior = transformed(@Prior, transf) # the transformed prior

@model function MSM(m)
    θt ~ transformed_prior
    m ~ MvNormal(θt, V)
end

length = 1000
nchains = 4
# This does not work
chain = sample(MSM(m),
        MH(:θt => AdvancedMH.RandomWalkProposal(MvNormal(zeros(k), V))),
        MCMCThreads(), length, nchains;
        init_params=m)
#
# single thread, this works
chain = sample(MSM(m),
        MH(:θt => AdvancedMH.RandomWalkProposal(MvNormal(zeros(k), V))),
        length;
        init_params=m)
5

Update: the error is occurring with Turing v0.21.1, but it does not occur with v0.20.2. I will try to narrow down exactly which version provokes the error, and then I suppose I should file an issue with Turing.jl

UPDATE: it works with Turing 0.20.4, but not with 0.21.0.

in both cases, the other Turing related packages are at
AdvancedMH 0.6.7
MCMCChains 5.1.0

1 Like
6

I tend to agree. It sounds like there might be a bug, or an intended change in behavior. However, I did not see anything in the closed issues or merged pull requests that immediately suggest the later.

1 Like
7

I have filed an issue, at https://github.com/TuringLang/Turing.jl/issues/1810

The same problem occurs with NUTS, which is what is used in the example code in the issue.

1 Like
8

From the AbstractMCMC docs:

  • init_params (default: nothing): if set to init_params !== nothing, then the ith element of init_params is used as initial parameters of the ith chain. If one wants to use the same initial parameters x for every chain, one can specify e.g. init_params = Iterators.repeated(x) or init_params = FillArrays.Fill(x, N).

So when sampling a single chain, you pass a single set of initial parameters. When sampling multiple chains, you pass an iterator of initial parameters. e.g.

julia> using Turing

julia> @model function mod()
           x ~ filldist(Normal(), 10)
       end;

julia> x0 = randn(10);

julia> sample(mod(), NUTS(), 1_000; init_params=x0);
┌ Info: Found initial step size
└   ϵ = 1.6
Sampling 100%|█████████████████████████████████████████████████████████████████████████| Time: 0:00:00

julia> sample(mod(), NUTS(), MCMCThreads(), 1_000, 4; init_params=Iterators.repeated(x0));
┌ Warning: Only a single thread available: MCMC chains are not sampled in parallel
└ @ AbstractMCMC ~/.julia/packages/AbstractMCMC/0eT8o/src/sample.jl:291
┌ Info: Found initial step size
└   ϵ = 1.6
┌ Info: Found initial step size
└   ϵ = 1.6
┌ Info: Found initial step size
└   ϵ = 1.6
┌ Info: Found initial step size
└   ϵ = 1.6
Sampling (1 threads) 100%|█████████████████████████████████████████████████████████████| Time: 0:00:00
3 Likes
9

Thanks, that solved it. I will close the issue I opened.