Third-order derivatives of `x / y` are `Inf` or `NaN` for small `y`, e.g. in gradients of Laplace approximations

Chiusa
  • Qualcosa non va?
#866 0 commenti 0 reazioni 0 assegnatari Vedi su GitHub

I maintainer di solito rispondono entro 1 giorno

Una pull request collegata è già stata integrata.

  • #867 di @devmotion — integrata

Valutazione

Difficoltà
3/5
Tempo stimato
1-2 giorni
Idoneità per principianti
65/100
Tipo di issue
Bug
Chiarezza
Specificata chiaramente
Stato di attività
Attiva
Stack tecnologico
julia
Ambito
backend, testing

Direzione di ricerca

Start by reading the division partials in src/partials.jl and the dual-number division handling in src/dual.jl, then reproduce the reported third-derivative cases with ForwardDiff. The issue proposes computing partials as (a - (aval / bval) * b) / bval; check the results against the finite derivatives shown in the report. Done means small-denominator derivatives remain finite and accurate in these examples.

Scritto dal modello di indicizzazione a partire dal testo della issue.

Descrizione

Gradients of Laplace approximations of marginal likelihoods require third derivatives of the log-density: the Laplace approximation contains the log-determinant of the Hessian with respect to the latent variables, and its gradient with respect to the parameters is computed by differentiating that Hessian once more. With ForwardDiff, these third derivatives are Inf or NaN when the log-density contains a division by a small number, even though the derivatives themselves are finite and well within the range of Float64.

For example, in a latent-variable model with an exponential decay and multiplicative noise, the scale of the noise becomes small for observations at later times:

using ForwardDiff, Distributions

a, ω, σ = 30.0, 1.0, 0.3
t = [1.0, 11.0]
y = [0.0014, 4.7e-47]

decay(k, z) = @. a * exp(-k * exp(z) * t)
function logjoint(k, z)
    μ = decay(k, z)
    return sum(logpdf.(Normal.(μ, σ .* μ), y)) + logpdf(Normal(0, ω), z)
end

# Laplace approximation of log ∫ p(y, z | k) dz at the mode ẑ of p(y, z | k)
∂z(k, z) = ForwardDiff.derivative(z -> logjoint(k, z), z)
∂²z(k, z) = ForwardDiff.derivative(z -> ∂z(k, z), z)
laplace(k, ẑ) = logjoint(k, ẑ) - log(-∂²z(k, ẑ)) / 2 + log(2π) / 2

function mode(k)
    lo, hi = 3.0, 3.5 # ∂z(k, lo) > 0 > ∂z(k, hi)
    for _ in 1:60
        mid = (lo + hi) / 2
        ∂z(k, mid) > 0 ? (lo = mid) : (hi = mid)
    end
    return (lo + hi) / 2
end

julia> k = 0.4;

julia> ẑ = mode(k)
3.2202943922395297

julia> laplace(k, ẑ)
102.67519619445666

julia> ForwardDiff.derivative(k -> laplace(k, ẑ), k)
-Inf

julia> h = 1e-6; (laplace(k + h, ẑ) - laplace(k - h, ẑ)) / 2h
-386.2937856666804

Cause

logpdf(Normal(μ, σ), y) computes (y - μ) / σ, here with σ ≈ 1.3e-47 for the second observation. For dual numbers x and y, the partials of x / y are computed as inv(bval) * a - (aval / (bval * bval)) * b:

https://github.com/JuliaDiff/ForwardDiff.jl/blob/v1.4.6/src/partials.jl#L99-L101

With nested duals, bval is itself a dual number, so aval / (bval * bval) divides by a dual number again, and the n-th derivative divides by y^(2^n). In the example, the third derivative divides by σ^8 ≈ (1.3e-47)^8, which underflows to 0.

The same happens without any statistics:

julia> f(x) = x / (x^2 + x^3); # = 1 / (x + x^2)

julia> g(x) = 1 / (x + x^2);

julia> d3(f, x) = ForwardDiff.derivative(x -> ForwardDiff.derivative(x -> ForwardDiff.derivative(f, x), x), x);

julia> for x in (1e-20, 1e-30, 1e-40)
           println("x = $x: ", d3(f, x), " vs ", d3(g, x))
       end
x = 1.0e-20: -6.000712508240443e80 vs -6.000000000000001e80
x = 1.0e-30: -Inf vs -5.999999999999998e120
x = 1.0e-40: NaN vs -6.000000000000002e160

Even at x = 1e-20 the third derivative of f has a relative error of 1e-4. The derivative of g is correct since a non-dual numerator is handled without the square (-(divv / v) * partials(y)):

https://github.com/JuliaDiff/ForwardDiff.jl/blob/v1.4.6/src/dual.jl#L539-L543

Computing the partials as (a - (aval / bval) * b) / bval, i.e., dividing by bval twice instead of by bval * bval, gives -386.2936383633966 in the first example and -6.000000000000001e80, -6.000000000000002e120 and -6.000000000000005e160 in the second one.

ForwardDiff v1.4.6, Distributions v0.25.132, Julia 1.13.2.

Lingua principale
Julia
Stelle
1k
Fork
165
Merge medio
8g 22h
PR unite (30g)
3

Preparare l'ambiente

Questo progetto non fornisce container di sviluppo, Dockerfile né guida per i contributori, quindi l'ambiente è a tuo carico: parti dal suo README e consulta la nostra guida al primo contributo per i passaggi generali.

Come iniziare

  1. Leggi tutta la issue e poi la guida ai contributi del progetto.
  2. Commenta sulla issue per dire che te ne occupi tu — evita che due persone facciano lo stesso lavoro.
  3. Fai un fork del repository e lavora su un branch.
  4. Apri una pull request che faccia riferimento al numero della issue.

Altre issue di JuliaDiff/ForwardDiff.jl

Tutte le issue di JuliaDiff/ForwardDiff.jl

Issue simili

Altre issue su Julia

Ricevi le nuove issue nella tua casella

Un breve riepilogo di issue GitHub adatte ai principianti.