Hacktoberfest 2026: die Issues, die Maintainer für den Oktober markiert haben – offen und einsteigerfreundlich. Hacktoberfest-Issues durchsuchen

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

Offen
  • Stimmt etwas nicht?
#866 0 Kommentare 0 Reaktionen 0 zugewiesene Personen Auf GitHub ansehen

Maintainer antworten meist innerhalb von 1 Tag

@devmotion arbeitet bereits daran.

Seit 25.10.2026.

  • #867 von @devmotion — offen

Bewertung

Schwierigkeit
3/5
Geschätzter Aufwand
1-2 Tage
Anfängerfreundlichkeit
65/100
Issue-Typ
Bug
Klarheit
Klar beschrieben
Aktivitätsstatus
Aktiv
Tech-Stack
julia
Bereich
backend, testing

Rechercherichtung

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.

Vom Indexierungsmodell aus dem Issue-Text verfasst.

Beschreibung

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.

Vorherrschende Sprache
Julia
Sterne
1k
Forks
165
Ø Merge
6 Std. 34 Min.
Gemergte PRs (30 T.)
1

Entwicklungsumgebung

Dieses Projekt bietet weder Dev-Container noch Dockerfile noch Beitragsleitfaden – die Einrichtung liegt bei Ihnen. Beginnen Sie mit der README; die allgemeinen Schritte stehen in unserem Leitfaden für den ersten Beitrag.

Erste Schritte

  1. Lesen Sie das ganze Issue und danach den Beitragsleitfaden des Projekts.
  2. Schreiben Sie ins Issue, dass Sie es übernehmen — das erspart doppelte Arbeit.
  3. Forken Sie das Repository und arbeiten Sie in einem Branch.
  4. Öffnen Sie einen Pull Request, der die Issue-Nummer nennt.

Mehr aus JuliaDiff/ForwardDiff.jl

Alle Issues in JuliaDiff/ForwardDiff.jl

Ähnliche Issues

Weitere Issues zu Julia

Neue Issues direkt in Ihr Postfach

Eine kurze Übersicht über anfängerfreundliche GitHub-Issues.