Enable NaN-safe mode by default

.github/workflows/ci.yml

@@ -15,7 +15,7 @@ permissions:
 
 jobs:
   test:
-    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - NaN-safe mode ${{ matrix.nansafe && 'enabled' || 'disabled' }} - ${{ github.event_name }}
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - unsafe mode ${{ matrix.unsafe && 'enabled' || 'disabled' }} - ${{ github.event_name }}
     runs-on: ${{ matrix.os }}
     strategy:
       fail-fast: false
@@ -29,7 +29,7 @@ jobs:
           - ubuntu-latest
           - windows-latest
           - macOS-latest
-        nansafe:
+        unsafe:
           - true
           - false
     steps:
@@ -38,10 +38,10 @@ jobs:
         with:
           version: ${{ matrix.version }}
       - uses: julia-actions/cache@v3
-      - name: Set NaN-safe mode preference
+      - name: Set unsafe mode preference
         run: |
           open(joinpath(pwd(), "Project.toml"), "a") do io
-            println(io, "\n[preferences.ForwardDiff]\nnansafe_mode = ${{ matrix.nansafe }}\n")
+            println(io, "\n[preferences.ForwardDiff]\nunsafe_mode = ${{ matrix.unsafe }}\n")
           end
         shell: julia --color=yes {0}
       - uses: julia-actions/julia-buildpkg@v1

Project.toml

@@ -1,6 +1,6 @@
 name = "ForwardDiff"
 uuid = "f6369f11-7733-5829-9624-2563aa707210"
-version = "1.4.7"
+version = "1.5.0"
 
 [deps]
 CommonSubexpressions = "bbf7d656-a473-5ed7-a52c-81e309532950"

docs/src/user/advanced.md

@@ -117,41 +117,36 @@ Note that it is usually best to pick a chunk size which divides evenly into the
 dimension. Otherwise, ForwardDiff has to construct and utilize an extra "remainder" chunk to
 complete the calculation.
 
-## Fixing NaN/Inf Issues
+## NaN-Safe Mode
 
-ForwardDiff's default behavior is to return `NaN` for undefined derivatives (or otherwise
-mirror the behavior of the function in `Base`, if it would return an error). This is
-usually the correct thing to do, but in some cases can erroneously "poison" values which
-aren't sensitive to the input and thus cause ForwardDiff to incorrectly return `NaN` or
-`Inf` derivatives. For example:
+A dual number whose perturbation component is zero is unperturbed: it must not propagate
+derivative information, regardless of the derivative of the function it is passed to.
+Naively propagating derivatives by multiplying them with the perturbation component
+violates this property whenever the derivative is not finite, since `0.0 * Inf` and
+`0.0 * NaN` are `NaN`. ForwardDiff therefore checks whether each perturbation component is
+zero before propagating derivative information ("`NaN`-safe mode"):
 
 ```julia-repl
-# the dual number's perturbation component is zero, so this
-# variable should not propagate derivative information
 julia> log(ForwardDiff.Dual{:tag}(0.0, 0.0))
-Dual{:tag}(-Inf,NaN) # oops, this NaN should be 0.0
-```
-
-Here, ForwardDiff computes the derivative of `log(0.0)` as `NaN` and then propagates
-this derivative by multiplying it by the perturbation component. Usually, ForwardDiff can
-rely on the identity `x * 0.0 == 0.0` to prevent the derivatives from propagating when
-the perturbation component is `0.0`. However, this identity doesn't hold if `isnan(y)`
-or `isinf(y)`, in which case a `NaN` derivative will be propagated instead.
+Dual{:tag}(-Inf,0.0)
 
-It is possible to fix this behavior by checking that the perturbation component is zero
-before attempting to propagate derivative information, but this check can noticeably
-decrease performance (~5%-10% on our benchmarks).
+julia> ForwardDiff.gradient(v -> sqrt(v[1]) + v[2], [0.0, 1.0])
+2-element Vector{Float64}:
+ Inf
+  1.0
+```
 
-In order to preserve performance in the majority of use cases, ForwardDiff disables this
-check by default. If your code is affected by this `NaN` behavior, you can enable
-ForwardDiff's `NaN`-safe mode by using the
+Similar to `@fastmath`, you can trade correctness for performance by disabling `NaN`-safe
+mode: derivatives may then be `NaN` even though they are well-defined (see the example
+below). Unsafe mode can make code faster, in particular code that is dominated by cheap
+arithmetic operations. To enable it, use the
 [Preferences.jl](https://github.com/JuliaPackaging/Preferences.jl) API to set
-the `nansafe_mode` preference to true, for example via:
+the `unsafe_mode` preference to true, for example via:
 
 ```julia-repl
 julia> using ForwardDiff, Preferences
 
