Metropolis updates for Diagrammatic MonteCarlo algorithms from SchwingerDyson equations
Abstract
We describe a general recipe for constructing Metropolis updates for Diagrammatic MonteCarlo (DiagMC) algorithms, based on the SchwingerDyson equations in quantum field theory. This approach bypasses explicit duality transformations, enumeration or classification of diagrams and can be used for lattice quantum field theories with unknown or complicated dual representations (such as nonAbelian lattice gauge theories). DiagMC algorithms constructed in this way can still be plagued by the sign problem, which is, however, completely different from the sign problem in conventional MonteCarlo simulations and has its origin in cancellations between diagrams with positive and negative weights. To test the presented approach, we apply DiagMC to calculate the first 7 orders of expansion in the quartic matrix model and find good agreement with analytic results, with the exception of the close vicinity of the critical coupling where the critical slowing down sets in.
Metropolis updates for Diagrammatic MonteCarlo algorithms from SchwingerDyson equations
P. V. Buividovich^{†}^{†}thanks: Speaker. ^{†}^{†}thanks: This work was supported by the S. Kowalevskaja award from the Alexander von Humboldt foundation.
Institute for Theoretical Physics, Regensburg University, D93053 Regensburg, Germany
\abstract@cs
1 Introduction
Diagrammatic MonteCarlo algorithms have turned out to be extremely useful for the firstprinciple studies of quantum systems for which ordinary configurationspace MonteCarlo simulations based on the sampling of field configurations become impossible due to the sign problem. Typically, Diagrammatic MonteCarlo algorithms are constructed by using some explicit diagrammatic representation of strong or weakcoupling expansion series. One constructs certain ergodic set of transformations on the diagram space and accepts each of the transformations in a Metropolis algorithm according to the ratio of the weights of the initial and transformed diagrams [1, 2]. Explicit forms of perturbative expansions are typically easy to construct for simple condensed matter systems with onsite or static interelectron interactions [1, 2] or for Abelian gauge theories [3]. For nonAbelian lattice gauge theories or principal chiral models, however, both strong and weakcoupling expansions become quite complicated at high orders, and it is difficult to classify and parameterize all the admissible diagrams and to calculate their weights (see [4] for some recent work in this direction).
On the other hand, it is known that weak or strongcoupling expansions in quantum field theories can be obtained by iteratively solving the SchwingerDyson equations. In contrast to explicit dual representations, SchwingerDyson equations are easy to derive in any QFT with continuous field variables. E.g. for nonAbelian gauge theories these are the MigdalMakeenko loop equations [5]. In these Proceedings we demonstrate how a suitable set of Metropolis updates for a DiagMC algorithm can be obtained directly from SchwingerDyson equations. After presenting the general prescription in Section 2, in Section 3 we discuss the practical implementations of the algorithm. Finally, in Section 4 we illustrate the presented approach on the example of the SchwingerDyson equations in the quartic Hermitian matrix model, and test it by calculating the first 7 orders of expansion of the free energy. In contrast to the stochastic solution of SchwingerDyson equations presented by the author in [6, 7], the algorithm presented here works even very close to the critical coupling at which the expansions being sampled diverge (of course, exhibiting the expected critical slowdown). Application of this algorithm to principal chiral models and nonAbelian lattice gauge theories will be described in detail in forthcoming publications.
2 Stochastic solution of linear equations
SchwingerDyson equations can always be represented as infinitedimensional linear equations on the space of field correlators (which include both connected and disconnected contributions) with all possible values of and :
(2.0) 
where is some theorydependent linear operator and represents the “source” terms in the SchwingerDyson equations (typically, the contact deltafunction term in the lowestorder SchwingerDyson equation). In Section 4 we give a particular example of such representation of SchwingerDyson equations for the quartic matrix model.
The solution of (2) can be written as formal geometric series in :
(2.0) 
In case one truncates the strong or weakcoupling expansion at some finite order, the above series contain only a finite number of terms, which correspond to a finite number of diagrams contributing up to given expansion order.
We propose to evaluate the series (2) by stochastically sampling the sequences of variables with arbitrary in (2) with probability
(2.0) 
In order to simplify the notation in what follows, let us also define the quantities
(2.0) 
The solution to the system (2) can be obtained as a histogram of the last variable in the sequence , where each occurrence of is weighted with the sign . This sign reweighting puts certain limitations on the use of the method, since we are effectively sampling the series in which the linear operator is replaced by and which typically have a smaller radius of convergence. In the worst case, the expectation value of the reweighting sign can decrease exponentially with , thus limiting the maximal expansion order which can be sampled by the algorithm.
In order to sample the sequences with probability (2), we use the MetropolisHastings algorithm with the following updates:
 Add element:

