Issue with specific coupling strength for time-dependent dynamics.

[Sion2002-create]

Hi,

I am having some issues for the coupling strength of alpha = 0.01 (well actually 0.01 * 2.2^2) when trying to implement the time-dependent dynamics using OQuPy, I have attached the following code below:

import numpy as np
import oqupy
from oqupy import operators as op
import time

omega_cutoff = 2.2
temperature = 4 * 0.1309
dt = 0.01
max_time = 5.0

alphas_to_test = [0.02, 0.03, 0.01]

for alpha in alphas_to_test:
    alpha_scaled = alpha * omega_cutoff**2
   
    correlations = oqupy.PowerLawSD(
        alpha=alpha_scaled, zeta=3, cutoff=omega_cutoff,
        cutoff_type='gaussian', temperature=temperature)
    bath = oqupy.Bath(oqupy.operators.sigma("z"), correlations)
   
    tempo_params = oqupy.TempoParameters(dt=dt, tcut=5.0, epsrel=1e-10)
   
    print(f"--- Testing alpha = {alpha} ---")
    start = time.time()
   
    # The problematic core computation
    process_tensor = oqupy.pt_tempo_compute(
        bath=bath,
        start_time=0.0,
        end_time=max_time,
        parameters=tempo_params
    )
   
    end = time.time()
    print(f"Result: alpha={alpha} took {end-start:.2f} seconds.\n")

In this case I found the following outcome after running the code:

Result: alpha=0.02 took 966.12 seconds.
Result: alpha=0.03 took 1200.45 seconds.
Result: alpha=0.01 took 50548.38 seconds.

For whatever reason for higher coupling strengths I have no issues running the code. However, for the particular value of alpha = 0.01, the code becomes very slow which I wouldn't expect compared to higher values of alpha. I would appreciate any help you could give me with this.

[gefux]

Hi @Sion2002-create,
Thank you for reaching out to us on this. Before responding, I'd like to first point out something physical that might confuse this conversation if not clarified first: One might naively expect that for α → 0 the Born-Markov approximations are valid and thus the process is Markovian, reflected by a process tensor of low bond dimension. While indeed the master equation (i.e. the evolution of the density matrix) in this case indeed approaches the GKSL equations, the full process (including all multi-time correlations) is not fully Markovian, even in the α → 0 limit. Thus there is in general no expectation that the process tensor would factorize, or become easy to compute.
Having said that, the general observation I, personally, have made so far is that lower coupling resulted in quicker computations and lower process tensor bond dimensions. (@piperfw: Have you ever observed something contrary, as described above?) While I don't think we can exclude a physical reason for your observation (i.e. an inherent strong non-Markovianity of the process at low coupling and finite temperature), I can also imagine that this might be some numerical issue. To that end, I wonder about the following questions:

1. What are the bond dimensions of the resulting process tensors?
2. Do the corresponding system dynamics appear physically plausible?
3. Do you find the same tendency of longer computation times for smaller alpha also at lower epsrel, say 1.0e-7, 1.0e-8, 1.0e-9?

[piperfw]

Hi both,
Thanks @Sion2002-create for raising this issue. I'm away touring without access to a computer so can't easily check the codebase, but would be interested to look into this further once I'm back late next week if still unresolved then.

My initial thoughts are this does seem pathological, because calculation of the correlation integrals should in fact be independent of alpha (@gefux please correct me if wrong), so I wonder if this is some numerical issue e.g. when performing svd sweeps at that epsrel. I think gefux's suggestions are good, although we can probably determine what's going wrong without system dynamics if you have not set up / have physical expectations for that part of the calculation yet.

To sanity check my above claim, could you try the plot_correlations_with_parameters and guess_tempo_parameters helper functions on your bath for different values of alpha? (You can see example usage of both functions in the Parameters Tutorial). I am expecting the output to be the same in each case.

