Technical Questions on Performing Imaginary-Time Evolution Using TDVP for Extended Hubbard Model

DMRG and Numerical Methods
,
1

I am currently solving the ground state of the “extended” Hubbard model, where the attractive potential is induced by electron-phonon interaction and is dynamically generated in another process synchronized with the imaginary-time evolution in TDVP. I have a few technical questions regarding this:

  1. Imaginary-Time Evolution Operator: Strictly speaking, the imaginary-time evolution should be e^{-(H-\langle H\rangle)\tau}, but in Example01 under ITensorTDVP.jl, the evolution operator is just H. However, when the spectrum is sufficiently negative (e.g. in a large system), KrylovKit seems to report an error: ERROR: LoadError: operator does not appear to be hermitian: 0.0 vs NaN. Subtracting the ground state energy can somewhat avoid this issue. However, how should we subtract \langle H\rangle? My current approach is:
E_electron = inner(psi', H, psi)
Hτ = let os = OpSum()
    os += (-E_electron, "Id", 1)
    O = MPO(os, siteinds(psi))
    H + O
end

Based on my understanding of TDVP, which essentially involves splitting the tangent space projector, the generation of Hτ should be executed on each bond. However, this might be unnecessary as the energy generally doesn’t change significantly, and the operation H + O is expensive. Additionally, for the H + O operation, can we use alg = "directsum"? What would be the consequences of this?
2. Optimal Constant for Subtraction in Imaginary-Time Evolution: Continuing from the previous question, although the strict derivation of the imaginary-time evolution equation indicates that the Hamiltonian’s expectation value should be subtracted (essentially normalize the state), as long as the relative order of the spectrum doesn’t change, I can still monotonically find the ground state (since normalization is performed after each sweep). Is there a rule of thumb regarding how much of a constant should be subtracted? Or is subtracting the expectation value the most reliable method? In other words, how can we efficiently approximate e^{-A\tau} \cdot v in Krylov space? I am not very familiar with this aspect, but I have heard that the speed of convergence is determined by the norm of the matrix (operator), condition number, etc. If anyone has experience or can share useful references, I would greatly appreciate it!
3. Reverse Step in TDVP: From my observation in Example01 under ITensorTDVP.jl, reverse_step=false is set in TDVP. I believe this is understandable because the ground state problem is essentially about finding a fixed point, and this might be faster. However, what could be the possible consequences of this approach? For example, could the energy fail to monotonically decrease (even if the bond dimension is sufficient)? I apologize for my lack of familiarity with this topic, and would greatly appreciate any guidance!

2

Update about 1:
I have tried three different approaches:

  1. H + O: adding MPO with default truncation;
  2. add(H, O, alg="directsum"): adding MPO without truncation;
  3. Pass [H, O] as MPOs into tdvp where the ProjMPOSum will be called internally.

Some representative outputs are:

1. # adding MPO with truncation
E_electron = -52.13917886125441
linkdims(Hτ) = [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 37, 38, 39, 39, 38, 38, 37, 38, 38, 38, 38, 38, 39, 39, 39, 38, 37, 37, 36, 35, 34, 33, 32, 31, 30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7]
linkdims(psi) = [1, 1, 1, 1, 1, 3, 3, 4, 4, 12, 13, 16, 36, 47, 60, 97, 106, 142, 176, 215, 215, 284, 244, 355, 270, 397, 354, 382, 388, 371, 410, 282, 317, 249, 240, 187, 156, 112, 119, 145, 113, 115, 156, 160, 166, 225, 230, 234, 343, 338, 361, 390, 315, 322, 308, 315, 284, 271, 255, 194, 212, 169, 109, 92, 58, 48, 39, 16, 13, 11, 4, 4, 3, 3, 1, 1, 1, 1, 1]
After sweep 1: maxlinkdim=397 maxerr=9.82E-07 current_time=0.06 time=411.757
----------------------------------------------------------------------
E_electron = -54.80946995491526
linkdims(Hτ) = [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 37, 37, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38, 39, 39, 39, 38, 37, 37, 36, 35, 34, 33, 32, 31, 30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7]
linkdims(psi) = [1, 1, 1, 1, 1, 3, 3, 4, 4, 12, 13, 16, 36, 47, 60, 95, 104, 140, 172, 209, 210, 276, 240, 344, 262, 383, 346, 369, 377, 365, 397, 274, 309, 240, 234, 183, 152, 112, 116, 141, 110, 113, 152, 156, 163, 219, 226, 227, 333, 328, 348, 380, 307, 316, 300, 307, 275, 262, 249, 190, 204, 164, 107, 90, 58, 48, 39, 16, 13, 11, 4, 4, 3, 3, 1, 1, 1, 1, 1]
After sweep 1: maxlinkdim=385 maxerr=9.91E-07 current_time=0.07 time=382.639

2. # directsum MPO 
E_electron = -53.06328573894331
linkdims(Hτ) = [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 47, 46, 45, 44, 43, 42, 41, 40, 39, 38, 37, 36, 35, 34, 33, 32, 31, 30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9]
linkdims(psi) = [1, 1, 1, 1, 1, 3, 3, 4, 4, 12, 13, 16, 36, 47, 60, 92, 97, 130, 154, 159, 202, 266, 235, 324, 258, 364, 275, 368, 342, 344, 330, 244, 302, 208, 213, 118, 133, 103, 87, 108, 96, 99, 130, 114, 134, 156, 210, 245, 333, 331, 355, 364, 345, 370, 270, 339, 249, 281, 220, 220, 175, 140, 101, 96, 57, 47, 35, 16, 12, 10, 4, 4, 3, 3, 1, 1, 1, 1, 1]
After sweep 1: maxlinkdim=358 maxerr=9.93E-07 current_time=0.06 time=711.231
----------------------------------------------------------------------
E_electron = -55.788064558771566
linkdims(Hτ) = [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 47, 46, 45, 44, 43, 42, 41, 40, 39, 38, 37, 36, 35, 34, 33, 32, 31, 30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9]
linkdims(psi) = [1, 1, 1, 1, 1, 3, 3, 4, 4, 12, 13, 16, 36, 47, 60, 90, 95, 126, 150, 153, 197, 258, 229, 313, 251, 351, 268, 358, 328, 334, 316, 236, 294, 206, 204, 116, 129, 101, 87, 104, 94, 97, 127, 112, 132, 154, 206, 239, 320, 323, 345, 353, 336, 357, 265, 327, 243, 273, 214, 214, 170, 138, 99, 94, 57, 47, 35, 16, 12, 10, 4, 4, 3, 3, 1, 1, 1, 1, 1]
After sweep 1: maxlinkdim=348 maxerr=9.99E-07 current_time=0.07 time=692.349

3. # tdvp with ProjSumMPO
E_electron = -52.457098556393056
linkdims(H) = [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 45, 45, 44, 43, 42, 41, 40, 39, 38, 37, 36, 35, 34, 33, 32, 31, 30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7]
linkdims(psi) = [1, 1, 1, 1, 1, 3, 3, 4, 4, 12, 13, 16, 36, 47, 56, 80, 91, 132, 155, 195, 214, 268, 238, 323, 268, 360, 290, 364, 305, 357, 294, 263, 298, 257, 230, 145, 132, 112, 101, 100, 110, 85, 100, 111, 140, 164, 166, 185, 233, 210, 240, 260, 239, 281, 238, 275, 234, 236, 218, 180, 185, 149, 106, 94, 58, 47, 35, 16, 13, 10, 4, 4, 3, 3, 1, 1, 1, 1, 1]
linkdims(O) = [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
After sweep 1: maxlinkdim=353 maxerr=9.99E-07 current_time=0.06 time=283.385
----------------------------------------------------------------------
E_electron = -55.119375825532586
linkdims(H) = [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 45, 45, 44, 43, 42, 41, 40, 39, 38, 37, 36, 35, 34, 33, 32, 31, 30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7]
linkdims(psi) = [1, 1, 1, 1, 1, 3, 3, 4, 4, 12, 13, 16, 36, 47, 54, 78, 89, 129, 149, 190, 208, 259, 232, 314, 262, 347, 282, 353, 299, 349, 284, 257, 292, 249, 226, 143, 128, 112, 99, 100, 108, 83, 100, 109, 138, 160, 164, 181, 227, 205, 235, 253, 234, 274, 233, 268, 228, 230, 214, 175, 180, 144, 104, 92, 58, 47, 35, 16, 13, 10, 4, 4, 3, 3, 1, 1, 1, 1, 1]
linkdims(O) = [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
After sweep 1: maxlinkdim=344 maxerr=9.94E-07 current_time=0.07 time=474.753

Conclusions:

  1. Subtracting \langle H \rangle from H with default truncation can lead to MPO with smaller bond dimensions, which makes tdvp faster and the truncation overhead is negligible. But this is elusive to me, if someone can give a hint, I’d greatly appreciate it!
  2. The unusual fast performance of 3(time=283.385) may results from the fast convergence of KrylovKit, which is off-topic but need to be studied further by extracting the convergenceInfo of krylovkit.

BTW, the Hamiltonian is dynamically generated and dependent of another synchronized process dealing with phonon state. The whole process is essentially the iteration between tdvp and phonon-evolution.

3

More math on subtracting \langle H\rangle in imaginary-time evolution:

|\varphi(\tau)\rangle=\frac{e^{-H \tau}|\varphi(0)\rangle}{\sqrt{\left\langle\varphi(0)\left|e^{-2 H \tau}\right| \varphi(0)\right\rangle}} \Rightarrow \partial_\tau|\varphi(\tau)\rangle=-(H-\langle H\rangle)|\varphi(\tau)\rangle

But there’s more to consider. For a non-normalized state |\varphi(\tau)\rangle:

E(\tau) = \frac{\langle \varphi(\tau) | H | \varphi(\tau) \rangle}{\langle \varphi(\tau)|\varphi(\tau)\rangle} = \frac{\sum_{n}E_{n}(\tau)e^{-2(E_{n}(\tau) - E(0))\tau}|c_n|^2}{\sum_{n}e^{-2(E_{n}(\tau) - E(0))\tau}|c_n|^2} = -\frac{1}{2}\partial_{\tau}\ln\langle \varphi(\tau)|\varphi(\tau)\rangle + E(0)

Thus, \langle H\rangle only needs to be calculated once and subtract half of the lognorm of MPS during all the necessary normalization process.

Regarding question 3, I still cannot figure it out. I know Jutho mentioned the equivalence of DMRG and TDVP in his 2016 paper, and DMRG only contains forward step. However, this is because the forward step has already reached the eigenstate. As for finite-step TDVP, I’m not sure the potential consequences of skipping reverse step. It seems like evolving each sites (except the boundary) 2τ and eventually we will arrived at ground state.