Learning Feynman Diagrams with Tensor Trains

We use tensor network techniques to obtain high-order perturbative diagrammatic expansions for the quantum many-body problem at very high precision. The approach is based on a tensor train parsimonious representation of the sum of all Feynman diagrams, obtained in a controlled and accurate way with the tensor cross interpolation algorithm. It yields the full time evolution of physical quantities in the presence of any arbitrary time-dependent interaction. Our benchmarks on the Anderson quantum impurity problem, within the real-time nonequilibrium Schwinger-Keldysh formalism, demonstrate that this technique supersedes diagrammatic quantum Monte Carlo by orders of magnitude in precision and speed, with convergence rates $1/N^2$ or faster, where $N$ is the number of function evaluations. The method also works in parameter regimes characterized by strongly oscillatory integrals in high dimension, which suffer from a catastrophic sign problem in quantum Monte Carlo calculations. Finally, we also present two exploratory studies showing that the technique generalizes to more complex situations: a double quantum dot and a single impurity embedded in a two-dimensional lattice.

Popular Summary

Many problems in physics have a formal solution in terms of high-dimensional integrals or sums. Examples include partition functions in statistical physics, path integrals, or Feynman diagrams in quantum field theory and many-body physics. To evaluate these integrals numerically, there is a unique class of technique: Monte Carlo sampling. While Monte Carlo techniques have been extremely successful in solving many problems, they also fail spectacularly when the integrand oscillates strongly, a common situation for fermionic quantum many-body problems. Here, we propose an alternative approach based on a tensor network representation of the integrand.

Our technique “learns” the integrand and constructs a compressed representation from which the high-dimensional integrals can be obtained easily. As in image compression, which leverages the existence of regularities in real-life images, the success of our approach is linked to the existence of an underlying structure in the theory, “factorizability,” which the algorithm unveils. The approach is fully general and, importantly, the factorizability property is orthogonal to what is needed for the success of Monte Carlo—the absence of a “sign problem” or, for variational matrix-product-state approaches, a low entanglement.

We showcase our approach by calculating Feynman diagrams up to very high order in the context of three quantum impurity models. Our results indicate an unprecedentedly fast convergence, including in difficult regimes where a catastrophic sign problem would forbid the calculation with Monte Carlo.

Figure 1

Overview of the main results. (a) Slice of the corresponding integrand ${\tilde{Q}}_{10}$ (orange line) compared to the MPS approximation (blue dots). The values of $u_2,u_3,\dots ,u_{10}$ are arbitrarily fixed (vertical dashed lines). (b) Relative error of the $n$th coefficient (for $n=19$) in the perturbative expansion of the charge $Q$ of the Anderson quantum impurity model (compared with the exact Bethe ansatz solution) versus the number $N$ of evaluation of ${\tilde{Q}}_{19}$. (c) Average sign defined as ${\eta }_n=Q_n/\int |{\tilde{Q}}_n|$ for $n=10$ versus on-site energy ${\epsilon }_d$. (d) Relative error of the rank-50 MPS approximation [pivot error, as defined in Eq. (24), divided by the value of the function ${\tilde{Q}}_n$ at the first pivot] versus ${\epsilon }_d$.

Figure 2

Illustration of the CI of a matrix. The large red triangles indicate real pivots, and the smaller red triangles indicate automatically generated pivots. The right-hand side contains only small subparts of the matrix.

Figure 3

Step by step representation of a simple algorithm to factorize a multidimensional tensor into a tensor train. The blue squares represent the inverses of the pivot matrices. Summation is implicit over the indices connecting two tensors (green lines).

Figure 4

Pictorial representation of TCI formula for $A$ defined in Eq. (16). The $T_{\alpha }$ tensors and $P_{\alpha }^{-1}$ pivot matrices are represented by orange and blue squares, respectively. The green lines correspond to contracted discrete indices ${\mathcal{I}}_{\alpha }$ or ${\mathcal{J}}_{\alpha }$, as indicated by the arrows. The thin black lines correspond to the variables $u_{\alpha }$.

Figure 5

Pictorial representation of the ${\mathrm{\Pi }}_{\alpha }$ tensor and its cross interpolation. The notation is the same as in Fig. 4.

