{'fraction_non_alphanumeric': 0.03756314710148552, 'fraction_numerical': 0.00655229330265593, 'mean_word_length': 4.397948164146868, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sign problem in the Bethe approximation

We propose a message-passing algorithm to compute the Hamiltonian expectation with respect to an appropriate class of trial wave functions for an interacting system of fermions. To this end, we connect the quantum expectations to average quantities in a classical system with both local and global interactions, which are related to the variational parameters and use the Bethe approximation to estimate the average energy within the replica-symmetric approximation. The global interactions, which are needed to obtain a good estimation of the average fermion sign, make the average energy a nonlocal function of the variational parameters. We use some heuristic minimization algorithms to find approximate ground states of the Hubbard model on random regular graphs and observe significant qualitative improvements with respect to the mean-field approximation.


I. INTRODUCTION
Finding the ground-state of a quantum system can be recast as an optimization problem by minimizing the Hamiltonian expectation over the space of trial wave functions. In practice, it is important for the efficiency of the variational method, to have a succinct representation of the trial wave functions that accurately describe the ground-state of the quantum system [1 -3]. Nevertheless, finding the optimal variational parameters could be a hard computational task even in one dimension [4] since the objective function is an average quantity computed over the exponentially large Hilbert space of the physical system, let alone the complex landscape of the energy function induced by different sources of frustrations. And the problem is more serious for fermions due to the global nature of the fermion sign [5]. Still, the main strategy to deal with the above variational problem, is to use Monte Carlo (MC) method both in computing the Hamiltonian expectation and in optimizing over the variational parameters [6].
In this paper we further develop the variational quantum cavity method introduced in Ref. [7] to study the ground-state properties of an interacting fermion system. More precisely, for a given instance of the variational parameters, we connect the quantum expectations to average quantities in a classical system of interacting particles or spins, where the interactions are related to the variational parameters. Then, instead of MC sampling, we use the Bethe approximation [8], or cavity method in the replica-symmetric approximation [9,10] to estimate the classical expectations. Within the Bethe approximation, the probability marginals are obtained by an efficient and local message-passing (MP) algorithm [11,12]; the estimated marginals are good as long as the interaction graph is locally tree-like, the classical system is in a replica-symmetric phase and it is effectively mean-field [13]. Some applications of the cavity method in quantum systems can be found in Refs. [14][15][16][17][18][19]. One may find some connections among these papers, the (statistical) dynamical mean-field theory [20], and density-matrix renormalization group [21].
The trial wave functions can be characterized by the type of interactions included in the associated classical system. We usually start from the mean-field (MF) approximation considering only the one-body interactions, and improve on that by adding higher-order interactions to capture the relevant correlations. For bosons, a good estimation of the quantum expectations can be obtained by considering local interactions involving only a few number of particles [22]. As a result, the Hamiltonian expectation is a local function of the variational parameters and we can again utilize the Bethe approximation to estimate the optimal parameters by a higher level MP algorithm [7,23]. In the case of fermions, however, we have to work with global interactions involving an extensive number of particles to deal with the global nature of the fermion sign. Consequently, the average energy becomes a nonlocal function of the variational parameters and we can not exploit the local MP algorithms to optimize over the parameters anymore.
In this paper we take the Hubbard model and propose a class of trial wave functions with both local and global interactions, where the Hamiltonian expectation can be computed by an MP algorithm. Using some heuristic minimization algorithms we find approximate ground-states of the Hubbard model in random regular graphs of degree C = 3. The results are considerably better than the MF predictions, and are close to the exact solutions in small systems. For comparison we also present some results in one-and two-dimensional lattices.
In the next section we give some definitions and use the mean-field approximation to illustrate the main points that are relevant for the following discussions. In Sec. III we introduce the global ansatz of the wave functions and the machinery we need to deal with the global interactions within the Bethe approximation. The numerical data are presented in Sec. IV and finally we summarize the results in Sec. V.

Consider the Hubbard model with Hamiltonian
with index i = 1, . . . , N that labels the sites in the quantum interaction graph E q . The c † iσ and c iσ are creation and annihilation operators for a fermion of spin σ =↑, ↓ at site i. We will work in the occupation number representation |n ↓ ) n N↓ |0 and will take the following order of the sites and spins: (1 ↑, . . . , N ↑)(1 ↓, . . . , N ↓).
We assume that there is a path in E q connecting 1 → 2 → · · · → N representing the ordering backbone. Given a trial wave function |Ψ(P ) = n ψ(n; eP )|n depending on a set of variational parameters P , we write the Hamiltonian expectation as Ψ(P )|H|Ψ(P ) = n |ψ(n; P )| 2 [E 0 (n) + E 1 (n)] with Depending on the trial wave function we obtain different expressions for e ijσ (n) but we always have e i (n i↑ , n i↓ ) = U i n i↑ n i↓ − ν i (n i↑ + n i↓ ).
The goal is to consider µ(n; P ) = |ψ(n; P )| 2 as a probability measure in a classical system and to compute the above average quantities within the Bethe approximation. The classical measure is, in general, represented by µ(n; P ) ∝ a φ a (n ∂a ) with the set of classical interactions E c ≡ {φ a (n ∂a )|a = 1, . . . , A}. Here, ∂a is the subset of variables that appear in interaction a.
In a MF approximation, we take a factorized trial wave function including the Gutzwiller interactions [24], with complex parameters K i , B i↑ , and B i↓ . This results in the following classical measure µ(n; P ) ∝ i exp 2K R i n i↑ n i↓ + σ 2B R iσ n iσ ≡ i µ i (n i↑ , n i↓ ). By superscript R, we mean the real part of the parameters. Given the above measure, we find e ijσ (n) = −t ij Re δ n jσ n iσ ,01 e −∆nσψ + δ n iσ n jσ ,01 e +∆nσψ , where we defined ∆ nσ ψ = (K * i n iσ − K * j n jσ ) + (B * iσ − B * jσ ) with↑ =↓ and vice versa. Here we can easily compute the average local energies, e.g., There are a few points to mention here: First, the only difference with a bosonic system is the sign term (−1) i<k<j n kσ . It is clear that, in the absence of this sign and for t ij ≥ 0, we can minimize the average energies by setting the imaginary part of the parameters to zero.
When the sign term is present or the t ij take different signs, one can show that, starting from real parameters, one always remains in the real subspace of the parameters following a gradient descent algorithm. This is true not only for the MF wave function, but also for the class of wave functions that we consider in this paper. Second, in the MF approximation, the average of the sign term is exponentially small in the number of sites k between i and j. This suggests that smaller average energies are obtained by maximizing the overlap between the ordering chosen in the trial wave function and the quantum interaction graph E q . Moreover, the density profile would also depend on the ordering unless the parameters are constrained to respect the system's translational symmetries. This artifact of the MF approximation has to be cured by adding interactions to the classical interaction graph E c to correlate distant variables along the ordering backbone. And finally, due to the sign term the average energy is not a local function of the parameters. This sets some restrictions on the optimization algorithms that we can use in minimizing the Hamiltonian expectation.