Then, as a next step I would actually try repeating the three calculations with the values suggested by guess_parameters, which will likely be with a lower esprel (as in gefux's suggestions 3.), to see if the issue is indeed tied to precision or other comoutational parameter somehow.

[paulreastham]

Hi @Sion2002-create @piperfw @gefux, I had a look at this and found I can reproduce it using the TEBD code we have under development as it's quicker to run and so test. I played with various parameters and it seems to be to do with the interplay of dt and alpha. I used a much larger epsrel (1e-7) and find the same slowing down for dt=0.01 and alpha=0.01, but not for dt=0.005 or dt=0.02 with that alpha. It's also present for a different cut-off, and z=2, but not z=1 (so to do with long-time behaviour of correlation function). I think all this supports @piperfw's suggestion that this is a numerical pathology, perhaps an unfortunate cancellation of terms. It's probably hard to track down beyond that but perhaps doesn't need to be, as one could modify the parameters to avoid it perhaps?

[Sion2002-create]

Hi I must apologise for the slow response. I was away last week, I have managed to check the first and last question that @gefux suggested, I haven't got around to @piperfw's suggestion yet. The maximum bond dimension in each case found is summarised below:

Alpha = 0.01

epsrel = 1e-06 | max χ = 6
epsrel = 1e-07 | max χ = 22
epsrel = 1e-08 | max χ = 305
epsrel = 1e-09 | max χ = 338
epsrel = 1e-10 | max χ = 364
epsrel = 1e-11 | max χ = 916

Alpha = 0.02

epsrel = 1e-06 | max χ = 10
epsrel = 1e-07 | max χ = 98
epsrel = 1e-08 | max χ = 33
epsrel = 1e-09 | max χ = 41
epsrel = 1e-10 | max χ = 51
epsrel = 1e-11 | max χ = 228

Alpha = 0.03

epsrel = 1e-06 | max χ = 17
epsrel = 1e-07 | max χ = 250
epsrel = 1e-08 | max χ = 35
epsrel = 1e-09 | max χ = 37
epsrel = 1e-10 | max χ = 98
epsrel = 1e-11 | max χ = 310

Here χ is the bond dimension, it seems to blow up for lower SVD tolerances for alpha = 0.01, this is most likely responsible for the much longer runtimes. I also find the values for 1e-7 a little strange for the higher coupling strengths. Here are the runtimes for 1e-7, 1e-8 and 1e-9 SVD tolerances for building the process tensor

epsrel = 1.0e-9

Result: alpha=0.01 took 11338.90 seconds.
Result: alpha=0.02 took 524.21 seconds.
Result: alpha=0.03 took 675.28 seconds.

epsrel = 1.0e-8

Result: alpha=0.01 took 7526.12 seconds
Result: alpha=0.02 took 270.67 seconds.
Result: alpha=0.03 took 311.83 seconds

epsrel = 1.0e-7

Result: alpha=0.01 took 3800.68 seconds.
Result: alpha=0.02 took 7119.64 seconds.
Result: alpha=0.03 took 6783.17 seconds.

The trend continues at 1.0e-9 and 1.0e-8 for the slowing down at alpha = 0.01. However, 1.0e-7 has quite a different result I'm not sure what causes this as @paulreastham I'm guessing didn't find this, maybe I made some mistake. To me it seems like quite a difficult problem to track down I might try @piperfw's suggestion and looking more at comparing to another implementation of TEMPO I'm using namely uniTEMPO that does not seem to have this issue. Thanks for all the help.

[emile-cochin]

Hi,

I tried running the attached code with the modified numerical integration scheme for the correlation function of #150/#151 and found that the smaller alpha value took less time than the others as expected (1367s for 0.03, 889s for 0.02 and 566s for 0.01), unlike when I tried the current integrator.
My guess is that the current numerical integration of the correlation function has issues handling the small integrand for the smallest alpha value, resulting in numerical error in the $\eta$ coefficients that are larger than they should, in turn resulting in a longer computation time.

I will deal with the final changes of this PR in the next 2 weeks, but in the meantime the alternative integrator is available in this branch (using it requires setting alt_integrator=True in oqupy.PowerLawSD(..., alt_integrator=True)).

[Sion2002-create]

Hi,

I believe that @emile-cochin's suggestion was a good one as the code now seems to run as expected, the bond dimension seems more reasonable. I got the following bond dimensions:

Alpha = 0.01

epsrel = 1e-06 | max χ = 6
epsrel = 1e-07 | max χ = 22
epsrel = 1e-08 | max χ = 22
epsrel = 1e-09 | max χ = 30
epsrel = 1e-10 | max χ = 41

And the time taken seems more reasonable:

Result: alpha=0.01 took 609.15 seconds
Result: alpha=0.02 took 973.00 seconds
Result: alpha=0.03 took 1291.47 seconds.

These are for epsrel = 1e-10, I think the issue may now be solved?