Polynomial complexity despite the fermionic sign

TL;DR

This paper demonstrates that for certain Feynman diagrammatic Monte Carlo methods, the computational complexity grows only polynomially with inverse error, challenging the notion that the fermionic sign problem leads to exponential difficulty.

Contribution

It shows that specific Monte Carlo approaches to fermionic problems can have polynomial complexity despite the sign problem, offering a new perspective on computational feasibility.

Findings

Computational time increases polynomially with inverse error.

Challenges the belief that the sign problem causes exponential complexity.

Applicable to convergent Feynman diagrammatic series.

Abstract

It is commonly believed that in quantum Monte Carlo approaches to fermionic many- body problems, the infamous sign problem generically implies prohibitively large computational times for obtaining thermodynamic-limit quantities. We point out that for convergent Feynman diagrammatic series evaluated with the Monte Carlo algorithm of [Rossi, arXiv:1612.05184], the computational time increases only polynomially with the inverse error on thermodynamic-limit quantities.

Equations

$$t(ϵ)=O({ϵ}^{-\alpha }).$$

$$⟨A{⟩}_P=\frac{{\sum }_{\mathcal{C}}P(\mathcal{C})A(\mathcal{C})}{{\sum }_{\mathcal{C}}P(\mathcal{C})},$$

$$\frac{1}{N_{MC}}\sum _{i=1}^{N_{MC}}A({\mathcal{C}}_i),$$

$${ϵ}_{\mathrm{stat}}={\sigma }_A\sqrt{\frac{2{\tau }_{\mathrm{auto}}+1}{N_{MC}}},$$

$$a_n=\int dX\sum _{\mathcal{T}\in {\mathcal{S}}_n}\mathcal{D}(\mathcal{T},X),$$

$$a_n=⟨A_n{⟩}_{P_n},$$

