Hacktoberfest 2026:维护者为十月标记出来的 issue,仍然开放、适合新手。 浏览 Hacktoberfest issue

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

未关闭
#866 0 条评论 0 个 reaction 已指派 0 人 在 GitHub 查看

维护者通常 1 天内回复

@devmotion 已经在做这个了。

开始于 2026年10月25日。

  • #867 来自 @devmotion —— 未关闭

评估

难度
3/5
预计耗时
1-2 天
新手友好度
65/100
Issue 类型
缺陷
描述清晰度
描述清楚
活跃度
活跃
技术栈
julia
领域
backend, testing

调研方向

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.

由索引模型根据 Issue 内容生成。

描述

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.

主要语言
Julia
星标
1k
派生
165
平均合并
6 小时 34 分钟
30 天内合并 PR
1

环境准备

这个项目没有提供开发容器、Dockerfile 或贡献指南,环境需要你自己搭建:先看它的 README,通用步骤见我们的新手贡献指南。

从这里开始

  1. 先读完整个 Issue,再读项目的贡献指南。
  2. 在 Issue 下留言说明你要接手 —— 这能避免两个人做同样的事。
  3. Fork 仓库,在一个分支上完成修改。
  4. 提交 Pull Request,并在描述里引用这个 Issue 编号。

JuliaDiff/ForwardDiff.jl 的其他 Issue

查看 JuliaDiff/ForwardDiff.jl 的全部 Issue

相似的 Issue

更多 Julia Issue

把新 issue 发到你的邮箱

精选适合新手参与的 GitHub issue 摘要。