INTERACTIONS
The simplest interactions to improve the MF approximation are local two-body or Jastrow interactions J ijσ n iσ n jσ [25]. It is not difficult to guess that these interactions are not enough to capture the sign correlations. The interaction set could, of course, be enlarged by adding other many-body interactions also including different types of spins. Instead, here, we take another approach by introducing global variables ξ iσ ≡ (−1) i k=1 n kσ . Then the sign term (−1) i<k<j n kσ can be written as ξ iσ ξ jσ (−1) n jσ , which is a local function of the global variables [26,27]. Accordingly, we can have global one-body interactions Θ iσ ξ iσ and global two-body interactions Γ ijσ ξ iσ ξ jσ in the classical interaction graph. We call this set of trial wave functions the global ansatz. In general one could have interactions of type Notice that the above interactions do not necessarily respect the symmetries of the system.
However, by minimizing the Hamiltonian expectation over the variational parameters, we get closer to the ground state of the system and, therefore, minimizing the effect of these asymmetries.
In the following, we consider the global one-and two-body interactions, i.e.: where, for simplicity, we are going to assume E c = E q . As a result, we obtain with ∆ nσ ψ as given before and The average energy is computed with respect to the following classical measure: The indicator functions I ξ i ,ξ j check the constraints ξ iσ = (−1) n iσ ξ (i−1)σ on the global variables if j = i ± 1, and I ξ i fixes the boundary condition ξ 1σ = (−1) n 1σ when i = 1. To estimate the average energy, we resort to the Bethe approximation, writing the local marginals in terms of the cavity ones satisfying the following set of equations [13]: where ∂i denotes the neighborhood set of site i in E c . These are the belief propagation (BP) equations [11] that are solved by iteration starting from random initial cavity marginals. In the replica-symmetric approximation we assume there is a fixed point to the BP equations describing the single Gibbs state of the system. The average e i (n i↑ , n i↓ ) µ is simply computed after the local marginal µ i (X i ), which is computed like µ i→j (X i ) but taking all the neighbors into account. In the other part of the average energy, we need to compute expectations, such as ξ iσ ξ jσ (−1) n jσ δ n jσ n iσ ,01 exp (−∆ nσ ψ − ∆ ξσ ψ) µ . To get around the problem of computing the average of a global quantity we rewrite it as exp(F −F ) ξ iσ ξ jσ (−1) n jσ δ n jσ n iσ ,01 e −∆nσψ μ , introducing the complex measureμ(X) ∝ exp (−∆ ξσ ψ) µ(X) and the corresponding free energyF . The free energy difference in the Bethe approximation is given by where ∆F i and ∆F ij are the free-energy changes by adding site i and the interaction between sites i and j, respectively, that is, and similarly for the complex measure [13]. In this way, we can compute the Hamiltonian expectation for the above class of trial wave functions with a local message-passing algorithm in a time complexity of order N 2 for sparse classical and quantum interaction graphs. Note that a small error in estimating the classical free energies could result in a large error in estimating the average energy due to the exponential factor exp(F −F ).

