Summary
When ForwardDiff.jacobian! is called on a function containing ifelse(u[i] <= 0, …) evaluated at exactly u[i] = 0, the resulting Jacobian row depends on the chunk size. Different chunk sizes select different branches of the ifelse, producing structurally different Jacobians for identical inputs.
MWE
using ForwardDiff
function f!(out, u)
out[1] = u[3]
out[2] = ifelse(u[1] <= 0, u[1], u[2]) # kink at u[1] = 0
out[3] = -u[2] / 100
end
u = [0.0, 512000.0, -5.0]
out = zeros(3)
# Default (chunksize = 3 for length(u) == 3)
cfg_default = ForwardDiff.JacobianConfig(f!, out, u)
J_default = zeros(3, 3)
ForwardDiff.jacobian!(J_default, f!, out, u, cfg_default)
@show J_default[2, :] # [0.0, 1.0, 0.0]
# Chunksize = 1
cfg_1 = ForwardDiff.JacobianConfig(f!, out, u, ForwardDiff.Chunk{1}())
J_1 = zeros(3, 3)
ForwardDiff.jacobian!(J_1, f!, out, u, cfg_1)
@show J_1[2, :] # [0.0, 0.0, 0.0]Output:
J_default[2, :] = [0.0, 1.0, 0.0]
J_1[2, :] = [0.0, 0.0, 0.0]
Mechanism
<= on Dual{Tag, Float64, N} is lexicographic: primal first, then partials. At u[1] = 0.0:
| expression |
primal |
partials |
verdict |
Dual(0.0, (1.0,)) <= 0.0 |
tie |
1 > 0 |
false |
Dual(0.0, (0.0,)) <= 0.0 |
tie |
0 == 0 |
true |
Chunksize = 3 (single pass, multi-seed)
u = [Dual(0.0; 1,0,0), Dual(512000; 0,1,0), Dual(-5; 0,0,1)]
Dual(0.0; 1,0,0) <= 0.0 ⇒ false
⇒ out[2] = u[2] = Dual(512000; 0,1,0)
⇒ J[2, :] = [0, 1, 0]
Consistent: single evaluation, single branch, three partials from the same branch.
Chunksize = 1 (three passes, one seed each)
Pass 1 — seed u[1]:
u = [Dual(0.0; 1), Dual(512000; 0), Dual(-5; 0)]
Dual(0.0; 1) <= 0.0 ⇒ false (non-zero partial)
⇒ out[2] = u[2] = Dual(512000; 0)
⇒ J[2, 1] = 0
Pass 2 — seed u[2]:
u = [Dual(0.0; 0), Dual(512000; 1), Dual(-5; 0)]
Dual(0.0; 0) <= 0.0 ⇒ true (zero partial, equal to plain 0.0)
⇒ out[2] = u[1] = Dual(0.0; 0)
⇒ J[2, 2] = 0
Pass 3: 0.
Result: J[2, :] = [0, 0, 0]. The Jacobian is column-wise inconsistent — column 1 linearizes the false branch, column 2 linearizes the true branch.
Why this matters
The row J[2, :] = [0, 0, 0] is singular. If a downstream solver (Newton iteration on an implicit DAE, gradient-based optimizer, etc.) receives this Jacobian it loses the ability to correct the coordinate associated with that row. In the case that prompted this report, an OrdinaryDiffEq.jl composite implicit-Euler integrator gets stuck on the wrong branch of an algebraic constraint (SciML/OrdinaryDiffEq.jl#3482).
The chunksize-3 result [0, 1, 0] is also not the "mathematically correct" one-sided derivative at the kink — both are artifacts of the lexicographic <=. But [0, 1, 0] is at least non-singular and consistent across columns, so Newton/optimizer can still make progress and the failure is silent rather than catastrophic.
Question / proposal
Should chunk-size-1 Jacobian assembly make a best-effort to match the multi-seed chunk behavior at kink points? Possible approaches:
- Tie-break
<= on Dual with zero partials the same way as <= with non-zero partials (i.e., always lexicographic) — currently the zero-partial case reduces to plain-Real <= and returns true at 0==0, which is what introduces the inconsistency.
- Document that
ifelse(x <= 0, …) is non-differentiable at x = 0 and the chunk size affects which side's gradient ForwardDiff picks — and recommend users either shift the threshold or write smooth replacements.
(1) would unify the chunksize-1 and chunksize-N results and fix the silent-wrong-Jacobian failure mode. (2) is the minimum "at least document it so callers know" response.
Versions
- ForwardDiff v1.3.3
- Julia 1.12.4
Summary
When
ForwardDiff.jacobian!is called on a function containingifelse(u[i] <= 0, …)evaluated at exactlyu[i] = 0, the resulting Jacobian row depends on the chunk size. Different chunk sizes select different branches of the ifelse, producing structurally different Jacobians for identical inputs.MWE
Output:
Mechanism
<=onDual{Tag, Float64, N}is lexicographic: primal first, then partials. Atu[1] = 0.0:Dual(0.0, (1.0,)) <= 0.01 > 0Dual(0.0, (0.0,)) <= 0.00 == 0Chunksize = 3 (single pass, multi-seed)
Consistent: single evaluation, single branch, three partials from the same branch.
Chunksize = 1 (three passes, one seed each)
Pass 1 — seed u[1]:
Pass 2 — seed u[2]:
Pass 3: 0.
Result:
J[2, :] = [0, 0, 0]. The Jacobian is column-wise inconsistent — column 1 linearizes the false branch, column 2 linearizes the true branch.Why this matters
The row
J[2, :] = [0, 0, 0]is singular. If a downstream solver (Newton iteration on an implicit DAE, gradient-based optimizer, etc.) receives this Jacobian it loses the ability to correct the coordinate associated with that row. In the case that prompted this report, an OrdinaryDiffEq.jl composite implicit-Euler integrator gets stuck on the wrong branch of an algebraic constraint (SciML/OrdinaryDiffEq.jl#3482).The chunksize-3 result
[0, 1, 0]is also not the "mathematically correct" one-sided derivative at the kink — both are artifacts of the lexicographic<=. But[0, 1, 0]is at least non-singular and consistent across columns, so Newton/optimizer can still make progress and the failure is silent rather than catastrophic.Question / proposal
Should chunk-size-1 Jacobian assembly make a best-effort to match the multi-seed chunk behavior at kink points? Possible approaches:
<=onDualwith zero partials the same way as<=with non-zero partials (i.e., always lexicographic) — currently the zero-partial case reduces to plain-Real<=and returns true at 0==0, which is what introduces the inconsistency.ifelse(x <= 0, …)is non-differentiable atx = 0and the chunk size affects which side's gradient ForwardDiff picks — and recommend users either shift the threshold or write smooth replacements.(1) would unify the chunksize-1 and chunksize-N results and fix the silent-wrong-Jacobian failure mode. (2) is the minimum "at least document it so callers know" response.
Versions