JuliaSmoothOptimizers / JuliaSmoothOptimizers/ADNLPModels.jl

SparseReverseADHessian records its ReverseDiff tape on uninitialised memory

Open
#383 2 comments 0 reactions 0 assignees View on GitHub
Dominant language
Julia
Stars
44
Forks
19
PR merge metrics
No merged PRs in 30d

Description

# `SparseReverseADHessian` records its ReverseDiff tape on uninitialised memory

## Summary

`SparseReverseADHessian` (and the `:optimized` predefined backend that uses it)
allocates its tape input as `Vector{ForwardDiff.Dual{…}}(undef, nvar)` and then
runs the model function on it while building the ReverseDiff tape. Uninitialised
memory can legally hold non-finite `Float64`s (e.g. when the allocator hands back
a just-freed `-Inf` bounds array). If the objective/constraints contain a
domain-restricted primitive — `cos`/`sin` (via `Base.sincos`), `sqrt`, `log`, … —
tape construction throws instead of returning a backend:

```
DomainError with -Inf:
sincos(x) is only defined for finite x.
```

The model is never solved; `ADNLPModel!` / the backend constructor itself fails.

## Where

`src/sparse_hessian.jl`, both `SparseReverseADHessian` methods
(`v0.8.13`, lines ~193 and ~207):

```julia
z = Vector{ForwardDiff.Dual{tagf, T, 1}}(undef, nvar) # not initialised
f_tape = ReverseDiff.GradientTape(f, z) # executes f(z)
...
zψ = Vector{ForwardDiff.Dual{tagψ, T, 1}}(undef, nvar) # not initialised
yψ = fill!(similar(zψ, ncon), zero(T))
cfgψ = ReverseDiff.compile(ReverseDiff.GradientTape(ψ, (zψ, yψ), ψ_tape)) # executes ψ(zψ, yψ)
```

`SparseADHessian` (forward mode) is **not** affected: `ForwardDiff.GradientConfig`
only sizes buffers, it never runs the function.

## Reproducer

`ADNLPModels 0.8.13`, `ForwardDiff 1.4.5`, `ReverseDiff 1.17.0`, Julia 1.12.

```julia
using ADNLPModels, ForwardDiff, ReverseDiff, LinearAlgebra, SparseArrays

# min sum(x) s.t. cos(x[1]) - 1 = 0 , x ∈ Rⁿ (x unbounded: no finite fallback)
f(x) = sum(x)
c!(cx, x) = (cx[1] = cos(x[1]) - 1; cx)
n, m = 4, 1
Hpat = sparse(1:n, 1:n, trues(n), n, n)

# --- 1. deterministic: the exact internal snippet, with zψ holding a value that
# `undef` legally may hold (here: -Inf, as a recycled `-Inf` bounds array).
ψ = (x, u) -> (tmp = similar(x, length(u)); c!(tmp, x); dot(tmp, u))
tagψ = ForwardDiff.Tag{typeof(ψ), Float64}
zψ = Vector{ForwardDiff.Dual{tagψ, Float64, 1}}(undef, n)
fill!(zψ, ForwardDiff.Dual{tagψ}(-Inf, 0.0))
yψ = fill!(similar(zψ, m), 0.0)
ReverseDiff.GradientTape(ψ, (zψ, yψ), ReverseDiff.GradientConfig((zψ, yψ)))
# -> ERROR: DomainError with -Inf: sincos(x) is only defined for finite x.

# --- 2. through the public constructor, no touching of internals: just make the
# small-object pool contain the -Inf pattern first (as a freed bounds
# array would), then call it. Fires ~80-90% of the time.
NEGINFDUAL = ForwardDiff.Dual{tagψ}(-Inf, -Inf)
poison() = (for _ in 1:5000; b = fill(NEGINFDUAL, n); Base.donotdelete(b); end; GC.gc(); GC.gc())
fails = 0
for _ in 1:100
poison()
try ADNLPModels.SparseReverseADHessian(n, f, m, c!, Hpat) catch; global fails += 1 end
end
@show fails # e.g. 81 / 100

# --- 3. contrast: forward-mode Hessian on the same model, same poisoning: 0 / 100.
fwd = 0
for _ in 1:100
poison()
try ADNLPModels.SparseADHessian(n, f, m, c!, Hpat) catch; global fwd += 1 end
end
@show fwd # 0 / 100
```

## Real-world hit

CTDirect.jl builds NLP variable bounds as `-Inf .* ones(nvar)`; those arrays are
freed and their storage is recycled into the `undef` `zψ` above. A free-final-time
optimal-control problem with a rotation-matrix dynamics (`cos(θ)`, `sin(θ)`, `θ`
unbounded) then fails to build on CI — deterministically on Julia 1.12, flakily on
1.10 — depending only on which dependency versions were resolved that run.
(control-toolbox/CTDirect.jl CI, `test/ci/test_all_ocp.jl` `:moonlander`.)

## Suggested fix

Initialise the tape inputs before recording — from `x0` (already a kwarg,
defaulting to `rand(nvar)`), or at least `fill!` the primal parts with a finite
value:

```julia
z = Vector{ForwardDiff.Dual{tagf, T, 1}}(undef, nvar)
@inbounds for i in 1:nvar
z[i] = ForwardDiff.Dual{tagf}(x0[i], zero(T))
end
```

and likewise for `zψ` / `lz`. Recording on a well-defined point also makes tape
compilation deterministic.

Contributor guide

Open the contributing guide

Research direction

Start in src/sparse_hessian.jl at both SparseReverseADHessian methods, then run the provided reproducer with the domain-restricted cos model. Ensure the tape inputs z, zψ, and lz are initialized from the available starting values before tape construction; done means the constructor succeeds reliably with recycled -Inf memory and the existing forward-mode behavior remains unaffected.

Written by the indexing model from the issue text.

Assessment

Tech stack
julia
Domain
backend-api-design
Issue type
Bug
Difficulty
3/5
Estimated time
1-2 days
Activity status
Active
Clarity
Clearly specified
Newbie friendliness
84/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.