Figure 6

Real-time noninteracting lesser Green’s functions for the three models considered in this paper. Upper panel: single quantum dot (SIAM) and double quantum dot connected to leads in the flat-band limit. Lower panel: one interacting site in the 2D lattice, including the tail of $g^{<}(t)$.

Figure 7

(a) Upper: factorization error estimator as a function of the tensor rank $\chi$, for the SIAM model, for $n=10$, ${\epsilon }_d=0$, $t_M=15$, Gauss-Kronrod points with 63 points, using the error functions ${\epsilon }_{\mathrm{\Pi }}^{\mathrm{env}}$ (orange curve) and ${\epsilon }_{\mathrm{\Pi }}$ (blue curve); see the text. (b) Lower: error of the integral $Q_n$ versus $\chi$, measured using different estimators. In this case, $Q_{10}^{\mathrm{last}}$ is $Q_{10}$ for $\chi =100$. The red line is the error due to the integral discretization obtained by varying the number of integration points.

Figure 8

Comparison of the integrand ${\tilde{Q}}_7$ for ${\epsilon }_d=-2$ (red stars) with the TCI of rank $\chi =2$, 4, 40 (dash-dotted, dashed, and solid black lines, respectively). The function is plotted in the $w$ variable (56) with $t_M=5$. The maximum error point is (0.71,0.62,0.0,0.07,0.45,0.60,0.25) (blue circles). Each panel corresponds to the variation of one variable $w_i$ starting from this point.

Figure 9

Relative error of the coefficients $Q_n$ with respect to the exact solution for $n=3$ (blue), $n=5$ (orange), $n=10$ (green), $n=15$ (red), and $n=19$ (purple). All calculations are carried out using the env and alternate search variants of the algorithm. (a) Error versus tensor rank $\chi$. Large values of $t_M$ and $d$ are used so that the accuracy is limited only by $\chi$: $t_M=20$, ${\mathrm{CH}}_{255}$ ($n=3$), $t_M=25$, ${\mathrm{CH}}_{255}$ ($n=5$), and $t_M=30$, ${\mathrm{GK}}_{63}$ ($n=10$, $n=15$, $n=19$) The dashed lines are guides to the eye for $1/{\chi }^p$ scaling, which corresponds to $1/N^{p/2}$ in terms of the number $N$ of function calls. $p=1$ corresponds to Monte Carlo scaling. (b) Error versus $N$ in a low-precision calculation with $t_M=15$ and $d=15$ (${\mathrm{GK}}_{15}$). The colored circles indicate the CPU time of the calculation when performed on a single core. (c) Error versus $n$ for $\chi =10$ (blue), $\chi =20$ (orange), and $\chi =30$ (green) ($t_M=30$ and ${\mathrm{GK}}_{15}$). (d) $|Q_n|$ versus $n$, distinguishing positive and negative coefficients. The vertical lines indicate values of $n$ where $|Q_n|$ is very small, leading to larger relative errors.

Figure 10

$Q_n(t)$ after a quench at $t=0$ for the SIAM model at ${\epsilon }_d=\overline{\alpha }=0$, in equilibrium. The coefficients are normalized by their exact asymptotic values $Q_n(t=\infty )=Q_n^{\text{Bethe}}$. Upper: abrupt quench $\lambda (t)=\theta (t)$. Different curves correspond to different integration techniques, with indistinguishable results: simplex integration (55) (yellow continuous line), Fourier technique (h4) (red dotted line), quantum quasi-Monte-Carlo (QQMC) [33] (blue dashed line), and direct numerical integration of Eq. (42) for $n=2$, 3, 4 (green dashed-dotted line). The arrows indicate the value of $n$. Lower: $Q_n(t)$ for $n=3$ (orange) and $n=10$ (blue) using a continuous increase of the interaction strength $\lambda (t)=\sin (\pi t/2\tau )$ for $t<\tau$, and $\lambda (t)=1$, for $t\ge \tau$ (with simplex integration).

Figure 11

