search

gaurav-arya / StochasticAD.jl

Research package for automatic differentiation of programs containing discrete randomness.
MIT License
198 stars 16 forks source link

Infinitesimal derivative part of `StochasticTriples` end up as `NaN` for some instances of `MVNormal` #79

Closed frankschae closed 1 year ago

frankschae commented 1 year ago

MWE:

using StochasticAD, Distributions, LinearAlgebra, ForwardDiff
ϕ = 0.3
function mysample_q(ϕ)
    q_αβ = MvNormal(zeros(4), ϕ * Diagonal(rand(4)))
    return pdf(q_αβ, rand(4))
end

mysample_q(ϕ)
stochastic_triple(x -> mysample_q(x), ϕ) # 0.019046573788085913 + NaNε
ForwardDiff.derivative(x -> mysample_q(x), ϕ) # works

Note that MvNormal(zeros(4), ϕ * Array(Diagonal(rand(4)))) works fine.

gaurav-arya commented 1 year ago

This minimizes to the following issue:

julia> using StochasticAD

julia> st = stochastic_triple(0.5)
StochasticTriple of Float64:
0.5 + 1.0ε

julia> zero(1 / zero(st)) # should be 0.0 + 0.0ε
StochasticTriple of Float64:
0.0 + NaNε

This is triggered by https://github.com/JuliaStats/PDMats.jl/blob/fff131e11e23403931a42f5bfb3384f0d2b114c9/src/utils.jl#L63.

With #80 we should get this right.