-julia> set_preferences!(ForwardDiff, "nansafe_mode" => true)
+julia> set_preferences!(ForwardDiff, "unsafe_mode" => true)
 ```
 
 Note that Julia has to be restarted and ForwardDiff has to be reloaded after changing
@@ -162,16 +157,27 @@ Alternatively, you can set the preference before loading ForwardDiff, for exampl
 ```julia-repl
 julia> using Preferences, UUIDs
 
-julia> set_preferences!(UUID("f6369f11-7733-5829-9624-2563aa707210"), "nansafe_mode" => true)
+julia> set_preferences!(UUID("f6369f11-7733-5829-9624-2563aa707210"), "unsafe_mode" => true)
 
 julia> using ForwardDiff
 
 julia> log(ForwardDiff.Dual{:tag}(0.0, 0.0))
-Dual{:tag}(-Inf,0.0)
+Dual{:tag}(-Inf,NaN)
+
+julia> ForwardDiff.gradient(v -> sqrt(v[1]) + v[2], [0.0, 1.0])
+2-element Vector{Float64}:
+  Inf
+ NaN
 ```
 
+!!! compat "ForwardDiff 1.5"
+    Before ForwardDiff 1.5.0, `NaN`-safe mode was disabled by default and could be enabled by
+    setting the `nansafe_mode` preference to true. The `nansafe_mode` preference is
+    deprecated: setting it emits a warning, and `nansafe_mode = false` disables `NaN`-safe
+    mode only if the `unsafe_mode` preference is not set.
+
 In the future, we plan on allowing users and downstream library authors to dynamically
-enable [`NaN`-safe mode via the `AbstractConfig`
+toggle [`NaN`-safe mode via the `AbstractConfig`
 API](https://github.com/JuliaDiff/ForwardDiff.jl/issues/181).
 
 ## Hessian of a vector-valued function

src/ForwardDiff.jl

@@ -24,4 +24,18 @@ include("hessian.jl")
 
 export DiffResults
 
+function __init__()
+    if NANSAFE_MODE_PREFERENCE !== nothing
+        Base.depwarn(
+            "The `nansafe_mode` preference of ForwardDiff is deprecated. " *
+            "NaN-safe mode is enabled by default since ForwardDiff 1.5.0, and the preference to disable it is now called `unsafe_mode`. " *
+            "For now, `nansafe_mode = false` still disables NaN-safe mode unless `unsafe_mode` is set. " *
+            "Please remove the `nansafe_mode` preference with `Preferences.delete_preferences!(ForwardDiff, \"nansafe_mode\"; force = true)` and restart Julia.",
+            :__init__;
+            force = true,
+        )
+    end
+    return nothing
+end
+
 end # module

src/dual.jl

@@ -416,6 +416,10 @@ end
 
 Base.iszero(x::Dual) = iszero(value(x)) && iszero(partials(x))  # shortcut, equivalent to x == zero(x)
 
+@inline _isfinite_all(x::Dual) = _isfinite_all(x.value) & all(_isfinite_all, x.partials.values)
+@inline _is_strong_divisor(x::Dual) = _is_strong_divisor(x.value) & all(_isfinite_all, x.partials.values)
+@inline _is_strong_zero_type(::Type{<:Dual{<:Any,V}}) where {V} = _is_strong_zero_type(V)
+
 for pred in [:isequal, :(==)]
     @eval begin
         @define_binary_dual_op(
@@ -619,90 +623,37 @@ end
     calc_hypot(x, y, z, Tz),
 )
 
-# fma #
-#-----#
-
-@generated function calc_fma_xyz(x::Dual{T,<:Any,N},
-                                 y::Dual{T,<:Any,N},
-                                 z::Dual{T,<:Any,N}) where {T,N}
-    ex = Expr(:tuple, [:(_fma_partial(value(x), partials(y, $i), _fma_partial(value(y), partials(x, $i), partials(z, $i)))) for i in 1:N]...)
-    return quote
-        $(Expr(:meta, :inline))
-        v = fma(value(x), value(y), value(z))
-        return Dual{T}(v, $ex)
-    end
-end
-
-@inline function calc_fma_xy(x::Dual{T}, y::Dual{T}, z::Real) where T
-    vx, vy = value(x), value(y)
-    result = fma(vx, vy, z)
-    return Dual{T}(result, _mul_partials(partials(x), partials(y), vy, vx))
-end
+# fma/muladd #
+#------------#
 
