Seeking Tips for Solving an ODE with DifferentialEquations

Specific Domains Modelling & Simulations
,
1

Our package QuantumAnnealing solves ODEs that arise in design of quantum annealing algorithms. The package includes two solution methods simulate (a hand crafted solution) and simulate_de (a generic solution approach using DifferentialEquations).

Source code here.

We have noticed that there can be correctness issues with our simulate_de implementation when running long time horizon problems, which makes me suspect we need to tweak the solve settings we are using. To that end we are seeking advise on the following points,

  • Are there changes we can make to improve the solution accuracy for long time horizons?
  • Are there changes we can make to improve runtime performance?
  • Are there any ways to avoid the precompile time of DifferentialEquations when running in CI? Compilation of DifferentialEquations is currently about 80% of our CI runtime.

Here is a MWEs to highlight the correctness issue. Note that sum(z_measure_probabilities(ρ) should always be 1.0 in a perfect simulation.

using QuantumAnnealing

ising_model = Dict((1,) => 1.0, (1,2) => -1.0)

ρ = simulate_de(ising_model, 1.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 7.214235431263205e-9

ρ = simulate_de(ising_model, 10.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 8.508320654354584e-6

ρ = simulate_de(ising_model, 100.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 0.0003582731716551546

ρ = simulate_de(ising_model, 1000.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 0.003550520566418469

ρ = simulate_de(ising_model, 10000.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 0.03483614681737934

ρ = simulate_de(ising_model, 100000.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 0.2990146971495551


ρ = simulate(ising_model, 10000.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 1.3211653993039363e-14

ρ = simulate(ising_model, 100000.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 1.0436096431476471e-13


ρ_target = simulate(ising_model, 10000.0, AS_LINEAR); ρ_de = simulate_de(ising_model, 10000.0, AS_LINEAR); error = sum(abs.(ρ_target .- ρ_de)) # 0.0348400079153087

ρ_target = simulate(ising_model, 100000.0, AS_LINEAR); ρ_de = simulate_de(ising_model, 100000.0, AS_LINEAR); error = sum(abs.(ρ_target .- ρ_de)) # 0.2990175297085176

CC @zmorrell

2

I see from the source code that you are just calling solve(prob), thus letting DifferentialEquations automagically choose the solver for you.
While that works out well for most simple cases, I think you might want to have more control on that.

I didn’t look into the details of the problem here but I believe, in general, a sensible first step when you have an unstable solution is to just try different solvers and timesteps (adaptive vs fixed), depending on your problem requirements. There’s a lot to choose from: ODE Solvers · DifferentialEquations.jl. It can really make a difference.

Also, a quick plot of the solution vs time rather than just the endpoint might give you some more insight (e.g., is it just straight blowing up, doing crazy oscillations, showing “non-physical” values…?)

3

I assume you’re tweaking tolerance?

Is your f function non-allocating and have you optimized solver choice? Is it stiff? Have you chosen linear solvers etc. effectively?

4

@ccoffrin I tried lowering the relative tolerance while specifying a non-stiff method and that seems to be giving adequate results for the two qubit test model shown here as well as a 3 qubit test case I ran as well. I suggest that we use Vern9() and fiddle with the default tolerances in the simulate_de function.

2 Likes
5

Thanks all for out input on this. We were able to get the solve method to produce correct in all of the cases I listed above by adjusting the abstol and reltol. This was a great help, thanks!