Variational optimization algorithms for uniform matrix product states

We combine the density matrix renormalization group (DMRG) with matrix product state tangent space concepts to construct a variational algorithm for finding ground states of one-dimensional quantum lattices in the thermodynamic limit. A careful comparison of this variational uniform matrix product state algorithm (VUMPS) with infinite density matrix renormalization group (IDMRG) and with infinite time evolving block decimation (ITEBD) reveals substantial gains in convergence speed and precision. We also demonstrate that VUMPS works very efficiently for Hamiltonians with long-range interactions and also for the simulation of two-dimensional models on infinite cylinders. The new algorithm can be conveniently implemented as an extension of an already existing DMRG implementation.

Figure 1

Plot of energy density error $\mathrm{\Delta }e$ and gradient norm $\parallel B\parallel$ for VUMPS with a single-site or $N$-site unit cell: (a) TFI model in the gapped symmetry broken phase at $h=0.48$ and $D=25$, (b) gapped isotropic $S=1 XXZ$ antiferromagnet and $D=120$, (c) $S=1/2 XXZ$ antiferromagnet in the gapped symmetry broken phase at $\mathrm{\Delta }=2$ and $D=54$, (d) critical isotropic $S=1/2 XXZ$ antiferromagnet at $\mathrm{\Delta }=1$ and $D=137$, (e) critical Fermi Hubbard model at $U=10$ and $D=256$, (f) critical $S=1/2$ Haldane-Shastry model at $D=200$, and the isotropic $S=1/2$ Heisenberg antiferromagnet on a Cylinder of circumference $W=8$ (g) and $W=10$ (h). The uMPS ground-state approximations of (c), (e), (f), (g), and (h) are not translation invariant and have been obtained from algorithm 3 with a multisite unit cell. Regardless of the criticality of the model, VUMPS converges exponentially fast in gradient norm $\parallel B\parallel$. Notice that at the point where the energy has already converged to machine precision, the gradient is still quite far from zero, and the state thus still some distance from the variational optimum.

Figure 2

(Top) Schmidt spectrum of the $S=1$ Heisenberg antiferromagnet for $D=120$, converged to gradient norm $\parallel B\parallel <{10}^{-15}$. (Bottom) The table shows the first 30 Schmidt values in descending order. The degeneracies are reproduced to 15 digits of precision, without exploiting any symmetries.

Figure 3

Evolution of the Schmidt spectrum (top) and the gradient norm $\parallel B\parallel$ and various observables (bottom) with iteration number for the TFI model at $h=0.45$. Here we defined the deviation of an observable $O$ from its exact value as $\mathrm{\Delta }O=|\langle O\rangle -{\langle O\rangle }_{\mathrm{exact}}|$, similar to (33) for the energy. We used bond dimensions $D=[9,19,33,55]$, i.e., we increased the bond dimension three times during the optimization process as soon as $\parallel B\parallel$ dropped below ${10}^{-14}$ (at iterations 36, 60, and 84). It is apparent that while high lying Schmidt values converge quite quickly, the better part of the final iterations goes into converging low lying Schmidt values. Moreover, one can see that there is quite some rearrangement of exactly these low lying Schmidt values every time $\parallel B\parallel$ reaches a local maximum (e.g., around iterations 12, 40, 70, 88, and 92) during a nonmonotonous phase of gradient evolution (see also Sec. 3b2).

Figure 4

Comparison of the convergence rate of the gradient norm $\parallel B\parallel$ for the TFI model at $h=0.55$, with $D=33$ and 35. Convergence is roughly 4 times faster for $D=33$ as compared to $D=35$. The inset shows the Schmidt spectrum of the ground state (up to $D=43$). For $D=35$ the smallest Schmidt values form an incomplete degenerate multiplet, whereas for $D=33$ the multiplet is complete.

Figure 5

Scaling of the variational ground-state energy $e$ with truncations error ${\epsilon }_{\rho }$ and bond dimension $D$ for the isotropic $S=1/2$ Heisenberg antiferromagnet. We plot the energy $e$ vs truncation error ${\epsilon }_{\rho }$ on the left and vs inverse bond dimension $1/D$ on the right. The exact ground-state energy is given by $e_{\mathrm{exact}}=-0.4431471805599453$ to 16 digits of precision. We obtain estimates from a linear fit $e\left({\epsilon }_{\rho }\right)=e+a{\epsilon }_{\rho }$ (left) and a power law fit $e(1/D)=e+a{(1/D)}^b$ (right), where the estimate on the right has one more digit of precision and is roughly four times more accurate than the estimate on the left.

Figure 6

Example performance comparisons between the sequential and parallel algorithm presented in Sec. 2e. We show examples for the gapped $S=1/2 XXZ$ antiferromagnet at $\mathrm{\Delta }=2$ and $D=87$ with (a) a two-site unit cell, (b) a four-site unit cell, as well as (c) the critical Hubbard model at $U=5$ and $D=66$ with a two-site unit cell and (d) the critical Haldane-Shastry model at $D=60$ with a two-site unit cell. In each example, both approaches were initialized in the same (random) initial state. It can be seen that while the parallel approach usually takes less time for one unit-cell update iteration, the sequential approach shows faster convergence in the number of iterations. Both approaches show similar performance in the considered cases, except for example (d).

Figure 7

Comparative benchmark plots for VUMPS, IDMRG, and ITEBD. We plot the error $\mathrm{\Delta }e$ on the left and the gradient norm $\parallel B\parallel$ on the right vs total computing time $t$ in seconds for (a) the gapped isotropic $S=1 XXZ$ antiferromagnet at $D=120$, (b) the gapped $S=1/2 XXZ$ antiferromagnet at $\mathrm{\Delta }=2$ and $D=54$, (c) the critical isotropic $S=1/2 XXZ$ antiferromagnet at $D=70$, and (d) the critical Hubbard model at $U=5$ and $D=65$. The dashed lines are a guide to the eye and denote the minimum values of $\mathrm{\Delta }e$ and $\parallel B\parallel$ obtained by the respective algorithm. The insets show a plot of the entire ITEBD and/or IDMRG simulation with logarithmic time scale. It is obvious that VUMPS reaches convergence orders of magnitude faster than IDMRG or ITEBD, especially for critical systems. Notice also that, while $\mathrm{\Delta }e$ differs only by a few percent between the different algorithms [up to $\approx 10\%$ for (d)], VUMPS always manages to also converge $\parallel B\parallel$ essentially to machine precision, whereas IDMRG and ITEBD stagnate at some substantially higher values, remaining quite far from the variational optimum.