Third-order derivatives of `x / y` are `Inf` or `NaN` for small `y`, e.g. in gradients of Laplace approximations
维护者通常 1 天内回复
评估
调研方向
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,通用步骤见我们的新手贡献指南。
从这里开始
- 先读完整个 Issue,再读项目的贡献指南。
- 在 Issue 下留言说明你要接手 —— 这能避免两个人做同样的事。
- Fork 仓库,在一个分支上完成修改。
- 提交 Pull Request,并在描述里引用这个 Issue 编号。
JuliaDiff/ForwardDiff.jl 的其他 Issue
-
missing floatmin, floatmax, and issubnormal for Dual可能已有人在做 @devmotion 于 30 天前认领。 未关闭
难度 2/5 1-3 小时 新手友好度 68/100
JuliaDiff/ForwardDiff.jl#852 · 6 条评论 ·
维护者通常 1 天内回复
-
`min` does not return the minimum with respect to `isless` when values tie可能已有人在做 @devmotion 于 4 天前认领。 未关闭
难度 3/5 1-2 天 新手友好度 55/100
JuliaDiff/ForwardDiff.jl#868 ·
维护者通常 1 天内回复
-
`fma` and `muladd` ignore NaN-safe mode可能已有人在做 @devmotion 于 15 天前认领。 未关闭
难度 3/5 1-2 天 新手友好度 55/100
JuliaDiff/ForwardDiff.jl#860 ·
维护者通常 1 天内回复
-
Enable NaN-safe mode by default可能已有人在做 @devmotion 于 17 天前认领。 未关闭
难度 5/5 一周以上 新手友好度 25/100
JuliaDiff/ForwardDiff.jl#857 · 1 个 reaction ·
维护者通常 1 天内回复
-
难度 4/5 3-5 天 新手友好度 52/100
JuliaDiff/ForwardDiff.jl#854 · 1 条评论 ·
维护者通常 1 天内回复
查看 JuliaDiff/ForwardDiff.jl 的全部 Issue
相似的 Issue
-
documentation invalid
难度 1/5 1-3 小时 新手友好度 88/100
JuliaPhysics/BeamletOptics.jl#190 ·
维护者通常 2 天内回复
-
mul! with a Symmetric or Hermitian factor throws for a matrix view whose columns are not contiguous未关闭
难度 2/5 1-3 小时 新手友好度 72/100
JuliaLang/LinearAlgebra.jl#1829 ·
维护者通常 1 天内回复
-
难度 2/5 1-3 小时 新手友好度 78/100
MagneticResonanceImaging/MRIReco.jl#320 ·
维护者通常 1 天内回复
-
难度 2/5 1-3 小时 新手友好度 78/100
rfordatascience/tidytuesday#1098 ·
-
documentation
难度 2/5 1-3 小时 新手友好度 70/100
trixi-framework/TrixiParticles.jl#1364 ·
维护者通常 2 天内回复