Metropolis updates for Diagrammatic Monte-Carlo algorithms from Schwinger-Dyson equations
We describe a general recipe for constructing Metropolis updates for Diagrammatic Monte-Carlo (DiagMC) algorithms, based on the Schwinger-Dyson 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 non-Abelian 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 Monte-Carlo 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 Monte-Carlo algorithms from Schwinger-Dyson 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, D-93053 Regensburg, Germany
Diagrammatic Monte-Carlo algorithms have turned out to be extremely useful for the first-principle studies of quantum systems for which ordinary configuration-space Monte-Carlo simulations based on the sampling of field configurations become impossible due to the sign problem. Typically, Diagrammatic Monte-Carlo algorithms are constructed by using some explicit diagrammatic representation of strong- or weak-coupling 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 on-site or static inter-electron interactions [1, 2] or for Abelian gauge theories . For non-Abelian lattice gauge theories or principal chiral models, however, both strong- and weak-coupling 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  for some recent work in this direction).
On the other hand, it is known that weak- or strong-coupling expansions in quantum field theories can be obtained by iteratively solving the Schwinger-Dyson equations. In contrast to explicit dual representations, Schwinger-Dyson equations are easy to derive in any QFT with continuous field variables. E.g. for non-Abelian gauge theories these are the Migdal-Makeenko loop equations . In these Proceedings we demonstrate how a suitable set of Metropolis updates for a DiagMC algorithm can be obtained directly from Schwinger-Dyson 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 Schwinger-Dyson 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 Schwinger-Dyson 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 slow-down). Application of this algorithm to principal chiral models and non-Abelian lattice gauge theories will be described in detail in forthcoming publications.
2 Stochastic solution of linear equations
Schwinger-Dyson equations can always be represented as infinite-dimensional linear equations on the space of field correlators (which include both connected and disconnected contributions) with all possible values of and :
where is some theory-dependent linear operator and represents the “source” terms in the Schwinger-Dyson equations (typically, the contact delta-function term in the lowest-order Schwinger-Dyson equation). In Section 4 we give a particular example of such representation of Schwinger-Dyson equations for the quartic matrix model.
The solution of (2) can be written as formal geometric series in :
In case one truncates the strong- or weak-coupling 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.
In order to simplify the notation in what follows, let us also define the quantities
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 Metropolis-Hastings 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 .
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 . For the three updates defined above we find the following acceptance probabilities:
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
3 Practical implementation of the stochastic solution of Schwinger-Dyson equations
Typically, for Schwinger-Dyson 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 re-use 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 Schwinger-Dyson 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 Schwinger-Dyson equations as some transformations on the external legs of the diagrams being sampled. The Monte-Carlo 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 Monte-Carlo 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
where we integrate over Hermitian matrices . A full set of linear Schwinger-Dyson equations of the form (2) for this model can be written in terms of the multi-trace expectation values expanded into power series in :
The Schwinger-Dyson equations can be derived in the conventional way by expanding the full derivative in the path integral, and take the following form:
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 right-hand side of equations (4), (4) and (4) are positive, thus no sign reweighting is necessary.
A subtle point in the implementation of Diagrammatic Monte-Carlo for the Schwinger-Dyson 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 Monte-Carlo 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 Schwinger-Dyson 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 right-hand sides of equations (4). These Metropolis updates provide a stochastic implementation of the so-called topological recursion  for the matrix model (4).
In order to test the performance of our Diagrammatic Monte-Carlo algorithm, we consider the expectation values which can be directly related to the expansion of the free energy in the matrix model (4):
On Fig. 1 we compare the results for the coefficients obtained in the above described DiagMC algorithm with the analytic results obtained in . 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 Monte-Carlo 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 Monte-Carlo data is significantly smaller than the exact result. Analysis of the histograms of for these data points suggests that the origin of under-sampling lies in the heavy-tailed 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.
-  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.
-  N. Prokof’ev, B. Svistunov, Worm Algorithms for Classical Statistical Models, Phys.Rev.Lett. 87 (2001), 160601.
-  C. Gattringer, New developments for dual methods in lattice field theory at non-zero density, PoS LATTICE2013 (2013), 002, ArXiv:1401.7788.
-  C. Gattringer, C. Marchis, Abelian color cycles: a new approach to strong coupling expansion and dual representations for non-Abelian lattice gauge theory (2016), ArXiv:1609.00124.
-  Y. Makeenko, A. A. Migdal, Quantum chromodynamics as dynamics of loops, Nucl. Phys. B 188 (1981), 269 – 316.
-  P. V. Buividovich, Schwinger-Dyson equations in large-N quantum field theories and nonlinear random processes, Phys.Rev.D 83 (2011), 045021, ArXiv:1009.4033.
-  P. V. Buividovich, A method for resummation of perturbative series based on the stochastic solution of Schwinger-Dyson equations, Nucl. Phys. B 853 (2011), 688 – 709, ArXiv:1104.3459.
-  W. K. Hastings, Monte-Carlo sampling methods using Markov chains and their applications, Biometrika 57 (1970), 97 – 109.
-  B. Eynard, A short overview of the “topological recursion” (2014), ArXiv:1412.3286.
-  M. Mariño, R. Schiappa, M. Weiss, Nonperturbative effects and the large-order behavior of matrix models and topological strings, Comm.Num.Theor.Phys. 2 (2008), 349 – 419, ArXiv:0711.1954.