IV. NUMERICAL SIMULATIONS
Having the Hamiltonian expectation for an arbitrary instance of the variational parameters, we need an optimization algorithm to find the optimal parameters. This is a computationally hard problem, and we have to resort to some heuristic algorithms to find an approximate ground state for the system. Let us start from a local minimization algorithm where, in each step, we fix all the parameters except in a small region of the system and minimize the associated energy contribution. For instance, in case Γ ijσ = 0, we take the subset {K i , B i↑ , B i↓ , Θ i↑ , Θ i↓ } and minimize the following energy: The energy function is chosen as the sum of the average energies that explicitly depend on the subset of the parameters. The index i is selected randomly, and the corresponding parameters are updated. The process ends when no local update can decrease the average energy.
In another algorithm, we use a population of the parameters {P a |a = 1, . . . , N p } and update the population in order to find smaller average energies. More precisely, in each step, we select two sets of parameters (P a , P b ) and find the set P ab minimizing the average energy along the line λP a + (1 − λ)P b for λ ∈ [0, 1]. Then, we replace the maximal member of the population with P ab and change the position of points a and b to somewhere between P ab and the minimal member of the population P min . That is, P a = (P ab + P min )/2 + R a and P b = (P ab + P min )/2 + R b for some random vectors (R a , R b ).
And finally, in a gradient descent algorithm, the parameters are updated as δP = −η ∂ H µ ∂P for some small and positive η's. This means that we need to have the cavity susceptibilities , which can be obtained by taking the derivative of the BP equations.
For instance, we have if I ξ i = 1, otherwise, χ K l i→j (X i ) = 0. Here we defined and C K l i→j is obtained by normalization X i µ i→j (X i )χ K l i→j (X i ) = 0. These are called the susceptibility propagation equations [28]. Similarly, we can write the equations for the cavity susceptibilities for the complex measure defined in the previous section.
We can use a combination of the above algorithms to approach the optimal parameters, for example, the local minimization algorithm followed by the gradient descent algorithm.
In the following, we always start from zero initial parameters or a population of parameters distributed randomly around zero in the case of population dynamics. The algorithms are repeated a number of times to get the best outcome for different realizations of the update process.
Let us consider the Hubbard model on a homogeneous quantum interaction graph where U i = U, ν i = ν, and t ij = t. In the following, we set t = 1. For simplicity, we only present the results obtained with real parameters; in fact by the above wave functions and algorithms we find the same behaviors when we take the imaginary parameters into account.
Starting from the MF approximation, in Fig. 1, we compare the average charge density in a random regular graph of degree C = 3 with that in one dimension; the MF approximation qualitatively reproduces the expected phase transitions as the chemical potential ν increases for fixed U. For U < U c , we observe, in turn, an empty phase, a partially filled metallic phase, and a completely filled insulator phase. The figure only displays densities smaller  Figure 1 shows how the optimal parameters change by the chemical potential for the paramagnetic solution.
Going beyond the MF approximation, we always obtain smaller energies by adding global interactions other than the local two-body interactions. Moreover, even with the global onebody interactions, we obtain results that are comparable with those obtained by the global two-body interactions where it is more difficult to minimize the average energy, and for small and loopy graphs, the approximation errors cancel out the gain from the global two-body interactions. Here, we present the results for the case Γ ijσ = 0 where the classical interaction graph has no loops, thus, the Bethe estimation of the Hamiltonian expectation is exact and an upper bound for the ground-state energy; this is the case even if we had the global on-site interactions Υ i ξ i↑ ξ i↓ in the classical interaction graph. The global interactions in these wave functions can be considered as effective ones representing higher-order global interactions.
The associated wave functions are expected to describe the physical state of the system for U ≫ t well. In Fig. 2 we compare the minimum energies with the exact solutions on small systems computed by the power method using an infinitesimal imaginary time evolution 1 − Hτ to minimize the Hamiltonian expectation [29]. Indeed, the power method can also be implemented within the variational formalism where, in each step, one has to project the change in the wave function onto the space of the variational parameters; see Ref. [30] for more advanced methods. limit. And in random regular graphs, we observe a shift to smaller chemical potentials for the empty-metal transition point and a shift to larger chemical potentials for the metalinsulator transition point, with respect to the one-dimensional model. One can attribute these shifts to the larger connectivity (here, C = 3) of the random graphs that provide more Finally, we report some preliminary results for the Hubbard model on the 2D square lattice. Figure 6 shows the average energy and density obtained by the MF approximation when we allow for a phase transition from the paramagnetic (P) phase to the ferromagnetic (F) and antiferromagnetic (AF) phases. The former transition happens before the latter one, and the difference between the two energies decreases as one approaches the half-filling density where the AF phase has a smaller energy. Adding the global one-body interactions, we find a transition from the P phase to a mixed phase of F and AF regions, see Fig. 7; the system is more ferromagnetic (antiferromagnetic) for smaller (larger) chemical potentials. Figure 8 displays the charge and spin profiles close to the half-filling density. We observe tendencies towards the holes condensation [31] where a ferromagnetic phase of small density is separated from an antiferromagnetic phase of higher density.
In Fig. 9, we display the time t e we need to compute the Hamiltonian expectation and the local average quantities in the 2D square lattice with the global one-body interactions, given the set of parameters K i , B iσ , and Θ iσ . The computation time for updating the population of the parameters (one sweep) is N p t e , and, in practice, we need a few hundred sweeps of the updates to reach the stationary state. We are reminded that, as long as one only considers the global one-body interactions, the algorithm gives the exact Hamiltonian expectation and, therefore, a good upper bound for the ground-state energy. For comparison, we also show the computation time for a diluted lattice where a small fraction of the global two-body interactions are present in addition to the ones in the ordering backbone.

V. DISCUSSION
The Hamiltonian expectation can be computed by an efficient and distributive messagepassing algorithm, which is asymptotically exact on random and sparse interaction graphs.
To obtain a good estimation of the average fermion sign, we have to work with the global interactions in the classical system. Unfortunately, this makes the average energy a nonlocal function of the variational parameters, resulting in an optimization problem that is not amenable anymore to local message-passing algorithms. Moreover, the performance of any optimization algorithm strongly depends on the quality of the estimated average energy or the approximation method that is used to take the average of the energy function. The study can systematically be improved in both directions by considering more accurate inference algorithms using generalized Bethe approximations [32,33] or incorporating replica-symmetry-breaking and more sophisticated optimization algorithms to find the optimal variational parameters.

ACKNOWLEDGMENTS
We are grateful to A. Montorsi, M. Müller and S. Sorella for helpful discussions. Support from ERC Grant No. OPTINF 267915 is acknowledged.