Skip to content

gradient! leaves part of the result unwritten when f does not depend on its argument #847

Description

@devmotion

When f does not depend on its argument but carries a perturbation from an enclosing differentiation, ForwardDiff.gradient! writes only part of the result and leaves the rest exactly as it found it. The gradient is zero in this situation, so every entry should be written.

using ForwardDiff, DiffResults
x = [1.0, 2.0, 3.0]

ForwardDiff.derivative(1.0) do a
    out = fill(a * 111.0, 3)                    # a buffer with recognisable junk
    ForwardDiff.gradient!(out, z -> a * 2.0, x)
    @show ForwardDiff.value.(out)               # [0.0, 111.0, 111.0], should be [0.0, 0.0, 0.0]
    return zero(a)
end

The allocating form has the same defect against the buffer it just allocated, so it returns uninitialized memory. That is often zero, which hides the bug, but not always:

julia> ForwardDiff.derivative(1.0) do a
           @show ForwardDiff.gradient(z -> a * 2.0, x)
           return zero(a)
       end
gradient = Dual{}[Dual(0.0,0.0), Dual(2.3611913415e-314,2.361521164e-314), Dual(5.0e-324,2.260386021e-314)]

The DiffResult form errors instead:

ForwardDiff.derivative(1.0) do a
    r = DiffResults.GradientResult(fill(a * 0.0, 3))
    ForwardDiff.gradient!(r, z -> a * 2.0, x)
    return zero(a)
end
# ERROR: MethodError: Cannot `convert` an object of type Dual{…} to an object of type Float64

The StaticArray path is correct and returns [0.0, 0.0, 0.0], because its @generated extract_gradient is built from length(x) and uses the indexed partials(T, y, i).

Cause

src/gradient.jl#L65-L72:

extract_gradient!(::Type{T}, result::AbstractArray, y::Real) where {T} = fill!(result, zero(y))
function extract_gradient!(::Type{T}, result::AbstractArray, dual::Dual) where {T}
    idxs = structural_eachindex(result)
    for (i, idx) in zip(1:npartials(dual), idxs)
        result[idx] = partials(T, dual, i)
    end
    return result
end

A Dual{S,V,N} with S ≺ T carries no T-perturbation, so its gradient is zero everywhere and the ::Real method above is what it needs — but it is a Dual, so it takes the second method, where npartials reports N for the wrong layer. Whenever that N is smaller than length(x) the tail of the result is never written.

L59-L63 shares the dispatch and then passes partials(T, dual) to DiffResults.gradient!; for this case that is zero(dual), a scalar Dual rather than a Partials, which is the same shape problem as in #846.

Suggested fix

Dispatch on Dual{T} rather than on Dual, so a value carrying no T-perturbation falls through to the ::Real method, which already fills the whole result with zeros.

Version: ForwardDiff master (v1.4.5), Julia 1.12.7, StaticArrays 1.9.19.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions