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?