# How to troubleshoot ForwardDiff

**URL:** <https://discourse.julialang.org/t/how-to-troubleshoot-forwarddiff/99728>\
**Category:** Optimization (Mathematical)\
**Created:** [June 1, 2023, 1:57pm UTC](https://discourse.julialang.org/t/how-to-troubleshoot-forwarddiff/99728 "2023-06-01T13:57:09Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![johnb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnb/32/44115_2.png) [@johnb](https://discourse.julialang.org/u/johnb)\
**Post date:** [June 1, 2023, 1:57pm UTC](https://discourse.julialang.org/t/how-to-troubleshoot-forwarddiff/99728/1 "2023-06-01T13:57:10Z")

</div>

I’ve been trying to write a loss function for a PDE that I’m solving in order to identify parameters and I would like to use the gradient of my loss function in parameter optimization. Unfortunately, ForwardDiff returns a zero vector regardless of what the parameter values are. Numerical gradient determination is very slow, so I would prefer not to use that if I can avoid it.

Is there a way to get ForwardDiff to return where it’s having issues, or do I need to kind of guess and check potential problem areas?

Thanks!

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [June 1, 2023, 1:59pm UTC](https://discourse.julialang.org/t/how-to-troubleshoot-forwarddiff/99728/2 "2023-06-01T13:59:05Z")

</div>

Do you have an if statement that can cause your loss to return something that does not depend on the parameters?

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 1, 2023, 2:58pm UTC](https://discourse.julialang.org/t/how-to-troubleshoot-forwarddiff/99728/3 "2023-06-01T14:58:26Z")

</div>

What do you mean by “numerically”? Querying the gradient with finite differences or with ForwardDiff has a similar computational complexity

---

<div class="post-metadata">

**Author:** ![johnb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnb/32/44115_2.png) [@johnb](https://discourse.julialang.org/u/johnb)\
**Post date:** [June 1, 2023, 3:06pm UTC](https://discourse.julialang.org/t/how-to-troubleshoot-forwarddiff/99728/4 "2023-06-01T15:06:07Z")

</div>

I have a couple layers of error checking in case of solver errors as well as a parameter rescaling system.

For the parameter rescaling, at function compilation, I generate a list of parameters which I want to fit on a log scale, and use an if statement to either take `exp10()` of the input parameter or leave it as is before using it to remake my problem. I have tried bypassing this code but still obtained a zero vector.

In terms of the problem error handling, first I use a try/catch around the solution generation, and then after a solution return from the solver, I use the solution retcode to determine whether the solve was successful. In both cases, infinite loss is returned

```julia
try
    sol = solve(my_problem)
    loss = my_loss_function(sol)
catch e
    loss = Inf
end

```

and my loss function contains

```julia
if sol.retcode != :Success
    loss = Inf
else
    loss = rest_of_loss_code      

```

where the actual loss code is basically just a convoluted version of sum of squares since I can only obtain experimental data for one boundary of my PDE system.

Removing the first checkpoint does not change the result.

The second checkpoint is recommended here:  
[https://docs.sciml.ai/DiffEqParamEstim/stable/methods/optimization\_based\_methods/#Note-About-Loss-Functions](https://docs.sciml.ai/DiffEqParamEstim/stable/methods/optimization_based_methods/#Note-About-Loss-Functions)

While I am not using DiffEqParamEstim, I could not think of a better way to ensure I do not unintentionally optimize for an improperly generated solution.

However, removing that check revealed what may be the issue. My loss code begins by manipulating the continuous solution into a vector of discrete average values in order to match my experimental data collection strategy. It would seem that this approach, which relies on the code below, is incompatible with forward diff.

```julia
for j = 1:maxfracs
    CV = startCV[j]
    push!(FractionData, mean(sol(CV2time(CV):CV2time(CV+fracvol[j]), 10.0, dv = c(t,x))))
end

```

In particular, retrieval of the vector of solution values by `sol(t, x, dv=dv)` may be causing the issue. I will update once I figure out what specifically is causing the error.

Thank you for the tip to check if/else statements!

---

<div class="post-metadata">

**Author:** ![johnb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnb/32/44115_2.png) [@johnb](https://discourse.julialang.org/u/johnb)\
**Post date:** [June 1, 2023, 3:19pm UTC](https://discourse.julialang.org/t/how-to-troubleshoot-forwarddiff/99728/5 "2023-06-01T15:19:31Z")

</div>

By “numerically”, I mean using `FiniteDiff.finite_difference_gradient(error_startup,default_params)` or `grad(central_fdm(5,1),error_startup,default_params)`.

These two return different results, and the fifth order central fdm result which I trust more takes approximately a minute to complete vs 0.3s for my error function. Given that I have 14 parameters, this makes sense, but I am hoping that a proper ForwardDiff gradient will be a bit faster.

Given what I now know about the ForwardDiff result, it seems like it was not actually “doing” anything, which may have been why the time to return something was so low. I will need to fix that issue before I figure out how long each gradient evaluation actually takes to run.

---

<div class="post-metadata">

**Author:** ![johnb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnb/32/44115_2.png) [@johnb](https://discourse.julialang.org/u/johnb)\
**Post date:** [June 1, 2023, 3:49pm UTC](https://discourse.julialang.org/t/how-to-troubleshoot-forwarddiff/99728/6 "2023-06-01T15:49:39Z")

</div>

Update: It seems like the PDE solve is failing, which causes the problem with solution retrieval.

This is strange, because `error_(default_params)` returns a reasonable value and when I `@show sol.retcode` it shows `sol.retcode = SciMLBase.ReturnCode.Success`

When I include `@show sol.retcode` in the error function passed through ForwardDiff, it returns `sol.retcode = SciMLBase.ReturnCode.DtLessThanMin`.

This explains the zero vector error at least.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 1, 2023, 10:52pm UTC](https://discourse.julialang.org/t/how-to-troubleshoot-forwarddiff/99728/8 "2023-06-01T22:52:49Z")

</div>

I’m curious about the edit. I was going to take a look and then the code was deleted. Should I not?

---

<div class="post-metadata">

**Author:** ![johnb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnb/32/44115_2.png) [@johnb](https://discourse.julialang.org/u/johnb)\
**Post date:** [June 1, 2023, 11:05pm UTC](https://discourse.julialang.org/t/how-to-troubleshoot-forwarddiff/99728/9 "2023-06-01T23:05:29Z")

</div>

I found a mistake in my example code that fixed the error and made the toy model solve properly. I am still having the issue with my real system, but I have not been able to come up with a MWE yet. I’m still working on it though!

So far I’ve ruled out:

- Vector equations
- Averages and interpolations in loss function
- Multiple solves (experiments) per loss function evaluation
- Interpolations used as boundary conditions
- Solver choice
- Use of push!() to construct vectors
- Loss calculation split into multiple functions
- Solver tolerance (lowered to very low value with no change)  
updates:
- PDAE system using mass matrix
- Actually stiff set of equations
- jac = true and sparse = true
- Data pulled from spreadsheet
- Chebyshev discretization

I still need to test the following features present in my real system:

- ~~PDAE system using mass matrix (my current guess)~~
- ~~Data pulled from spreadsheet (seems unlikely)~~
- ~~Actually stiff set of equations (maybe)~~
- ~~Chebyshev discretization (seems unlikely)~~

I’m not sure what else I should test, but if none of these gives a real example, I may get new ideas from the results I’ve gotten so far.

I do know that it has to do with the PDE solve before any actual loss calculations are performed since I can tell that the solver halts immediately since the only value present in `sol.t` is 0.0.

Updated with newly tested items.

---

<div class="post-metadata">

**Author:** ![johnb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnb/32/44115_2.png) [@johnb](https://discourse.julialang.org/u/johnb)\
**Post date:** [June 2, 2023, 4:21pm UTC](https://discourse.julialang.org/t/how-to-troubleshoot-forwarddiff/99728/10 "2023-06-02T16:21:23Z")

</div>

Got a working example! It seems to have something to do with this equation.

```julia
 0.0 ~ 0.1 * u(t,x)*max(0,(u(t,x) - v(t,x)))^r - (u(t,x) - v(t,x)) * ((1.1+sin(t))^r)]

```

and in particular, this term `max(0,(u(t,x) - v(t,x)))^r`

```julia
using DifferentialEquations, ModelingToolkit, MethodOfLines, DomainSets, ForwardDiff, FiniteDiff, StatsBase
# Parameters, variables, and derivatives

@parameters t, x, p, q, r
@variables u(..), v(..)
#@variables u[1:n_comp](..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2
params = Symbolics.scalarize(reduce(vcat,[p => 1.5, q => 1.2, r .=> 1.3]))
# 1D PDE and boundary conditions
eqs = [Dt(u(t, x)) ~ -q*Dx(u(t,x)) + p * Dxx(u(t, x)),
        0.0 ~ 0.1 * u(t,x)*max(0,(u(t,x) - v(t,x)))^r - (u(t,x) - v(t,x)) * ((1.1+sin(t))^r)]
bcs = [u(0, x) ~ 0.0,
        v(0,x) ~ 0.0,
        u(t, 0) ~ 0.5*exp((-(t^2))/0.01),
        v(t,0) ~ 0.5*exp((-(t^2))/0.01),
        Dx(u(t,0)) ~ 0.0,
        Dx(v(t,0)) ~ 0.0]

# Space and time domains
domains = [t ∈ Interval(0.0, 1.0),
        x ∈ Interval(0.0, 1.0)]

# PDE system

@named pdesys = PDESystem(eqs, bcs, domains, [t,x], [u(t,x), v(t,x)], params)

dx = 0.1
# Method of lines discretization
order = 2
discretization = MOLFiniteDifference([x=>dx], t; approx_order = order)
# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)
# Solve ODE problem
sol = solve(prob, KenCarp47(), saveat=0.02);
function lossfun(sol)
        output = []
        for i = 1:10
                push!(output,mean(sol((i-1)/10:i/10, 1.0, dv = u(t,x))))
        end
        outputvec = stack(output)'
        diffs = sum(abs.(outputvec.-0))
        return diffs
end
function errortest(param_vec)
    loss = []
    newprob = remake(prob, p = param_vec)
    try
            sol = solve(newprob, saveat = 0.1, KenCarp47(), reltol=10^-8, abstol=10^-8)
            if sol.retcode != :Success
                    println("anomalous solution generated")
                    loss = Inf
            else
                    diffs = lossfun(sol)
                    loss = diffs
                    
            end
    catch
            println("solve failed")
            push!(loss, Inf)
    end
    return loss
end 

default_params = [1.5 1.2 1.3]
#This works
errortest(default_params)

#This does not
ForwardDiff.gradient(errortest,default_params)

#This works
FiniteDiff.finite_difference_gradient(errortest,default_params)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 2, 2023, 4:28pm UTC](https://discourse.julialang.org/t/how-to-troubleshoot-forwarddiff/99728/11 "2023-06-02T16:28:28Z")

</div>

Try running with ForwardDiff in NaNSafe mode.

---

<div class="post-metadata">

**Author:** ![johnb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnb/32/44115_2.png) [@johnb](https://discourse.julialang.org/u/johnb)\
**Post date:** [June 2, 2023, 4:42pm UTC](https://discourse.julialang.org/t/how-to-troubleshoot-forwarddiff/99728/12 "2023-06-02T16:42:54Z")

</div>

NaNSafe mode did not seem to fix things here. I got the same warnings as before.

```julia
 Warning: dt(2.220446049250313e-16) <= dtmin(2.220446049250313e-16) at t=0.00019963280237591814, and step error estimate = 0.001174877026175027. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase .julia\packages\SciMLBase\qp2gL\src\integrator_interface.jl:599
┌ Warning: Solution has length 1 in dimension t. Interpolation will not be possible for variable u(t, x). Solution return code is DtLessThanMin.
└ @ MethodOfLines .julia\packages\MethodOfLines\dLkkv\src\interface\solution\solution_utils.jl:17
┌ Warning: Solution has length 1 in dimension t. Interpolation will not be possible for variable v(t, x). Solution return code is DtLessThanMin.
└ @ MethodOfLines .julia\packages\MethodOfLines\dLkkv\src\interface\solution\solution_utils.jl:17
anomalous solution generated
1×3 Matrix{Float64}:
 0.0 0.0 0.0

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 2, 2023, 4:44pm UTC](https://discourse.julialang.org/t/how-to-troubleshoot-forwarddiff/99728/13 "2023-06-02T16:44:48Z")

</div>

> [@johnb](#):
>
> `solve(newprob, saveat = 0.1, KenCarp47(), reltol=10^-8, abstol=10^-8)`

`solve(newprob, KenCarp47(), saveat = 0.1, reltol=10^-8, abstol=10^-8)` ?

What about another solver like Rodas5P() or FBDF?

---

<div class="post-metadata">

**Author:** ![johnb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnb/32/44115_2.png) [@johnb](https://discourse.julialang.org/u/johnb)\
**Post date:** [June 2, 2023, 4:58pm UTC](https://discourse.julialang.org/t/how-to-troubleshoot-forwarddiff/99728/14 "2023-06-02T16:58:18Z")

</div>

With the toy model, Rodas5P gives the same error as before, but FBDF is able to successfully return a gradient vector.

Unfortunately, my real system has crazy oscillations when using FBDF, so I am unable to take advantage of this solution.

I’ve tried most of the solvers on the list below, and only KenCarp47 (which isn’t on the list) avoids the oscillatory behavior.

[https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae\_solve/#OrdinaryDiffEq.jl-(Mass-Matrix)](https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/#OrdinaryDiffEq.jl-(Mass-Matrix))

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 2, 2023, 5:02pm UTC](https://discourse.julialang.org/t/how-to-troubleshoot-forwarddiff/99728/15 "2023-06-02T17:02:58Z")

</div>

> [@johnb](#):
>
> I’ve tried most of the solvers on the list below, and only KenCarp47 (which isn’t on the list) avoids the oscillatory behavior.

Is the mass matrix removed by structural simplification? Check `prob.f.mass_matrix`

---

<div class="post-metadata">

**Author:** ![johnb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnb/32/44115_2.png) [@johnb](https://discourse.julialang.org/u/johnb)\
**Post date:** [June 2, 2023, 5:08pm UTC](https://discourse.julialang.org/t/how-to-troubleshoot-forwarddiff/99728/16 "2023-06-02T17:08:57Z")

</div>

It does not seem to be removed. I’m getting an 18x18 matrix where the first nine diagonal elements are 1 for the toy model, and a larger version of the same thing for my real system.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 2, 2023, 5:19pm UTC](https://discourse.julialang.org/t/how-to-troubleshoot-forwarddiff/99728/17 "2023-06-02T17:19:03Z")

</div>

A singular mass matrix shouldn’t work with KenCarp47 and it should be throwing in that case?

---

<div class="post-metadata">

**Author:** ![johnb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnb/32/44115_2.png) [@johnb](https://discourse.julialang.org/u/johnb)\
**Post date:** [June 2, 2023, 5:22pm UTC](https://discourse.julialang.org/t/how-to-troubleshoot-forwarddiff/99728/18 "2023-06-02T17:22:33Z")

</div>

I’m not sure what’s going on then, because it solves the problem without issues, it just doesn’t work with ForwardDiff.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 2, 2023, 5:23pm UTC](https://discourse.julialang.org/t/how-to-troubleshoot-forwarddiff/99728/19 "2023-06-02T17:23:19Z")

</div>

@xtalax could you take a look?

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [June 12, 2023, 10:05pm UTC](https://discourse.julialang.org/t/how-to-troubleshoot-forwarddiff/99728/20 "2023-06-12T22:05:45Z")

</div>

I have (finally) seen this, I will take a look soon