# Discrepancy of ODE Sensitivity Analysis paper results with benchmarks

**URL:** <https://discourse.julialang.org/t/discrepancy-of-ode-sensitivity-analysis-paper-results-with-benchmarks/129853>\
**Category:** Numerics\
**Tags:** benchmark, ode, sciml\
**Created:** [June 13, 2025, 9:18am UTC](https://discourse.julialang.org/t/discrepancy-of-ode-sensitivity-analysis-paper-results-with-benchmarks/129853 "2025-06-13T09:18:10Z")\
**Posts on this page:** 12\
**Page:** 2

<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 16, 2025, 12:37pm UTC](https://discourse.julialang.org/t/discrepancy-of-ode-sensitivity-analysis-paper-results-with-benchmarks/129853/22 "2025-06-16T12:37:20Z")

</div>

Can you make a PR? That’s the easiest way to explain the change

---

<div class="post-metadata">

**Author:** ![mra](https://avatars.discourse-cdn.com/v4/letter/m/d78d45/32.png) [@mra](https://discourse.julialang.org/u/mra)\
**Post date:** [June 16, 2025, 12:42pm UTC](https://discourse.julialang.org/t/discrepancy-of-ode-sensitivity-analysis-paper-results-with-benchmarks/129853/23 "2025-06-16T12:42:44Z")

</div>

Yes, I can. Here I have also a comparison for N=2:8.

Unchanged code:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/3/8/3894cee0e23cb1cd348b6bc14c47d6f81f32e8a6.png)

Fixed code:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/f/7/f7e73668554d7a44065f0c90c7ee1e9f249a7919.png)

Code fix:

```julia
csa = map(csan) do n
  bfun, b_u0, b_p, brusselator_jac, brusselator_comp = makebrusselator!(PROBS, n)
  @time ts = map(ADJOINT_METHODS[1:2end÷3]) do alg
    @info "Running $alg ($(SciMLSensitivity.alg_autodiff(alg) ? "AD" : "user")-jac) for N = $n"
    f = SciMLSensitivity.alg_autodiff(alg) ? bfun : ODEFunction(bfun, jac=brusselator_jac)
    solver = Rodas5(autodiff=false)
    @time diffeq_sen_l2(f, b_u0, tspan, b_p, bt, solver; sensalg=alg, tols...)
    t = @elapsed diffeq_sen_l2(f, b_u0, tspan, b_p, bt, solver; sensalg=alg, tols...)
    return t
  end
  @show n,ts
  ts
end

```

The difference is huge, since now Quadrature CASA user-Jacobian is the fastest method and previously was the 3rd slowest. And it make sense that it should be the fastest, since computing user-Jacobian is much easier than running autodiff to get the Jacobian.

---

<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 16, 2025, 12:47pm UTC](https://discourse.julialang.org/t/discrepancy-of-ode-sensitivity-analysis-paper-results-with-benchmarks/129853/24 "2025-06-16T12:47:45Z")

</div>

> [@mra](#):
>
> since computing user-Jacobian is much easier than running autodiff to get the Jacobian.

No, it shouldn’t be in general because it requires constructing the whole Jacobian is required for the user one vs just the vjp. But the key here is that previous versions would not have specialized on bandedness and sparsity, which it would not have done before the DI change

---

<div class="post-metadata">

**Author:** ![mra](https://avatars.discourse-cdn.com/v4/letter/m/d78d45/32.png) [@mra](https://discourse.julialang.org/u/mra)\
**Post date:** [June 16, 2025, 1:53pm UTC](https://discourse.julialang.org/t/discrepancy-of-ode-sensitivity-analysis-paper-results-with-benchmarks/129853/25 "2025-06-16T13:53:41Z")

</div>

> <https://github.com/SciML/SciMLBenchmarksOutput/pull/9>
>
> Function f was not used and because of that user-defined brusselator\_jac was not… used either, which yielded misleading results.

I have created the PR

---

<div class="post-metadata">

**Author:** ![mra](https://avatars.discourse-cdn.com/v4/letter/m/d78d45/32.png) [@mra](https://discourse.julialang.org/u/mra)\
**Post date:** [June 16, 2025, 2:12pm UTC](https://discourse.julialang.org/t/discrepancy-of-ode-sensitivity-analysis-paper-results-with-benchmarks/129853/26 "2025-06-16T14:12:42Z")

</div>

> [@ChrisRackauckas](#):
>
> But the key here is that previous versions would not have specialized on bandedness and sparsity, which it would not have done before the DI change

Could you please expand on that? I am not sure if I understand it correctly, especially what do you mean by “previous versions”?

To my understanding, currently in SciMLBenchmarks there is no specialization regarding sparsity since `jac_prototype` is not provided i.e. everything is done on dense matrices. Furthermore adding sparsity with the use of `jac_prototype` yields better results, which means that the code where `jac_prototype` was not set, did not specialize.

Here is my implementation of adding sparsity to the SciMLBenchmark:

```julia
function jac_sparsity_adtypes(u0, func, p)
    du0 = similar(u0)
    return ADTypes.jacobian_sparsity((du, u) -> func(du, u, p, 0.0), du0, u0, TracerSparsityDetector())
end

function auto_sen_l2(f, u0, tspan, p, t, alg=Tsit5(); diffalg=ReverseDiff.gradient, kwargs...)
    jac_sparsity = jac_sparsity_adtypes(u0, f, p)
    f_ode = ODEFunction(f, jac_prototype=jac_sparsity)
  test_f(p) = begin
    prob = ODEProblem{true, SciMLBase.FullSpecialize}(f_ode,convert.(eltype(p),u0),tspan,p)
    sol = solve(prob,alg,saveat=t; kwargs...)
    sum(sol.u) do x
      sum(z->(1-z)^2/2, x)
    end
  end
  diffalg(test_f, p)
end

csa = map(csan) do n
  bfun, b_u0, b_p, brusselator_jac, brusselator_comp = makebrusselator!(PROBS, n)
  @time ts = map(ADJOINT_METHODS[1:2end÷3]) do alg
    @info "Running $alg ($(SciMLSensitivity.alg_autodiff(alg) ? "AD" : "user")-jac) for N = $n"
    jac_sparsity = jac_sparsity_adtypes(b_u0, bfun, b_p)
    f = SciMLSensitivity.alg_autodiff(alg) ? ODEFunction(bfun, jac_prototype=jac_sparsity) : ODEFunction(bfun, jac=brusselator_jac, jac_prototype=jac_sparsity)
    solver = Rodas5(autodiff=false)
    @time diffeq_sen_l2(f, b_u0, tspan, b_p, bt, solver; sensalg=alg, tols...)
    t = @elapsed diffeq_sen_l2(f, b_u0, tspan, b_p, bt, solver; sensalg=alg, tols...)
    return t
  end
  @show n,ts
  ts
end

```

And here are the results with the sparsity and bugfix implemented:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/6/c/6c373119dc6bc1164a7a8c3bf07a2c1944fb3d30.png)

It seems that ForwardDiff gained the most from the change, but still ForwardDiff is not even close to being faster than adjoint methods (for N\<=9), as it is the case in my benchmark.

---

<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 17, 2025, 7:26pm UTC](https://discourse.julialang.org/t/discrepancy-of-ode-sensitivity-analysis-paper-results-with-benchmarks/129853/27 "2025-06-17T19:26:56Z")

</div>

> [@mra](#):
>
> I have created the PR

It’s on the output repo, not the generator file [SciMLBenchmarks.jl/benchmarks/AutomaticDifferentiation/BrussScaling.jmd at master · SciML/SciMLBenchmarks.jl · GitHub](https://github.com/SciML/SciMLBenchmarks.jl/blob/master/benchmarks/AutomaticDifferentiation/BrussScaling.jmd). I’ll see if I can move it over this next weekend.

---

<div class="post-metadata">

**Author:** ![mra](https://avatars.discourse-cdn.com/v4/letter/m/d78d45/32.png) [@mra](https://discourse.julialang.org/u/mra)\
**Post date:** [June 18, 2025, 7:33am UTC](https://discourse.julialang.org/t/discrepancy-of-ode-sensitivity-analysis-paper-results-with-benchmarks/129853/28 "2025-06-18T07:33:18Z")

</div>

I have created a second pull request (now on a correct repo), if that helps.

> <https://github.com/SciML/SciMLBenchmarks.jl/pull/1276>
>
> Function \`f\` was not used and because of that user-defined \`brusselator\_jac\` was… not used either.
> 
> \## Additional context
> 
> Bug was found during the discussion on Julia Discourse (\[link to the discussion\](https://discourse.julialang.org/t/discrepancy-of-ode-sensitivity-analysis-paper-results-with-benchmarks/129853/26))

---

<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 18, 2025, 8:31am UTC](https://discourse.julialang.org/t/discrepancy-of-ode-sensitivity-analysis-paper-results-with-benchmarks/129853/29 "2025-06-18T08:31:15Z")

</div>

That helps thanks. Though I still need to do one last thing to the `remake_zero!` PR.

---

<div class="post-metadata">

**Author:** ![mra](https://avatars.discourse-cdn.com/v4/letter/m/d78d45/32.png) [@mra](https://discourse.julialang.org/u/mra)\
**Post date:** [June 18, 2025, 10:46am UTC](https://discourse.julialang.org/t/discrepancy-of-ode-sensitivity-analysis-paper-results-with-benchmarks/129853/30 "2025-06-18T10:46:07Z")

</div>

I have produced two additional plots that seem to explain the origin of the rest of the differences between my benchmark and SciMLBenchmark.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/6/7/671f69af2cd231ada21a5235a3bd1dabd439f5fb.png)  
Plot 1: SciMLBenchmark code with added sparsity and `brusselator_2d_loop` function substituted with my implementation of Brusselator

 ![image](https://global.discourse-cdn.com/julialang/original/3X/1/3/139893fb6b1da0ba032dece9ba42964c3f092703.png)

Plot 2: My benchmark with the parameters (timespan, tolerances, saveat etc.) set to the same values as in SciMLBenchmark. Entries with ADJ are using `SciMLSensitivity.adjoint_sensitivities` instead of `DifferentiationInterface.gradient` to get the sensitivities.

To summarize:

1. ForwardDiff in my benchmark was scaling better due to provided sparse `jac_prototype` and ForwardDiff was the algorithm that gained the most performance from this change.
2. In my benchmark ForwardDiff was faster for N\<10 because:  
2.1 My implementation of brusselator was slower for all algorithms except ForwardDiff (the speed remained the same) and Interpolating CASA AD-v^T seeding (was a bit faster in comparison to SciMLBenchmark implementation).  
2.2 My choice of parameters (timespan, tolerances, etc.) was promoting 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 18, 2025, 10:51am UTC](https://discourse.julialang.org/t/discrepancy-of-ode-sensitivity-analysis-paper-results-with-benchmarks/129853/31 "2025-06-18T10:51:49Z")

</div>

> [@mra](#):
>
> ForwardDiff in my benchmark was scaling better due to provided sparse `jac_prototype` and ForwardDiff was the algorithm that gained the most performance from this change.

We should add a ForwardDiff (Sparse) to the SciMLBenchmark version. Definitely an interesting observation of something that has changed in the ecosystem since 2019.

> [@mra](#):
>
> 2.1 My implementation of brusselator was slower for all algorithms except ForwardDiff (the speed remained the same) and Interpolating CASA AD-v^T seeding (was a bit faster in comparison to SciMLBenchmark implementation).  
> 2.2 My choice of parameters (timespan, tolerances, etc.) was promoting ForwardDiff.

Yeah this can be hard to account for in a simple plot. That’s why I generally say the cutoff is somewhere between like 10 to 100, because it’s problem-dependent. The scaling difference is clear it’s just unclear where you should switch, and this is really just one of many examples to try to justify the switch point.

---

<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 18, 2025, 12:42pm UTC](https://discourse.julialang.org/t/discrepancy-of-ode-sensitivity-analysis-paper-results-with-benchmarks/129853/32 "2025-06-18T12:42:24Z")

</div>

I don’t understand how providing the Jacobian prototype can change benchmark results if the sparsity detection step is not measured anyway?

---

<div class="post-metadata">

**Author:** ![mra](https://avatars.discourse-cdn.com/v4/letter/m/d78d45/32.png) [@mra](https://discourse.julialang.org/u/mra)\
**Post date:** [June 18, 2025, 2:01pm UTC](https://discourse.julialang.org/t/discrepancy-of-ode-sensitivity-analysis-paper-results-with-benchmarks/129853/33 "2025-06-18T14:01:19Z")

</div>

Below is the link to the `ODEFunction`, which takes as parameter `jac_prototype`

> **[ODE Problems · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/types/ode_types/#SciMLBase.ODEFunction)**
>
> Documentation for DifferentialEquations.jl.

To my understanding providing sparsity pattern via `jac_prototype` accelerates computation as not all entries of the jacobian need to be computed and all operations on the jacobian operate on a sparse matrix, shaving of computation.

Sparsity detection step is just to automatically determine the sparsity pattern of the jacobian and it could be done ahead of time, since it is done only once (in particular when pattern is known it can be provided manually).

Therefore the pattern detection itself is not a part of the benchmark. But all the operations done by sensitivity analysis and ODE solver are. If `jac_prototype` is provided, some of those operations will be done on sparse matrices instead of dense ones, therefore speeding up the sensitivity analysis.

[Previous page](https://discourse.julialang.org/t/discrepancy-of-ode-sensitivity-analysis-paper-results-with-benchmarks/129853.md?page=1)