Why does ForwardDiff.jl give me all-zero Jacobian matrix?

General Usage
, , , ,
1

Hi guys,

I am using one Julia-based model from https://github.com/LandEcosystems/Sindbad.jl, and tries to use ForwardDiff.jl to calculate the Jacobian of the cost function to the optimized parameters. I run the pipeline (simulation, optimization, and Jacobian calculation) on different sites (different locations). but I don’t know why ForwardDiff.jl gives me all-zero Jacobian matrix at some sites (not all sites run, but some sites…like 5 out of 40 sites in total). below are some details. Any suggestions/feedbacks would be greatly appreciated!

Thank you!

=======

I tried to calculate the parameter uncertainty by calculate the Hessian matrix (J^T*J, as the covariance matrix), then the std of the diagonal elements as the uncertainty. J is of the shape as (n_simulation, n_parameters). For example, if we optimize 10 parameters to simulate a variable of time series length 1000, then the size of J would be (1000, 10). We build a function taking input as the optimized parameter, returning the time series of simulation (size (n_parameter,1)). then use ForwardDiff.jacobian to calculate J matrix. The problem is at some sites, J is all zero…but why?

If a minimal example is required, I can try to write a pseudo code here…but due to the requirements of usage of other input data in the simulations, I am afraid a minimal reproducible example is not feasible here…

2

Without a reduced example it is hard to say definitively, but an all-zero ForwardDiff Jacobian usually points to one of two possibilities:

  1. the model output is genuinely locally insensitive to those parameters at those sites, or
  2. the ForwardDiff.Dual numbers are being dropped somewhere inside the simulation pipeline.

A useful first check is to compare ForwardDiff against a simple finite-difference perturbation at one of the problematic sites. For example, perturb one parameter at a time:

y0 = f(popt)

j = 1
ϵ = 1e-4
p1 = copy(popt)
p1[j] += ϵ

y1 = f(p1)

norm(y1 - y0) / ϵ

You could repeat this for several parameters. If these finite-difference sensitivities are also zero or extremely small, then the site may be in a flat or saturated regime, or the optimized parameters may not actually be affecting the simulated output for that site.

If the finite-difference sensitivity is nonzero but

ForwardDiff.jacobian(f, popt)

is all zeros, then I would suspect an AD issue. Common causes are parameters being converted back to Float64, copied into an Array{Float64}, assigned into fields typed as Vector{Float64}, or otherwise losing their ForwardDiff.Dual type inside the model.

A quick diagnostic is to check whether the dual type survives through the function:

J = ForwardDiff.jacobian(p -> begin
    @show eltype(p)
    y = f(p)
    @show eltype(y)
    y
end, popt)

If eltype(p) is a ForwardDiff.Dual type but eltype(y) is Float64, then derivative information is being dropped somewhere in the simulation.

I would also check for allocations like

zeros(Float64, ...)
Array{Float64}(undef, ...)

inside the simulation code. For ForwardDiff compatibility, these often need to be written generically, for example using eltype(p) or similar, so that arrays can store dual numbers.

One additional point: if you are using J'J for parameter uncertainty, make sure the scaling and statistical assumptions are appropriate. In many least-squares settings, the covariance approximation is proportional to inv(J'J), often scaled by an estimate of the residual variance, rather than just J'J itself. Also, if J is zero or nearly rank-deficient, the corresponding uncertainty estimate will be ill-conditioned or non-identifiable.