-@generated function calc_fma_xz(x::Dual{T,<:Any,N},
-                                y::Real,
-                                z::Dual{T,<:Any,N}) where {T,N}
-    ex = Expr(:tuple, [:(_fma_partial(y, partials(x, $i), partials(z, $i))) for i in 1:N]...)
-    return quote
-        $(Expr(:meta, :inline))
-        v = fma(value(x), y, value(z))
-        Dual{T}(v, $ex)
-    end
+# `f` is `fma` or `muladd`
+@inline function calc_muladd_xyz(f::F, x::Dual{T}, y::Dual{T}, z::Dual{T}) where {F,T}
+    vx, vy = x.value, y.value
+    return Dual{T}(f(vx, vy, z.value), _muladd_partials(f, y.partials, x.partials, z.partials, vx, vy))
 end
 
-@define_ternary_dual_op(
-    Base.fma,
-    calc_fma_xyz(x, y, z),                         # xyz_body
-    calc_fma_xy(x, y, z),                          # xy_body
-    calc_fma_xz(x, y, z),                          # xz_body
-    Base.fma(y, x, z),                             # yz_body
-    Dual{Tx}(fma(value(x), y, z), partials(x) * y), # x_body
-    Base.fma(y, x, z),                              # y_body
-    Dual{Tz}(fma(x, y, value(z)), partials(z))      # z_body
-)
-
-# muladd #
-#--------#
-
-@generated function calc_muladd_xyz(x::Dual{T,<:Any,N},
-                                    y::Dual{T,<:Any,N},
-                                    z::Dual{T,<:Any,N}) where {T,N}
-    ex = Expr(:tuple, [:(_muladd_partial(value(x), partials(y, $i), _muladd_partial(value(y), partials(x, $i), partials(z, $i)))) for i in 1:N]...)
-    return quote
-        $(Expr(:meta, :inline))
-        v = muladd(value(x), value(y), value(z))
-        return Dual{T}(v, $ex)
-    end
+@inline function calc_muladd_xy(f::F, x::Dual{T}, y::Dual{T}, z::Real) where {F,T}
+    vx, vy = x.value, y.value
+    return Dual{T}(f(vx, vy, z), _mul_partials(x.partials, y.partials, vy, vx))
 end
 
-@inline function calc_muladd_xy(x::Dual{T}, y::Dual{T}, z::Real) where T
-    vx, vy = value(x), value(y)
-    result = muladd(vx, vy, z)
-    return Dual{T}(result, _mul_partials(partials(x), partials(y), vy, vx))
+@inline function calc_muladd_xz(f::F, x::Dual{T}, y::Real, z::Dual{T}) where {F,T}
+    return Dual{T}(f(x.value, y, z.value), _muladd_partials(f, x.partials, z.partials, y))
 end
 
-@generated function calc_muladd_xz(x::Dual{T,<:Any,N},
-                                   y::Real,
-                                   z::Dual{T,<:Any,N}) where {T,N}
-    ex = Expr(:tuple, [:(_muladd_partial(y, partials(x, $i), partials(z, $i))) for i in 1:N]...)
-    return quote
-        $(Expr(:meta, :inline))
-        v = muladd(value(x), y, value(z))
-        Dual{T}(v, $ex)
-    end
+for f in (:fma, :muladd)
+    @eval @define_ternary_dual_op(
+        Base.$f,
+        calc_muladd_xyz($f, x, y, z),                  # xyz_body
+        calc_muladd_xy($f, x, y, z),                   # xy_body
+        calc_muladd_xz($f, x, y, z),                   # xz_body
+        Base.$f(y, x, z),                              # yz_body
+        Dual{Tx}($f(x.value, y, z), x.partials * y),   # x_body
+        Base.$f(y, x, z),                              # y_body
+        Dual{Tz}($f(x, y, z.value), z.partials)        # z_body
+    )
 end
 
-@define_ternary_dual_op(
-    Base.muladd,
-    calc_muladd_xyz(x, y, z),                         # xyz_body
-    calc_muladd_xy(x, y, z),                          # xy_body
-    calc_muladd_xz(x, y, z),                          # xz_body
-    Base.muladd(y, x, z),                             # yz_body
-    Dual{Tx}(muladd(value(x), y, z), partials(x) * y), # x_body
-    Base.muladd(y, x, z),                             # y_body
-    Dual{Tz}(muladd(x, y, value(z)), partials(z))      # z_body
-)
-
 # sin/cos #
 #--------#