ForwardDiff type error?

[Red-Portal]

in _unsafe_setindex!, in typeassert, expected Float64, got ForwardDiff.Dual{Nothing,Float64,10}

I'm getting this error while using HMC.

    @turing.model LDA(x, K, α, β) = begin
        N, D = size(x)
        L = 5 # ordinal levels

        ϕ = Array{Float64, 3}(undef, L, D, K);
        for i in 1:D
            for j in 1:K
                ϕ[:, i, j] ~ Dirichlet(L, α)
            end
        end

        θ = Matrix{Float64}(undef, K, N)
        for i in 1:N
            θ[:,i] ~ Dirichlet(K, β)
        end

        marginals = reshape(reshape(ϕ, (L * D, K)) * θ,
                            (L, D, N))
        for i in 1:N
            for j in 1:D
                x[i, j] ~ CategoricalReal(marginals[:,j,i])
            end
        end
        return ϕ, θ
    end

This is the model I'm currently using.
The error is caused here.

ϕ[:, i, j] ~ Dirichlet(L, α)

Also, I experienced problems while using the vectorized notation.

ϕ = Array{Float64, 3}(undef, L * D, K);
ϕ ~ [Dirichlet(K, alpha)]

This causes a dimension error.
How can I correctly used the vectorized notation?

[cpfiffer]

The issue here comes from the fact that the values inside θ = Matrix{Float64}(undef, K, N) and ϕ = Array{Float64, 3}(undef, L, D, K) are constrained to be of type Float64. The automatic differentiation converts those values from Float64 to ForwardDiff.Dual{Nothing,Float64,10}.

You can fix this by changing Float64 to Real in your array/matrix definition; use θ = Matrix{Real}(undef, K, N) and ϕ = Array{Real, 3}(undef, L, D, K) instead.

As to the list comprehension, I'm not sure how exactly it works with matrices like that. You might benefit from changing Array{Float64, 3} to Vector{Vector{Real}} sinceDirichlet is vector-valued.

Someone else might be more helpful (@trappmartin?) on the list comprehension side. I know there's been some discussion about improving it in the past (#604).

[Red-Portal]

@cpfiffer Thanks that fixed the issue. However after changing to Real instead of Float64 I noticed serious performance degradation. Seems like the matrix-matrix product below is not performed in parallel.

reshape(ϕ, (L * D, K)) * θ

I do not have deep enough insight to judge this issue but seems like OpenBLAS gemm is not called when Real is used

[trappmartin]

You are correct, using Real will cause performance issues. I would not recommend this if not necessary.

You might be able to get away with something like this:

θ = Vector{Vector}(undef, N)
for i in 1:N
    θ[i] ~ Dirichlet(K, β)
end

[Red-Portal]

That answers my question so I'll close this. Hope documentation on performance tips can be imrpoved

[v-a-s-a]

Hi @Red-Portal, I'm running into the same issue. Were you able to find a suitable workaround?

[cpfiffer]

The performant fix is in the performance tips section of the docs. Does that help?

[v-a-s-a]

That does indeed solve it, thank you @cpfiffer!