Tip

Run this Jupyter Notebook locally: Download Notebook

Performance of evolution

Cutoff in the coefficients

The time evolution operator for a Hamiltonian \(H\)

\[O[H, \Delta t] = \text{exp}\left(-i H \Delta t\right)\]

uses a KPM expansion as a function of time \(\Delta t\), instead of the usual Fermi energy \(\varepsilon\) of spectral functions.

In the expansion

\[O[H, \Delta t] = \sum_{m=0}^{\infty} T_m(H) c_m(\Delta t)\]

the coefficients \(c_m\) are proportional to the Bessel functions of the first kind of order \(m\), so they decay exponentially fast after a certain threshold \(m^*\).

Thus, is it extremely efficient to cutoff the expansion for a certain target accuracy

\[O[H, \Delta t] \simeq \sum_{m=0}^{m^*} T_m(H) c_m(\Delta t)\]

The exact value of \(m^*\) depends on the time step times the bounds of the spectrum \(\Delta t\times\Delta E\), and on the desired accuracy. To find \(m^*\) for all values of \(\Delta t\times\Delta E\), we define a few typical target accuracies: \(10^{-8}, 10^{-16}, 10^{-32}\), and numerically obtain the answer. Then, we will use a numerical estimate of the function that fits \(m^*\).

import numpy as np
import holoviews as hv
hv.extension()

import scipy
import scipy.optimize
from scipy.special import jv

from kpm_tools.evolution import coef_evo

numerical evaluation of the coefficients

First, we obtain an ndarray of the coefficients vs the time delta \(\Delta t\) and the number of moments num_moments.

Here, we use \(\Delta E = 1\), such that dt is the rescaled time delta \(\Delta t \times \Delta E\).

dt_array = np.logspace(-3,2,1000)
num_moments = 500

coefs_array = np.array([
    coef_evo(num_moments, 0, ab=(2,0), delta_time=dt)
    for dt in dt_array
])
coefs_array = coefs_array[:,1:,0]

num_moments_array = np.arange(1, num_moments)

make a fit

Numerical estimation of the best num_moments

Then we obtain the best num_moments for a given dt and accuracy.

accuracy = 1e-16

# get the transition point in `m` such that coefs are smaller than accuracy

def max_m_array(coefs_array, accuracy):
    return np.count_nonzero(np.abs(coefs_array) > accuracy, axis=1)

(
    hv.QuadMesh((
        num_moments_array,
        dt_array,
        np.log(np.abs(coefs_array))/np.log(10)
    ),
        kdims=['M', '$dt$'],
        vdims=hv.Dimension('z', label='${log}_{10}|c_n|$', range=(np.log(accuracy)/np.log(10), 0))
    ).opts(colorbar=True, fig_size=200)
    *
    hv.Curve((
        max_m_array(coefs_array, accuracy),
        dt_array,
    ))
)[:30,:4]

/tmp/ipykernel_5521/785646435.py:5: RuntimeWarning: divide by zero encountered in log
  np.log(np.abs(coefs_array))/np.log(10)

We can already see that only a handful of moments is needed for a fairly large vale of the time step dt.

Analitical fit of the best num_moments

Since we don’t want to store arrays of best num_moments, we make an analytical approximation to the best num_moments for a given accuray and as a function of dt.

def func(dt, A, B, C, D):
    return (A * np.log(1/accuracy) * (dt) **  B) + C * dt + D * np.log(1/accuracy)

p_acc16, _ = scipy.optimize.curve_fit(func,
                         dt_array,
                         max_m_array(coefs_array, accuracy),
                         p0=(0.5, 0.25, 2, 0),

#                                  bounds=([0,0,0,2], [1,1,1,2+1e-8])
                        )

hv.QuadMesh((
    num_moments_array,
    dt_array,
    np.log(np.abs(coefs_array+np.finfo(float).eps))/np.log(10)
),
    kdims=['M', '$dt$'],
    vdims=hv.Dimension('z', label='$|c_n|$', range=(np.log(accuracy)/np.log(10), 0))
).opts(logy=True, logz=False, colorbar=True, fig_size=200) * \
hv.Curve((np.array([func(dt, *p_acc16) for dt in dt_array]),
          dt_array)).opts(logy=True, logx=True)

Result

Finally, these parameters give the best choice of num_moments as a function of dt for typical accuracy goals: 1e-8, 1e-16, 1e-32.

'p_acc8 = [0.42464678, 0.33500138, 2.00044861, 0.07265139]'
'p_acc16 = [0.34310197, 0.33722835, 2.01232381, 0.08868359]'
'p_acc32 = [0.28132266, 0.32958046, 2.04428473, 0.09146841]'

'p_acc32 = [0.28132266, 0.32958046, 2.04428473, 0.09146841]'

From the scaling analysis we can see that large time steps are possible, and require a proportionally larger amount of moments in the KPM expansion.

Smaller timesteps have a slight overhead in the number of moments: \(m^*=4\) for very small time steps of \(10^{-3}\) but only \(m^*=6\) for a ten times larger time step \(10^{-2}\).

Thus, it is more convenient -in terms of performance- to make one big time step, rather than many subsequent smaller time steps.

Nevertheless, it is still fairly efficient to make many time steps for visualizing the time evolution.