Charge $Q(U,t)$ on the dot for the SIAM model after switching on the interaction abruptly (yellow solid line, $\tau =0$) or smoothly (blue dashed line, $\tau =2$) with $\lambda (t\le \tau )=\sin (\pi t/2\tau )$ and $\lambda (t\ge \tau )=1$. Both curves include terms $Q_n(t)$ up to $n=14$. The black dotted line corresponds to the resummed series using the Euler transform of the first 20 $Q_n^{\text{Bethe}}$ coefficients, as in Ref. [25]. In all calculations, we use ${\epsilon }_d=\overline{\alpha }=0$ with $U=2,t_M=20$.

Figure 12

Integrals $Q_n=\int {\tilde{Q}}_n$ and $\int |{\tilde{Q}}_n|$ for $n=10$ versus ${\epsilon }_d$ for the SIAM model. The green and orange colors correspond $Q_{10}<0$ and $Q_{10}>0$, respectively.

Figure 13

(a) Upper: error of $Q_{10}$ versus $\chi$ in the SIAM model, for different ${\epsilon }_d$. Gauss-Kronrod integration is used with $d=31$ points for ${\epsilon }_d=2,0,-0.3$ and $d=255$ points for ${\epsilon }_d=-2,-4$. (b) Lower: the same error, but for fixed $\chi =50$, varying ${\eta }_{10}$.

Figure 14

Relative error with respect to the Bethe ansatz for $Q_{10}(t_M)={\int }_0^{t_M}{\tilde{Q}}_{10}$ and ${\int }_0^{15}{\tilde{Q}}_{10}^{\text{extra}}$, versus $t_M$. ${\tilde{Q}}_{10}^{\text{extra}}$ is the extrapolation of the TCI of ${\tilde{Q}}_{10}$ obtained for $v_i\le t_M$. We used Gauss-Kronrod rule with 63 points for the integration. The error saturates at ${10}^{-5}$ due to the discretization error.

Figure 15

Convergence with respect to the tensor rank $\chi$ of the ten-dimensional integral $Q_{10}$, for the SIAM model. All panels show the relative error ${\epsilon }_{10}^Q=|Q_{10}(\chi )/Q_{10}^{\text{Bethe}}-1|$ measured against the Bethe ansatz solution versus the bond dimension $\chi$. (a) Comparison of the choices of variables $u$, $v$, and $w$. (b) Comparison of different one-dimensional quadrature rules of $d$ nodes: Gauss Kronrod (${\mathrm{GK}}_d$) and Chebyshev (${\mathrm{CH}}_d$). (c) Comparison of different pivot selection algorithms, full search and alternate search, defined in Sec. 3b3, and the weighted learning variant defined in Sec. 3b5. (d) Effect of using ${\epsilon }_{\mathrm{\Pi }}$ (23) or ${\epsilon }_{\mathrm{\Pi }}^{\mathrm{env}}$ (30) as the error function in the pivot selection. (e) Saturation of the error due to the maximum time cutoff $t_M$. (f) Absence of a dependence on the pivot acceptance condition; see the text.

Figure 16

One interacting site in an infinite two-dimensional lattice. Upper: integrand ${\tilde{Q}}_5(v_1,v_2,v_3,v_4,v_5)$ versus $v_3$ for a random choice of the values $v_1$, $v_2$, $v_4$, and $v_5$. The right panel enlarges the tail. Actual integrand (orange curve) and TTD approximation at $\chi =150$ (blue circles) are both shown. Lower: $|Q_n|$ versus $n$ for three values of the maximum time, $t_M=10$ (stars), $t_M=100$ (pluses), and $t_M=10000$ (circles). Blue (red) symbols correspond to positive (negative) values of $Q_n$.

Figure 17

Three different factorizations used for the double quantum dot problem. From top to bottom: vertex factorization, time factorization, and full factorization.

Figure 18

Double quantum dot. Upper: error of $Q_{10}$ versus $\chi$ for the three algorithms introduced in the text. The reference value $Q_{10}^{\mathrm{ref}}$ is obtained using the time factorization algorithm. Lower: $|Q_n|$ versus $n$, obtained using the vertex factorization algorithm. Blue (red) symbols correspond to positive (negative) values of $Q_n$. We use ${\gamma }_0=0.5\mathrm{\Gamma }$ and $t_M=15$.