With probability add a new element to the sequence , where the probability distribution of is . Multiply the sign variable by .
 Remove element:

If the sequence contains more than one element, with probability remove the last element , thus transforming the sequence into . Multiply the sign variable by .
 Restart:

If the sequence contains only one element , replace it by with probability . The probability distribution of is . Set the sign variable to .
Since for the above updates the detailed balance condition for transition probabilities between the sequences and is in general not satisfied, they should be then accepted or rejected with the probability [8]. For the three updates defined above we find the following acceptance probabilities:
(2.0) 
After some algebra one can express the overall acceptance rate as . In order to reach the optimal performance of the algorithm with acceptance close to unity, one can start the simulations with some initial value of , calculate the expectation values for several trial values and then select the value which maximizes the acceptance.
It is also useful to calculate the normalization factor which relates the (sign weighted) histogram of the last element in the sequences and the actual solution of the system (2). From (2) and (2) one can obtain a linear equation for : , where is the probability distribution of the last element in the sequences , with any sign . From the above equation we can finally express as
(2.0) 
3 Practical implementation of the stochastic solution of SchwingerDyson equations
Typically, for SchwingerDyson equations of the general form (2) the coefficients and are different from zero for only a few values of and . Moreover, the space of variables is in most cases infinite. Therefore it is not possible but also not necessary to store and entirely in computer memory. Rather, it is convenient to implement them as a user defined functions which return only nonzero values of and for a given . Given all nonzero values of , the algorithm should calculate and the acceptance probability for the Add index transition. The calculated values can be stored as an ordered sequence (which, in fact, has the stack structure). Since at each step of the random process either a new element is attached to the sequence or the last element is removed, and all other elements are not changed, one can reuse the previously calculated values of for the calculation of the acceptance probability of the Remove index transition.
Let us now turn to the optimal practical way of saving the sequences . In practice, the variables are themselves the sequences of some variables (e.g. lattice coordinates or momenta), and saving them entirely in computer memory is quite impractical. Let us note, however, that for a given sequence the next proposed element only depends on . Moreover, as a consequence of the sparseness of the coefficients and , only a rather small set of new elements might be proposed for a given . Thus in practice it is possible to enumerate all possible updates for a given set of SchwingerDyson equations and characterize them by some small set of parameters (in most cases, only a few integer or real numbers). Instead of keeping the whole sequence in memory, one can store only the topmost element as well as the history of updates which led to a given sequence. If the Remove index is proposed and accepted by the algorithm, the last update should be “undone”.
Sometimes it can be also useful to graphically interpret different terms in the SchwingerDyson equations as some transformations on the external legs of the diagrams being sampled. The MonteCarlo process described above can be then thought of as a graphical editor for diagram drawing which has several commands for updating the diagrams, such as drawing a line or a vertex. At each MonteCarlo step one either chooses a random update command with some probability, or undoes the previous update. It is clear then that the state of the algorithm can be completely described by the current diagram and the history of commands which were used to draw it.
4 Test case: high orders of expansion for the matrix model
In this Section we test the DiagMC algorithm described above on the simplest example of the quartic Hermitian matrix model, which is specified by the partition function
(4.0) 
where we integrate over Hermitian matrices . A full set of linear SchwingerDyson equations of the form (2) for this model can be written in terms of the multitrace expectation values expanded into power series in :
(4.0) 
The SchwingerDyson equations can be derived in the conventional way by expanding the full derivative in the path integral, and take the following form:
(4.0)  
(4.0)  
(4.0) 
Correspondingly, the states of our DiagMC algorithm (variables in the notation of the previous Sections) can be described by sequences of positive integers of any length , labelled by the variable (which has a geometric interpretation in terms of diagram genus). Note that for positive all the terms on the righthand side of equations (4), (4) and (4) are positive, thus no sign reweighting is necessary.
A subtle point in the implementation of Diagrammatic MonteCarlo for the SchwingerDyson equations (4)(4) is that the quantities are not normalizable due to factorial growth of the number of contributing diagrams with the genus and hence cannot be directly interpreted as statistical weights. On the other hand, at fixed the number of contributing diagrams grows only exponentially with the number of vertices and external legs. For fixed this growth is compensated by the powers of the coupling if is smaller than the critical value . To deal with the growth of at large , and we assume that are proportional to the probability to find the state in the MonteCarlo process times the normalization factor which compensates for this growth. After such a redefinition of variables, we can directly identify the Metropolis updates and their weights from the SchwingerDyson equations (4), (4) and (4):
 Create single trace:

