High Order Numerical Stability

[ccoffrin]

This seems to work without issue

ρ = simulate_flexible_order(one_spin_model, 100.0, AS_LINEAR, 2, 4)
2×2 Matrix{ComplexF64}:
 0.963726+0.0im         0.177282-0.0594113im
 0.177282+0.0594113im  0.0362744+3.46945e-18im

Around order 12 there appears to be some issues.

ρ = simulate_flexible_order(one_spin_model, 100.0, AS_LINEAR, 2, 13)
2×2 Matrix{ComplexF64}:
 2.08621e11+0.0im         8.62285e10+1.04498e10im
 8.62285e10-1.04498e10im  3.61639e10-9.53674e-7im

ρ = simulate_flexible_order(one_spin_model, 100.0, AS_CIRCULAR, 2, 13)
2×2 Matrix{ComplexF64}:
  0.000163706+6.77626e-21im  -0.000392954+3.23259e-5im
 -0.000392954-3.23259e-5im    0.000949614-5.42101e-20im

 ρ = simulate_flexible_order(one_spin_model, 100.0, AS_LINEAR, 2, 15)
2×2 Matrix{ComplexF64}:
 NaN+NaN*im  NaN+NaN*im
 NaN+NaN*im  NaN+NaN*im

ρ = simulate_flexible_order(one_spin_model, 100.0, AS_CIRCULAR, 2, 15)
2×2 Matrix{ComplexF64}:
 NaN+NaN*im  NaN+NaN*im
 NaN+NaN*im  NaN+NaN*im

[ccoffrin]

The issue is amplified in cases with more qubits.

im = Dict((1,) => -1, (1,2) => -1, (1,3) => -1, (1,4) => -1, (1,5) => -1, (2,3) => 1, (4,5) => 1)
ρ = simulate_flexible_order(im, 100.0, AS_LINEAR, 2, 10)

[ccoffrin]

It seems like long annealing times can also generate an issue,

ising_model = Dict((1,) => 0.1, (1,2) => 1.0, (2,3) => 1.0, (1,3) => 1.0)
ρ = simulate_flexible_order(ising_model, 100000.0, AS_CIRCULAR, 2, 5)

Outcome from in person discussion: this is an unavoidable behavior, requires increasing the number of steps.