$$\left\{\begin{array}[]{ll}\displaystyle P_{n}(\mathcal{T},X)=\left|\mathcal{D}(\mathcal{T},X)\right|&{\rm(DiagMC)}\\ \displaystyle P_{n}(X)=\left|\sum_{\mathcal{T}\in\mathcal{S}_{n}}\mathcal{D}(\mathcal{T},X)\right|&{\rm(CDet)}\end{array}\right.$$

$$A_{n}=\left\{\begin{array}[]{ll}\displaystyle z_{n}~{}{\rm sign}[\mathcal{D}(\mathcal{T},X)]&{\rm(DiagMC)}\\ \displaystyle z_{n}~{}{\rm sign}\left[\sum_{\mathcal{T}\in\mathcal{S}_{n}}\mathcal{D}(\mathcal{T},X)\right]&{\rm(CDet)}\end{array}\right.$$

$$z_{n}=\left\{\begin{array}[]{ll}\displaystyle\int dX\sum_{\mathcal{T}\in\mathcal{S}_{n}}|\mathcal{D}(\mathcal{T},X)|&{\rm(DiagMC)}\\ \displaystyle\int dX\left|\sum_{\mathcal{T}\in\mathcal{S}_{n}}\mathcal{D}(\mathcal{T},X)\right|&{\rm(CDet)}\;.\end{array}\right.$$

$$|a_n|\underset{n\to \infty }{=}O(R^{-n})$$

$$z_n$$

$${\sigma }_{A_n}=z_n\sqrt{1-⟨\mathrm{sign}{⟩}_{P_n}^2}\underset{n\to \infty }{\sim }{z}_n,$$

$${ϵ}_{\mathrm{stat}}(n)\underset{n\to \infty }{\sim }{z}_n\sqrt{\frac{2{\tau }_{\mathrm{auto}}(n)+1}{N_{MC}(n)}}.$$

$${ϵ}_{\mathrm{sys}}(N)=\sum _{n=N+1}^{\infty }{a}_n,$$

$${ϵ}_{\mathrm{sys}}(N)\underset{N\to \infty }{=}O(R^{-N}).$$

$$t_{n}\underset{n\to\infty}{\sim}\left\{\begin{array}[]{ll}\displaystyle\frac{1}{\epsilon^{2}}\,\frac{(n!)^{2}}{R_{D}^{2n}}&{\rm(DiagMC)}\\ \displaystyle\frac{1}{\epsilon^{2}}\,\left(\frac{3}{R_{C}^{2}}\right)^{n}&{\rm(CDet)}\;,\end{array}\right.$$

$$t(ϵ)$$

$$\alpha =2+\frac{\ln (3/{R_C}^2)}{\ln R}$$

$${ϵ}_{\mathrm{stat}}(L)\underset{L\to \infty }{\sim }\frac{e^{\#\beta L^d}}{\sqrt{t}},$$

$${ϵ}_{\mathrm{sys}}(L)\underset{L\to \infty }{\sim }{e}^{-\#L},$$

$$t(ϵ)\sim e^{\#\beta (\ln {ϵ}^{-1}{)}^d}.$$

Full text

Quantum Monte Carlo methods Lattice fermion models (Hubbard model, etc.)

The notion of fermion sign problem (FSP) was originally formulated in the context of auxiliary-field, path-integral and diffusion quantum Monte Carlo (QMC) methods [1, 2, 3, 4]. There, it was observed that the computational time required for calculating properties of the fermionic system to a given accuracy scales exponentially with the system volume. Later, the notion of FSP was implicitly extended to an arbitrary QMC approach dealing with interacting fermions, referring to the stochastic sampling of a non-sign-definite quantity with a near cancellation between positive and negative contributions. Sign-free Monte Carlo (MC) algorithms were emerging only as exceptions confirming the rule: In each such case, the absence of FSP was due to some special property of the simulated model (see, e.g., [5, 6, 7, 8, 9, 10, 11] and Refs. therein). Nowadays the FSP is generally perceived as one of the most important unsolved problems in the field of numerical studies of interacting fermionic systems in dimensions $d>1$.111Frustrated spin and frustrated bosonic (with restricted on-site Hilbert space) lattice models can be mapped to a system of interacting fermions and thus are part of the present discussion.

The main message of this letter follows from a simple observation. Suppose some quantity $Q$ is computed from a limit $Q=\lim_{n\to\infty}Q_{n}$, with an exponentially fast convergence, $|Q-Q_{n}|\sim e^{-\#n}$ (where $\#$ denotes some positive constant). Then in order to compute $Q$ up to an error $|Q-Q_{n}|=\epsilon$ it is sufficient to take $n\sim\ln\epsilon^{-1}$. Hence, even if the computational time increases exponentially as a function of $n$, $t\sim e^{\#n}$, the increase of $t$ as a function of $\epsilon^{-1}$ is only polynomial, $\ln t\sim\ln\epsilon^{-1}$.

This observation applies to the simulation of interacting fermions by the algorithm introduced in Ref. [12], denoted hereafter by the acronym CDet, for Connected Determinant Diagrammatic MC. This algorithm works directly in the thermodynamic limit since it evaluates the series of connected Feynman diagrams. It exploits two advantages of the fermionic sign: First, for fermions on a lattice at finite temperature, the series has a finite radius of convergence, so that the convergence as a function of diagram order $n$ is exponential; second, a factorial number of connected Feynman diagrams can be evaluated in exponential time using determinants (and a recursive formula).

In general, sign-alternation of observables simulated by MC methods is neither sufficient nor necessary to state that the problem is intractable. One should rather focus on what we will call the “computational complexity problem” (CCP) instead of the FSP. The key question is the one that is most relevant practically222We closely follow ideas expressed by D. Ceperley at the Meeting of the Simons collaboration on the many-electron problem, New York, Nov. 19-20, 2015.: How easily can one indefinitely increase the accuracy of the computed thermodynamic-limit answer? This leads to the following definition of the CCP that can be applied to any numerical scheme. Let $Q$ be the intensive quantity of interest in the thermodynamic limit. A numerical scheme is said to have CCP if the computational time $t$ required to obtain $Q$ with an error $\epsilon$ diverges faster than any polynomial function of $\epsilon^{-1}\to\infty$. The CCP is considered to be solved if

[TABLE]

Note that we consider unbiased methods, i.e., $\epsilon\to 0$ is the difference between computed value and exact value.

In what follows, we show in some detail that CDet solves the CCP, at least at finite temperature and small enough interaction. In parallel, we also discuss the conventional diagrammatic Monte Carlo approach (hereafter denoted by DiagMC) in which the sum over diagram topologies is done stochastically [13]. We then show that the conventional FSP leads to a CCP for path-integral and auxiliary-field QMC.

In quantum Monte Carlo, one typically generates configurations $\mathcal{C}$ according to a conveniently chosen unnormalised probability distribution $P(\mathcal{C})$ that is positive. Any sign alternation is taken into account when collecting statistics, and is absorbed into the quantity $A(\mathcal{C})$ that is being measured. The average with respect to $P$,

[TABLE]

is estimated through

[TABLE]

with $N_{MC}$ the number of MC measurements and $\mathcal{C}_{i}$ the configuration at the $i$-th measurement. By the central limit theorem, the $1\sigma$ statistical error on (3) is given by

[TABLE]

with $\tau_{\rm auto}$ the integrated autocorrelation time and $\sigma_{A}^{\phantom{.}2}=\langle A^{2}\rangle_{P}-\langle A\rangle_{P}^{\phantom{.}2}$ the variance on individual measurements.

We now specialize to CDet and DiagMC. We consider the computation of an observable (e.g. density or double occupancy). For convenience, we make two simplifications regarding the Monte Carlo algorithm. We expect that this does not change the final CCP scaling. The first simplification is that a separate simulation is performed for each order, while the normalisation factors $z_{n}$ (see below) are known. The second simplification is that in DiagMC, rather than sampling the self-energy diagrams and then obtaining observables from the Dyson equation, we consider here sampling the diagrams for the observable (including one-particle reducible diagrams), so that external variables are simply fixed (to zero in space and imaginary-time representation). A DiagMC configuration is then defined by a Feynman diagram topology together with values of the internal variables $X$. In CDet a configuration is defined only by the internal variables $X$ (the space and imaginary-time coordinates of the interaction vertices), while the weight of a configuration is given by the sum over all possible connected diagram topologies connecting the internal and external vertices.

Let us denote the contribution of a diagram of topology $\mathcal{T}$ for fixed internal variables $X$ by $\mathcal{D}(\mathcal{T},X)$. Let $a_{n}$ be the sum of all Feynman diagrams of order $n$:

[TABLE]

with $\mathcal{S}_{n}$ the set of all diagram topologies at order $n$. This can be rewritten in the form of Eq. (2):

[TABLE]

with the unnormalised distribution to be sampled chosen to be

[TABLE]

and

[TABLE]

with the normalization factors

[TABLE]

So, in DiagMC the diagrams are sampled according to the distribution $\left|\mathcal{D}(\mathcal{T},X)\right|$, while in CDet diagrams are grouped together via determinants and in the MC part of the algorithm one samples $X$ according to the distribution $\left|\sum_{\mathcal{T}\in\mathcal{S}_{n}}\mathcal{D}(\mathcal{T},X)\right|$. In what follows we neglect the statistical error on the normalisation factors $z_{n}$ since they are obtained by sampling a sign-positive quantity.

Here we consider fermions on a lattice at finite temperature, so that the radius of convergence of the diagrammatic series is finite [14, 13, 12]. Assuming that we are inside the radius of convergence, the convergence is exponential,

[TABLE]

with $R>1$ a constant. Here and in what follows we omit multiplicative constants and power laws which do not affect the dominant scaling behavior.

The number of diagrams scales factorially with the order $n$. For CDet, however, one takes into account cancellations between different diagram topologies. More specifically, we expect that for fermions on a lattice at finite temperature,

[TABLE]

with $R_{D}$ and $R_{C}$ positive constants.333 Equation (12) is a natural conjecture given Eq. (10); its rigorous proof may be obtained using techniques similar to those of Ref. [15] (J. Magnen, private communication). Equation (11) is plausible since this is the generic large-order behavior for bosonic theories [16].

Let us discuss the behavor of the average sign, $\langle{\rm sign}\rangle:=\langle{\rm sign}\,A_{n}\rangle_{P_{n}}=a_{n}/z_{n}$, as a function of the order $n$. For DiagMC, we see that $\langle{\rm sign}\rangle$ tends to zero factorially, a manifestation of the near-compensation between different diagrams. For CDet, $\langle{\rm sign}\rangle$ tends to zero exponentially in the generic case where $R_{C}<R$.

As a result, the variance on individual measurements behaves as

[TABLE]

which together with Eq. (4) gives for the statistical error bar on the $n$-th order contribution $a_{n}$,

[TABLE]

Here $\tau_{\rm auto}(n)$ is expected to increase at most polynomially with $n$, which we have checked numerically for CDet; we will neglect this $n$-dependence of $\tau_{\rm auto}$ since it will not affect the final scalings. An appropriate dependence $N_{MC}(n)$ of the number of MC steps on order will be specified below.

Note that the relative statistical error $\epsilon_{\rm stat}(n)/a_{n}\propto z_{n}/a_{n}=1/\langle{\rm sign}\rangle_{P_{n}}$ diverges factorially for DiagMC and exponentially for CDet. This can be viewed as a sign problem for diagrammatic Monte Carlo methods, limiting the order that can be reached. On the other hand, the exponential convergence (10), which is only possible thanks to the fermionic sign, implies that reaching very high orders is not necessary, as we now quantify.

The assumption of exponential convergence implies that the systematic error due to the finite diagram-order cut-off $N$,

[TABLE]

decreases exponentially,

[TABLE]

To achieve a final error $\sim\epsilon$, it is then natural to work in a regime where systematic and statistical errors are on the same order. We thus choose $N$ such that $R^{-N}\sim\epsilon$, and we take a computational time $t$ such that the total statistical error is $\epsilon_{\rm stat}\sim\epsilon$. Neglecting correlations between different orders, we have $\epsilon_{\rm stat}^{2}\simeq\sum_{n=0}^{N}\epsilon_{\rm stat}(n)^{2}$, which leads us to choose $N_{MC}(n)$ such that $\epsilon_{\rm stat}(n)$ is $n$-independent. Equation (14) together with Eqs. (11,12) then yield

[TABLE]

where the factor $3$ for CDet comes from the fact that the computational time per MC-step is $\sim 3^{n}$, because of the recursive formula that needs to be evaluated in order to eliminate disconnected diagrams [12]. As a result, for DiagMC most time is spent sampling the highest order, while for CDet this is the case only for $R_{C}<\sqrt{3}$. Finally, we get

[TABLE]

Hence polynomial scaling is nearly reached with DiagMC, and is achieved with CDet. The exponent for CDet is given by

[TABLE]

if $R_{C}<\sqrt{3}$.

The case $R_{C}>\sqrt{3}$ is particularly instructive. Here, most time is spent sampling low diagram orders, and one has $\alpha=2$, which is the best scaling one can achieve in any Monte Carlo computation. We thus conclude that fermionic sign—all by itself—does not necessarily leads to any qualitative effect on the scaling of computational time with $\epsilon$.

The above scalings are not purely academic considerations, as we illustrate with an example for CDet. We analyse the computation, reported in Ref. [12], of the pressure of Fermi-Hubbard model in two dimensions. The diagrammatic scheme is a bare series with bare tadpoles taken into account through a shift of the chemical potential. The Hubbard parameters are: interaction $U=2$, chemical potential $\mu=0.55978$ and inverse temperature $\beta=8$ (with hopping $=1$); this corresponds to a density $n=0.87500(2)$. Figure 1(a) shows that $|a_{n}|$ approaches an exponential behavior $R^{-n}$ with $R=2.5(1)$, while Figure 1(b) shows that $z_{n}$ approaches $R_{C}^{-n}$ with $R_{C}=0.75(3)$.444For the density and the kinetic energy, we find the same value of $R_{C}$ within our error bars. We can make three important observations. First, the exponent $\alpha=3.8(2)$ is not too large. Second, we clearly reach the asymptotic regime where Eqs. (12,16), and therefore also Eq. (19), are valid. Third, $|a_{n}|$ at low $n$ is ${\sim}100$ times larger than the extrapolation to low orders of the large-order behavior shown by the straight line in Fig. 1(a). These three observations explain why it was possible to obtain a ${\sim}10^{-6}$ relative accuracy for the pressure in Ref. [12]. Interestingly, the second and third observations hold independently of $U$ (with $\mu(U)=\mu_{0}+Un_{0}/2$ as in [12]). In contrast, $\alpha$ diverges when $U$ tends to the critical value $U_{c}=2\,R(U=2)\simeq 5.1$ such that $R(U_{c})=1$. Divergent-series summation methods may allow to approach this point and even to go beyond it; we leave this for future study.

To avoid possible confusion, let us remark that we do not claim any connection between our results and the computational complexity theory of computer science. In this theory, a ‘problem instance’ is defined by $\mathcal{N}$ parameters, and the P complexity-class is defined by polynomial scaling of computational time with respect to $\mathcal{N}$ for $\mathcal{N}\to\infty$ (in the worst case with respect to all possible instances). In the spin-glass problem discussed in [17], an instance is a disorder realisation, and $\mathcal{N}$ is the number of random couplings (which happens to coincide with the system volume). In contrast, here we consider problems defined by a fixed (usually small) number of model-parameters (e.g. $T$, $U$ and $\mu$ for the Hubbard model in the thermodynamic limit).

We turn to ‘traditional’ QMC, by which we mean here path-integral or auxiliary-field QMC. We consider the computation of an intensive quantity at finite temperature. Generically the FSP leads to an exponential scaling with spatial volume and inverse temperature of the average sign, and hence, due to Eq. (4), of the statistical error (see, e.g., [1, 2, 4, 10]):

[TABLE]

where $L$ is the linear system size.

Beside the statistical error, we also need to take into account the systematic error $\epsilon_{\rm sys}(L)$ coming from the finite size $L$. The total error $\epsilon$ entering the CCP is $\epsilon\sim\epsilon_{\rm stat}(L)+\epsilon_{\rm sys}(L)$. We assume that finite-size corrections decrease exponentially,

[TABLE]

which is expected generically (away from second-order phase transitions and at finite temperature). For a given computational-time $t$, the optimal strategy is to choose $L$ so that $\epsilon_{\rm sys}\sim\epsilon_{\rm stat}$, which yields

[TABLE]

So for $d>1$ the scaling of $t$ with $\epsilon^{-1}$ is quasi-polynomial and there is a CCP. In one dimension there is no CCP, which is another illustration of the simple observation presented in the introduction.

Apart from these asymptotic scalings, there are also practical advantages of diagrammatic methods over traditional QMC. For traditional QMC, the condition for getting close to the thermodynamic limit is typically $L$ much larger than the correlation length. Equation (21) then yields a computational time $t\propto e^{\#2\beta L^{d}}$ which is often prohibitive, meaning that one cannot get close to the thermodynamic limit, and that one cannot even reach the asymptotic scaling regimes (22,23). The situation is very different in diagrammatic expansions, where as shown by the above example, the asymptotic regime is accessible, and moreover the lowest orders typically set the scale while higher-order contributions are merely corrections.

In conclusion, unbiased numerical methods for solving quantum many-fermion problems should be evaluated on the basis of their scaling of computational time with respect to the final error bar on thermodynamic-limit quantities. The presence of a fluctuating sign does not suffice to say that a problem is intractable by Monte Carlo. Nothing prevents in principle a polynomial scaling of the CPU-time versus the inverse error bar. We demonstrated that such polynomial complexity is indeed achieved by the recently proposed CDet method when inside the radius of convergence of the Feynman diagrammatic series. Since this method offers the possibility to calculate properties of many-fermion systems in polynomial time, it is fair to say that the sign problem has become irrelevant here and a numerical solution to the many-fermion problem is available, at least in some region of parameter space.