I had a model that was running until about two weeks ago, but suddenly stopped. I thought I’d post here before filing an issue - it’s not clear to me if this is more of a reversediff issue or a Turing issue or something else.
Here’s a MWE:
using Turing, BenchmarkTools, ReverseDiff
@model function efa(P)
λ ~ filldist(Normal(0, 10), P, 1)
Ψ ~ filldist(truncated(Normal(0, 5), 0, Inf), P)
f_var ~ filldist(truncated(Normal(0, 5), 0, Inf), 1)
Σ1 = λ * f_var * λ' + Diagonal(Ψ) # Fails
# Σ1 = λ * f_var * λ' + diagm(Ψ) # Fails
# Σ1 = randn(P, P) + Diagonal(Ψ) # Works
# Σ1 = λ * f_var * λ' + randn(P, P) # Works
# Σ1 = λ * f_var * λ' + Diagonal(ones(P)) # Works
end
Turing.setadbackend(:reversediff)
efa_test = efa(3)
@time m1 = sample(efa_test, NUTS(1, 0.65), MCMCThreads(), 1, 4)
When I run the model above, I get an error in reversediff (in ReverseDiff’s increment_deriv! function, specifically). I narrowed it down to the line with Σ1. The first term on that line, λ * f_var * λ', creates a PxP matrix, to which I want to add the second term, a P-variate Diagonal matrix. Trying to do so leads to the error.
What makes things weirder is that when I replace either of the two terms with a random PxP matrix, the model runs fine. Also, if I replace the Ψ vector with ones(P), it runs fine. It’s not an issue with the first term, and it’s not the second term, but seems to be something specifically involving the combination of the two.
Am I missing something obvious?