, weight .
 Create two traces:

, takes random value in the range , weight . The genus is increased by one. To “undo” this update, one should store in memory the position at which the second trace was inserted.
 Insert line:

, weight .
 Merge two traces:

, weight . To “undo” this update, one should store in memory either or .
 Create vertex:

, weight .
 Split single trace:

, takes random value in the range , takes random value in the range , the genus is increased by one, weight . To “undo” this update, one should store in memory the position of the trace which is being split.
The initial states of the algorithm are either or with probabilities being proportional to and , corresponding to the first terms in the righthand sides of equations (4). These Metropolis updates provide a stochastic implementation of the socalled topological recursion [9] for the matrix model (4).
In order to test the performance of our Diagrammatic MonteCarlo algorithm, we consider the expectation values which can be directly related to the expansion of the free energy in the matrix model (4):
(4.0) 
On Fig. 1 we compare the results for the coefficients obtained in the above described DiagMC algorithm with the analytic results obtained in [10]. Each data point was obtained by averaging over Metropolis updates, which took several hours on a single CPU core. While the absolute value of spans almost 30 orders of magnitude for and (see left plot on Fig. 1), the ratios between the outcome of MonteCarlo sampling and analytic results are equal to one within statistical errors (see right plot on Fig. 1), except for the data points at higher genera and the values of close to , for which the MonteCarlo data is significantly smaller than the exact result. Analysis of the histograms of for these data points suggests that the origin of undersampling lies in the heavytailed distribution of observables, for which the expectation value should be saturated by rare extreme values. The required increase of statistics reflects the critical slowing down of the algorithm near the critical value , where the geometric series (2) are at the edge of convergence.
References
 [1] K. Van Houcke, E. Kozik, N. Prokof’ev, B. Svistunov, Diagrammatic Monte Carlo, in Computer Simulation Studies in Condensed Matter Physics XXI, Eds. D.P. Landau, S.P. Lewis, and H.B. Schuttler (Springer Verlag, Heidelberg, Berlin 2008) (2008), ArXiv:0802.2923.
 [2] N. Prokof’ev, B. Svistunov, Worm Algorithms for Classical Statistical Models, Phys.Rev.Lett. 87 (2001), 160601.
 [3] C. Gattringer, New developments for dual methods in lattice field theory at nonzero density, PoS LATTICE2013 (2013), 002, ArXiv:1401.7788.
 [4] C. Gattringer, C. Marchis, Abelian color cycles: a new approach to strong coupling expansion and dual representations for nonAbelian lattice gauge theory (2016), ArXiv:1609.00124.
 [5] Y. Makeenko, A. A. Migdal, Quantum chromodynamics as dynamics of loops, Nucl. Phys. B 188 (1981), 269 – 316.
 [6] P. V. Buividovich, SchwingerDyson equations in largeN quantum field theories and nonlinear random processes, Phys.Rev.D 83 (2011), 045021, ArXiv:1009.4033.
 [7] P. V. Buividovich, A method for resummation of perturbative series based on the stochastic solution of SchwingerDyson equations, Nucl. Phys. B 853 (2011), 688 – 709, ArXiv:1104.3459.
 [8] W. K. Hastings, MonteCarlo sampling methods using Markov chains and their applications, Biometrika 57 (1970), 97 – 109.
 [9] B. Eynard, A short overview of the “topological recursion” (2014), ArXiv:1412.3286.
 [10] M. Mariño, R. Schiappa, M. Weiss, Nonperturbative effects and the largeorder behavior of matrix models and topological strings, Comm.Num.Theor.Phys. 2 (2008), 349 – 419, ArXiv:0711.1954.