{'fraction_non_alphanumeric': 0.035746325305361946, 'fraction_numerical': 0.009661169001449175, 'mean_word_length': 4.160731244064578, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The LHMC Algorithm for Free Field Theory: Reexamining Overrelaxation

We analyze the autocorrelations for the LHMC algorithm in the context of free field theory. In this case this is just Adler's overrelaxation algorithm. We consider the algorithm with even/odd, lexicographic, and random updates, and show that its efficiency depends crucially on this ordering of sites when optimized for a given class of operators. In particular, we show that, contrary to previous expectations, it is possible to eliminate critical slowing down (z[int]=0) for a class of interesting observables, including the magnetic susceptibility: this can be done with lexicographic updates but is not possible with even/odd (z[int]=1) or random (z[int]=2) updates. We are considering the dynamical critical exponent z[int] for integrated autocorrelations rather than for the exponential autocorrelation time; this is reasonable because it is the integrated autocorrelation which determines the cost of a Monte Carlo computation.


Introduction
Stochastic overrelaxation, and especially its variant usually referred to as "hybrid overrelaxation," is generally considered to be the most efficient algorithm for generating lattice configurations in pure gauge theory. As such it is also frequently used for the purposes of quenched simulations, although at the currently studied (small) correlation lengths the relative improvement over the standard Metropolis or heatbath local algorithms is a priori rather modest.
The idea of generalizing overrelaxation methods to the stochastic case is due to Adler [2], who proposed an algorithm for multiquadratic actions such as free field theory (see also [3]). The initial expectations of improved performance were confirmed by studying the dynamical critical behaviour of the algorithm in this context analytically [4][5][6][7]. The main result of these studies was that the algorithm can be tuned to achieve the dynamical critical exponent z exp = 1 as compared to the generic value z exp = 2 for the standard local Metropolis or heatbath algorithms. Some additional theoretical insights were also obtained in references [8][9].
Several extensions of the algorithm to interacting field theories and subsequent partial modifications were introduced by various authors, including [5,[10][11][12][13][14]. These algorithms were studied in more detail in several numerical works such as [15][16][17][18]. While these almost invariably claim useful improvement, the truly systematic study (such that could make reasonable conclusions about dynamical critical exponents) is still missing.
The purpose of this paper is to analyze the free field behaviour of overrelaxation further. Apart from our original aim of evaluating certain autocorrelations, we obtained some interesting conceptual insights which we feel are worth communicating. It has become clear over the years that the performance of Monte Carlo algorithm should be assessed with respect to a given operator Ω (or set of operators): in particular, the most relevant characteristic then is the dynamical critical exponent z Ω , corresponding to the integrated autocorrelation for that quantity. While the dynamical critical exponent for the exponential autocorrelation time (which we denoted by z exp in previous paragraphs) usually gives us an upper bound, it often happens that z Ω < z exp for most observables.
If we have a situation where for operators of interest, z Ω can be made smaller than z exp , then the latter is somewhat irrelevant for the task at hand. In lattice field theory we actually have a specific group of operators we want to speed up, namely the ones for which the low energy (momentum) states are most important, because these operators are relevant for the continuum limit. We should thus attempt to optimize the overrelaxation algorithm for these quantities, rather than minimize z exp . It has been implicitly assumed in previous studies that zero momentum operators (such as the magnetic susceptibility) actually saturate the z exp -bound and so the above distinction is just a wishful thinking. Our main point in this paper is that while this is true for even/odd updates, where the dynamics of overrelaxation intrinsically couples the low and high frequency modes, it is not true for lexicographic updates, where the modes decouple at large volumes. In fact in the latter case we can tune the overrelaxation algorithm such that critical slowing down is completely eliminated for quantities that depend only on the zero momentum mode. We emphasize that this is a different tuning than the one leading to z exp = 1; in fact, with this tuning z exp = 2. However, in this case the z exp -bound is actually saturated by high frequency quantities most of which we are not interested in.
In the light of the above discussion, the "conventional wisdom" about the inability of local algorithms to perform better than z exp ≈ 1 (see for example [9]) can be rather misleading. While the statement as such might be true (but has not been proved even for free field theory), it does not tell us much about how the local algorithm will perform in the specific case of interest.
The second point we want to make concerns the introduction of the additional noise in the overrelaxation update. In particular we attempted the same optimization for a scheme in which the updated site is chosen at random from a uniform distribution. It turns out that in this case it is impossible to tune the overrelaxation parameter to reduce the critical slowing down for zero momentum quantities; in fact we cannot do better than z Ω = 2 for these operators. This means that overrelaxation can easily lose its magic if we introduce extra noise in the procedure.
In [1] it was shown that the overrelaxation algorithm for free field theory is a special case of the class of Hybrid Molecular Dynamics algorithms where the fields are updated one site at a time using an analytic solution of the equations of motion. This rather suprising and elegant connection is especially intriguing in the light of the fact that, as shown in [14], a corresponding algorithm can be found for gauge theories. Here we adopt this point of view and will frequently refer to free field overrelaxation as Local Hybrid Monte Carlo (LHMC).
We start in section 2 by formulating the LHMC algorithm and developing the techniques to calculate the integrated autocorrelations for typical zero momentum quantities, namely magnetization and susceptibility. In section 3 we analyze the performance of LHMC with even/odd, lexicographic, and random updates. We also give a proof that the integrated autocorrelation for the magnetic susceptibility vanishes at large volume with optimally tuned overrelaxation parameter when using the modified lexicographic update. In Appendix A we give an explicit asymptotic calculation of the update matrix for the lexicographic scheme, and in Appendix B we develop a formalism for bounding the finite volume corrections using methods of functional analysis.

Analytic Solution
We wish to consider a real scalar free field theory described by a functional integral with the usual action and a flat (Lebesgue) measure for the φ field. The lattice Laplacian operator is defined as whereμ is a unit vector in the µ direction. Since the theory is free it can be diagonalized using the Fourier transformed fields The action simplifies to 1 For notational convenience we will frequently omit the tildes and follow the convention that Fourier components are implied by subscripts p, q, . . . as opposed to x, y, . . ..

Local Hybrid Monte Carlo Updates
Consider the dependence of the action S upon one degree of freedom only, and the solution of the corresponding equations of motion is in terms of the overrelaxation parameter ζ ≡ 1 − cos ωt. Assembling the field variables and their conjugate fictitious momenta into the vectors we can write the elementary local update (2) in the form here M x and P x are L d × L d matrices, and the subscript refers to the site being updated. More explicitly, we have ω .
2 It is easy to get a factor of two wrong here! Definition (3) implicitly assumes that we have ordered the field variables in a certain way. Sweeping through the lattice in this order, the complete update is given by Note that the form of the matrices M and P depends upon the predetermined order in which we chose to sweep through the lattice. In fact, as we shall discuss, the efficiency of the algorithm (optimized for a given operator) can depend crucially on this order.
In addition to such fixed-order updates we will also analyze the algorithm in which the updated sites are chosen at random. The formulae (4-6) carry over unchanged to this case except that the updated site becomes an uniformly distributed random variable over d L .

Cost of Measuring the Operators
If we are interested in measuring the expectation value of some operator Ω, andΩ is the average over a sequence of T configurations generated by some Markov process, we expect Here Var Ω = Ω 2 − Ω 2 is the intrinsic variance of Ω, and A Ω is the integrated autocorrelation defined as where C c Ω is the connected autocorrelation function. The cost of measuring Ω to a given accuracy depends upon A Ω , and for local algorithms (such as LHMC) this is the only relevant characteristic which we need to compute in order to ascertain the performance of the algorithm. For a large system near criticality where z Ω is the dynamical critical exponent, corresponding to quantity Ω. 3 In this paper we are interested in studying the performance of LHMC when applied to measuring some interesting operators in the context of free field theory. To this end we shall mainly consider the magnetization M and magentization squared M 2 , where we define the magnetization as 3 Note that our definition of the dynamical critical exponent assumes that the large volume limit is taken first.
The expectation values of powers of M are given by In particular, M = 0, M 2 = 1 m 2 and M 4 = 3 m 4 . We thus expect that the measured values will satisfy

Autocorrelations
The general formalism for calculating the autocorrelations of any operator that is polynomial in the fields has been developed in [19]. Here we shall use the main ideas of that approach and derive explicit formulae for operators we are interested in.
For convenience we will work in momentum space. In particular, we assume that the momentum-space representation of matrices M and P that characterize the updating scheme of our choice is known (see equation (6)) and we express the autocorrelations in terms of the corresponding matrix elements.

Magnetization
If we are interested in the Monte Carlo evolution of quantities linear in field variables, we can average the relation (6) over the independent fictitious momenta, and write φ p π ≡ φ p (0) π → φ p (1) π = M p,q φ q π .
Consequently we have where we have used the relation An explicit representation for the normalized connected autocorrelation function C c M (t) may be obtained by introducing an auxilliary variable x and evaluating the generating function On the right hand side we have a rational expression Q(x)/R(x), where the degrees of polynomials satisfy the conditions 4 We can thus perform the partial fraction expansion of Q(x)/R(x) and write 5 Equating the coefficients of x t we can write the autocorrelation function in terms of the coefficients in the partial fraction expansion, namely Note that this also allows us to express the cumulative autocorrelation function in the form and the integrated autocorrelation as

Quadratic operators
Turning now to quadratic operators we consider a generic quantity of the form with the spectral density K q being some function of the momentum. The change of the quadratic monomial φ * p φ q after a single sweep is given by (M p,r φ r + P p,r π r ) * (M q,s φ s + P q,s π s ), which can be averaged over fictitious momenta to give We can express this as where φ Q p,q ≡ φ * p φ q , M Q pq,rs ≡ (M * ⊗ M) pq,rs = M * p,r M q,s , and P Q p,q ≡ r P * p,r P q,r .
The integrated autocorrelation is then given by Performing the sum over t and using the relation Note that the matrix elements of the identity matrix in the quadratic basis are Á pq,rs ≡ δ p,r δ q,s . In the case of M 2 we have K q = δ q,0 , which implies Note the formal similarity of this expression to that of equation (7) for the magnetization.
In what follows we will also refer to the integrated autocorrelation for the energy (action) E; in this case we have K q = 1 2 f 2 q , giving There is obviously a certain redundancy in the matrix M Q , since the basis elements we have used are not independent. In particular φ * p φ q = (φ * q φ p ) * . In practical calculations it is usually advantageous to reduce the basis to its independent subset, thus reducing the dimensionality of the matrices involved. Also, in case of M 2 we can use the basis {φ p φ q } instead of {φ * p φ q }. Notice also that the formalism for calculating autocorrelation functions which we described in detail for the magentization carries over unchanged; calculating the C c Ω (t) involves the partial fraction expansion of A Ω (x).

Three Updating Schemes
The exact analysis of autocorrelations for the linear stochastic update (6) is rather intractable in the general case, i.e., for an arbitrary ordering of the updated sites. However, in special cases the equations simplify and an exact analysis can be performed. In what follows, we will analyze the even/odd, lexicographic, and random updates.
From our discussion in previous sections it follows that our work will split into two steps: First we need to find the Fourier representation of matrices M and P corresponding to the updating scheme in question, and then we must evaluate the relevant formulae for autocorrelations. Note that the explicit form of matrix P is not needed for quantities we are interested in; this, however, is not true in general.

Even/Odd Updates
Consider splitting the lattice into odd and even sites, as on a checkerboard. If we choose to update all the sites of a given sublattice before the sites of the other sublattice, then the result does not depend on the order of sites within the sublattices. Indeed, updating all the even (odd) sites we may write xµ which is one for even (odd) sites and is zero otherwise. The function χ ± (x) has a simple Fourier representation, In matrix notation we have For an even update followed by an odd update, which is the fundamental ergodic Markov 6 Note that cos Note that with even/odd updates the Fourier modes of the field are only coupled in pairs. In particular, the lowest energy mode φ 0 is only coupled to the highest energy mode φ pc ; for this case we will abbreviate β ± 0 ≡ β ± = 1 2 ζ(−1 ± 2d/ω 2 ).

Autocorrelations for M
Using the formula (7) we can now trivially calculate the integrated autocorrelation for M, namely Since ζ ∈ [0, 2], this attains its minimum (as a function of ζ) for ζ = 2, where A M = − 1 2 , giving M = 0 ± 0. This is of course true even though, or rather because, ζ = 2 does not correspond to an ergodic algorithm. If we tune the overrelaxation parameter to ζ = ζ c ≡ (m 2 + 2d)/(m 2 + d) we have A M ≡ 0, implying z M = 0 while the algorithm is ergodic. It should be emphasized though that this does not necessarily mean that the algorithm generates an independent estimate of M after every sweep: indeed, let us calculate the autocorrelation function C c M (t). According to our general discussion of autocorrelations, this requires the partial fraction decomposition of If we are careful we observe that the equations of motion for the fictitious momenta are so it makes no difference whether we evolve the momenta or not. Of course, this just reflects the fact that the momenta on even and odd sites are independent random variables.
in particular, we have If x 1 and x 2 are the roots of the equation R(x) = 0, then Note that since thus the autocorrelation function has an oscillatory behaviour with regions of correlation and anticorrelation. It is only the integrated autocorrelation that sums to zero for any m, giving z M = 0.

Autocorrelations for M 2
The update matrix M is block-diagonal, so it suffices to consider the 4×4 block M Q p spanned by the basis The matrix elements of M and the above basis are real so this can be further reduced, and we use where we have abbreviated (M p ) ij ≡ M ij ; after a straightforward calculation we obtain .

is then given by
This equation indicates that without tuning of the overrelaxation parameter we have z M 2 = 2 as expected. To find an optimal tuning we need to minimize A M 2 with respect to ζ. Since we are only interested in the asymptotic behaviour as m → 0 (ξ → ∞), it suffices to consider just the leading terms as this limit is approached; in particular, we have which is minimized by With this tuning the integrated autocorrelation becomes giving z M 2 = 1. It is not possible to achieve z = 0 as in the case of magnetization. For completeness we give the expressions for the autocorrelation function. Following our general strategy we find The coefficients a 1 , b 1 , a 2 , b 2 are obtained by inserting these expressions into equation (13), while a 3 , b 3 are given by In Figure 1 we plot the autocorrelation function and the cumulative autocorrelation function for M 2 in two dimensional free field theory, together with numerical data for various sample sizes. The error bars on the numerical data were obtained by binning n measurements into √ n bins. The difficulty in estimating the integrated autocorrelation from small data sets is immediately apparent.

Lexicographic Updates
In this section we analyze updates in which the variable φ x−μ is always updated before the variable φ x . Strictly speaking such an updating scheme does not exist on a finite lattice with periodic boundary conditions; for example, if we start with variable φ 1 on a one-dimensional lattice with L sites, then the variable φ 1−1 ≡ φ L will certainly be updated after it. Allowing for these violations on the boundaries results in the class of updates that we call lexicographic after its most common implementation.
To find the update matrix in this case it is useful to write the linear stochastic update (6) in a different form. The local update (2) involves the updated variable and its nearest neighbours; some of the neighbours, however, might have already been updated in the current sweep. We may separate the dependence on the "old" and "new" variables and write note that the matrixP = Á 1 ω ζ(2 − ζ) is different from the matrix P of equation (6), in fact we have It is straightforward to find an explicit form of the matrices O and N for the lexicographic ordering. We have where we have separated the translationally noninvariant part The matrix R represents the violation of our ordering rule on the lattice boundaries; as such it constitutes only a correction to the translationally invariant part. To see this more explicitly and to take advantage of the translational invariance in the infinite volume limit, it is once again convenient to work in momentum space. We have where the translationally invariant parts are given by and where in the last equation we have defined δ µ p,q ≡ ν =µ δ pν ,qν . In momentum space the matrix elements of the translationally noninvariant contribution are explicitly suppressed by a factor of 1/L: this is expected since the violations have support only on the boundary whose size relative to the bulk volume is given by this factor.
A similar suppression should also be explicit in the form of the update matrix: in Appendices A and B we show that Here ζ/ω 2 < 1/d for m > 0, so that the neglected correction is exponentially small in lattice size. The matrix elements of G (l) have an L-independent bound and are of the form where {µ 1 , . . . , µ l } is a subset of the integers {1, . . . , d}. The leading term of equation (18) is translation invariant and corresponds to the situation with no violations of the ordering rule, namely The 1/L corrections arise from ordering violations on sites which have exactly one component having the boundary value (1 or L). We still have translational invariance in the remaining directions, which is expressed by the presence of δ ν p,q : explicitly we have withK µ defined in equation (43) of Appendix A. Similarly the 1/L l corrections originate with the sites with l components taking boundary values, and the product of δ-functions ensures the residual translation invariance in d − l directions.

A M and A M 2 at Infinite Volume
The analysis of autocorrelations becomes straightforward in an infinite volume because of the diagonal nature of the update matrix. For the magnetization we have Note that this result is identical to that of equation (11) for even/odd updates, and choosing we have again A M ≡ 0, giving z M = 0. Nevertheless there is an interesting difference: the autocorrelation function takes the form of a single exponential, namely giving the exponential autocorrelation time τ M vanishes at ζ = ζ c , consequently each sweep generates a new configuration in which the value of magnetization is completely decorrelated from the previous value.
For M 2 we have similarly and obviously τ M 2 = τ M /2. The integrated autocorrelation has a minimum at ζ c as defined in (23), where A M 2 vanishes, z M 2 = 0, and complete decorrelation is achieved. We have thus obtained the interesting result that with a lexicographic ordering of sites it is possible to tune the overrelaxation parameter ζ so that critical slowing down is completely eliminated. This is true for any quantity Ω(φ 0 ) which depends only on the zero momentum component of the field: the Monte Carlo evolution of φ 0 is not coupled to the other modes in the infinite volume limit, and for one sweep we have where (14)). Since M D 0,0 = 0 at ζ = ζ c we get thus the autocorrelation function as well as the integrated autocorrelation vanish. This is in contrast to the situation for even/odd updates: for lexicographic updates the Fourier modes decouple completely (in the infinite volume limit) and each mode is updated independently. Consequently, the autocorrelations of any quantity which depends on a single Fourier component can be tuned in an optimal fashion. For the even/odd case the modes are coupled in pairs and critical slowing down for, e.g., M 2 can only be reduced to z M 2 = 1 at best.

Other Quantities at Infinite Volume
Let us now try to understand in a little more detail what happens at ζ = ζ c . Consider the "staggered magnetization" squared, φ 2 pc , which is a function of only the highest frequency mode, at ζ = ζ c we have and thus z φ 2 pc = 2. In other words at ζ = ζ c overrelaxation does not help if we want to measure φ 2 pc : this is not surprising as we have not sped up the high frequency modes; we took advantage of the fact that lexicographic updates allowed us to ignore these modes in the situation where we were not interested in them. On the other hand, if we want to measure φ 2 pc we can tune ζ to ζ ′ c = (m 2 +2d)/(m 2 +3d) and eliminate critical slowing down completely for this quantity. The even/odd updating scheme does not have this flexibility.
If we are interested in measuring both φ 2 0 and φ 2 pc near criticality we are better off performing two lexicographic simulations without critical slowing down than one even/odd simulation. This is an extreme case of what we believe is a generic feature: if we want to optimize simulations with overrelaxation updates for some set of operators we can usually find a better solution with lexicographic ordering of sites than with even/odd ordering at large volumes. The reason is that the even/odd ordering introduces more constraints into the optimization than lexicographic update by coupling the modes in pairs.
It is not easy to find an interesting operator that cannot be optimized well with lexicographic overrelaxation. For the sake of simplicity let us stay in the domain of quadratic operators of the type given in equation (8). As we have already observed, if K(q) is peaked at any single momentum we expect z Ω ≈ 0. If the operator Ω couples to both low and high frequency modes, the contribution to autocorrelations of high frequency modes relative to low frequency ones will be suppressed by the f q factors in equation (10). To illustrate this, consider the simple quantity Ω = φ 2 0 + φ 2 pc ; using formula (10) we obtain We obviously get z Ω = 0 at ζ = ζ c because the the 1/m 2 behaviour of 1/[1 − (M D pc,pc ) 2 ] there is cancelled by the m 4 factor from its coefficient. Examples of important quantities of this type are two point functions.
According to the above argument, the most "dangerous" quantities are the ones for which K(q) favours high frequency modes but does not completely decouple from low frequency modes. The obvious such operator is the energy, where the spectral weights exactly cancel (see equation (1)). In one dimension we get z E = 1 without tuning, z E = 2 at ζ = ζ c , and ζ can be tuned to get z E = 1/2 for ζ = 2 − 2 √ 2m; optimization for even/odd updates also leads to z E = 1/2.
Note that according to the above discussion, we should get z int ≈ 0 for most of the quantities of interest for ζ = ζ c , which optimizes functions of the zero momentum Fourier component. Furthermore, if we insist on measuring many operators with different spectra in one simulation, we can interleave lexicographic sweeps with ζ = ζ c with sweeps at ζ corresponding to z exp = 1. This will give z int = 0 for most of the important operators while ensuring that z int ≤ 1 for all operators.

Autocorrelations in a Finite Volume
While the value of z int is defined by the infinite volume behaviour, it is useful to understand the finite size corrections to autocorrelations. Interestingly, for the magnetization there are no corrections to the infinite volume result: indeed, we have From equation (15) (22) is exact for any volume. In case of M 2 the situation is more complicated. (Á − M Q ) −1 is expected to have a 1/L correction to the translationally invariant part. We have not been able to evaluate the correction in the closed form and we have no reason to believe that the zero momentum matrix element vanishes; in fact, our numerical experiments in one dimension confirm that the correction is of order 1/L as is illustrated in Figure 2. The data shows a very small curvature at large values of L, which corresponds to the small negative intercepts for the linear fits; we do not understand the cause of this subtle effect.

Random Lexicographic Updates
On physical grounds it is obvious that the effect of surface terms on behaviour of A Ω for any quantity Ω will go to zero, typically as 1/L, in the limit of large volumes. A periodic lattice does not have any boundaries, it is only the updating scheme that introduces a boundary as the collection of points at which the lexicographic property is violated; for this reason expressions describing autocorrelations do not have exact translational invariance. An obvious way to enhance translational invariance is to randomize the choice of starting point for each sweep through the lattice. For this case we are able to provide a relatively simple direct proof that surface effects only give small corrections to the integrated autocorrelation even for quadratic operators.
The boundary is completely specified by the site that is updated first in the sweep; if x ≡ (x 1 , x 2 , · · · , x d ) is this first site then a site y belongs to the boundary if y j = x j or y j = x j+L−1 for at least one j. We thus need to give every lattice site an equal chance of being the starting point of the lexicographic chain.
Consider the representation of lexicographic updates in terms of the local update matrices (6) where the sequence x 1 , x 2 , . . . , x L d defines the lexicographic ordering. If we relabel the sites by a lattice translation x → x − a and apply the same lexicographic update in the new coordinates then in the old coordinates this corresponds to choosing the initial site to be x 1 + a and the update matrix to be The matrix M x+a can be written in the form Using the composition property S We now propose the updating scheme where before any given sweep, this relabelling of sites is performed at random. In other words, the translation vector a will be a random variable, uniformly distributed over d L . The linear update matrix averaged over these translations is given byM In Fourier space we have as expected, the linear update matrix retains only the diagonal part of M.
The quadratic update matrix in the {φ * p φ q } basis is similarly given bȳ (−a) ).
Performing the sum over translations yields M Q pq,rs = M * p,r M q,s δ p−q,r−s , implying that the quadratic monomial φ * p φ q couples only to monomials φ * r φ s satisfying p−q = r − s. Restricting ourselves to the sector with p − q = 0 and defining Φ p ≡ φ * p φ p we have hereM Q,0 is the corresponding L d × L d block of the quadratic update matrix given by the Hadamard product M * · M. We are interested in showing that Let us start by writing the matrix M explicitly as (18) where we have absorbed the exponentially small corrections into the matrices G (l) ; we then have a decomposition We now argue that the "error" matrix E is small at large volumes in the sense that its matrix norm tends to zero in that limit. To show this we use the maximum sum row norm A ≡ max p q |A p,q |. Consider first E (1) , which is a diagonal matrix since M D is diagonal; the matrix elements of M D and G (l) are bounded by L-independent constants, thus there is a constant C 1 such that E (1) ≤ C 1 /L for all L. Turning to E (2) , since there are an Lindependent number of terms in equation (32), it suffices to consider the norm of a general term in that sum. We have where k is the cardinality of the set {µ 1 , . . . , µ l } ∩ {ν 1 , . . . , ν j }; we have again used the fact that the matrix elements of the matricesĜ µ 1 ,...,µ l have an L-independent bound |Ĝ µ 1 ,...,µ l p,q | < C. The last inequality in (33) follows from the fact that k ≤ min(l, j) and consequently l + j − k ≥ 1. We have thus established that E (2) ≤ C 2 /L for all L. Putting the two bounds together we have As the final ingredient we use the standard formula for the error in the inverse, namely this is valid 8 for any matrix A and "error" E and any norm in which A −1 E < 1. Let A ≡ Á−M D * ·M D , then A −1 has an L-independent upper bound and lim L→∞ E = 0, so we have A −1 E ≤ A −1 · E < 1 for sufficiently large L. Formula (34) is thus applicable, and for sufficiently large L; the required result (30) follows.
In the infinite volume limit A 2 M is the same as for the standard lexicographic update (24). Whilst equation (35) tells us that the leading finite volume correction is at worst O(1/L), it is actually O(1/L 2 ) at ζ = ζ c . Indeed, the matrix elements of E (2) are supressed by 1/L 2 and consequently the O(1/L) part is determined solely by E (1) , which is diagonal. Using the explicit form (21) of matrix G (1) this leading correction can be evaluated, and the correction vanishes at ζ = ζ c . The 1/L 2 behaviour of the correction is in accord with numerical results, as shown in Figure 3.

Random Updates
Finally consider an updating scheme in which the sites to be updated are chosen at random. A sequence of L d such local updates will be called a "sweep".
For a single site update we use the update matrix M z given in equation (5) with the updated site z being a random variable, uniformly distributed over d L . In Fourier space we obtain (M z ) p,q = δ p,q + 1 L d e −2πi(p−q)z/L f (q), with For autocorrelations of quadratic operators we define the quadratic elementary update matrix

Autocorrelations for M
The linear update matrix for a random update sweep, M, averaged over the independent random variables z i is given by 8 The derivation involves the expansion of (A − E) −1 in a Neumann series. Explicitly, we have giving the integrated autocorrelation In the small mass limit we have indicating that z M = 2 regardless of the choice of ζ. Even though the modes have completely decoupled and the autocorrelation function is a single exponential, the random choice of the site to be updated completely destroys the coherence and the algorithm cannot be tuned to reduce the critical slowing down.

Autocorrelations for M 2
The quadratic update matrix M Q averaged over the random variables z i is Straightforward evaluation gives

Inserting this into equation (36) yields
we have not evaluated the O(1/L d ) correction in this case. In the infinite volume limit we obtain [Á − M Q ] −1 pq,rs = δ p,r δ q,s 1 − e f (r)+f (s) , and the integrated autocorrelation of M 2 is given by We therefore find that again that the algorithm with random updates cannot be tuned to reduce critical slowing down below z M 2 = 2.

Conclusions
The principal lesson which we may draw from this analysis is that it is possible to reduce the cost of a Local Hybrid Monte Carlo computation by a judicious choice of the order in which the variables are updated, as well as a careful tuning of the amount of randomness which is introduced into the system. While this differs from some previous claims, it is similar to the situation for deterministic Gauss-Seidel linear equation solvers. The relevant autocorrelation for the determination of the cost is the integrated autocorrelation, and it is in terms of this autocorrelation that we define the critical exponent z int . Furthermore, we define z int by A ∝ ξ z int as L, ξ → ∞ and ξ/L → 0. This is to be contrasted with a finite size scaling analysis in which the limits are taken in a different order. For a zero temperature quantum field theory our choice of limits seems more natural.
Most local algorithms are special cases of LHMC, so our analysis is of fairly general applicability if we assume that the behaviour of autocorrelations in asymptotically free interacting theories is similar to that of free field theory. We intend to study this empirically in a future publication.
Of course local algorithms are directly applicable only to theories with local bosonic actions, but they are used within some dynamical fermion algorithms such as that proposed by Lüscher [20].
One issue of practical importance which has not been addressed at all in this paper is the difficulty of implementing the lexicographic scheme on parallel computers. Nevertheless it may well be advantageous to use a local lexicographic scheme on parallel computers, as has been suggested for the case of preconditioning conjugate gradient linear equation solvers [21].
The full update matrix is the product (in lexicographic order) of the subspace update matrices. In particular, where S x,y ≡ Áδ x,y−1 represents a translation in the new coordinate d + 1.
In Fourier space it is convenient to introduce α ≡ e 2πi/L , so we have Note that we only show the indices corresponding to the new dimension. Writing T ≡ t(Á− N (d) ) −1 one can show by induction that the matrix elements of the lth power for l < L are given by which satisfies the following useful identities x j n P l−j (x 1 , . . . , x n−1 ) P l (x 1 , . . . , x n ) = x n P l−1 (x 1 , . . . , x n ) + P l (x 1 , . . . , x n−1 ) The Lth power has one extra term giving We can obtain a simpler recursion relation for (Á − N) −1 , which is related to M through multiplication by the diagonal matrix; indeed, using equations (29) and (15) we get from (37) Finally we need to transform P L−1 (α p , α q , T ) into an explicit form; straightforward algebraic manipulation leads to Recalling the definition of T we have also It is important to note here that the inverse on the right of this equation exists: we know that (Á− N (d) ) −1 exists for all |t| < 1/d, where t enters as a multiplicative factor in definition of N (d) (see equations (15)- (17)). However, in d + 1 dimensions t satisfies a more stringent constraint, namely t < 1/(d + 1), and consequently implying the existence of (Á − 1 1−tα −p N (d) ) −1 for all p. Inserting relations (39) and (40) into equation (38) we get an explicit recursion relation In principle this allows us to calculate the update matrix in any dimension. Consider building a 1-dimensional update matrix from single site updates (d = 0); Obviously, N (0) ≡ 0 and the recursion relation gives the exact result .
Note that the above result has the expected structure and the leading term exactly corresponds to [Á− N (1) D ] −1 with N (d) D denoting a diagonal part of N (d) , defined in (16). Using N (1) we can calculate N (2) etc. The only problem from the practical point of view is the calculation of T L = t L [(Á − N (d) ) −1 ] L . However, since t < 1 this term is exponentially small in L and can be neglected at large volumes. Similarly, the exponentially small terms will be generated (and can be neglected) from matrix multiplications in the recursion relation when using, e.g., the Poisson resummation formula for the matrix elements of the product. While the above program can be carried out to arbitrary number of dimensions, it is sufficient for our purposes to deduce the structure that follows from equations (41) and (42). Indeed, one can easily verify the following conclusions: 1. In any number of dimensions we have The matrix elements of K (l) are bounded by an L-independent bound and have the form where the modified delta function is defined in equation (19), and {µ 1 , . . . , µ l } represents any subset of integers {1, . . . , d} with l elements. This follows from the (42) and the fact that the above structure is invariant with respect to recursion relation (41). Note that the exponentially small terms could also be included in matrices K (l) , but we keep them separately for the practical reasons mentioned above.

The invariant form of the leading term is given bŷ
which corresponds to [Á − N D ] −1 as expected.
3. The invariant form of the first finite volume correction iŝ The corresponding forms for the update matrix M are given in equations (18)(19)(20)(21).

B Boundedness of Finite Volume Corrections
In this appendix we show that the surface corrections for the lexicographic updating scheme do indeed go to zero as the inverse of the lattice size. For the case of quadratic operators we shall assume that the infinite volume limit of det(Á − M * ⊗ M) = 0: we have not yet found a simple proof of this, although it would follow trivially from a proof of the ergodicity of the LHMC algorithm in an infinite volume. The basic idea of the proof is to expand about the infinite volume limit, and thus our proof is naturally expressed in the language of functional analysis. Our argument is essentially that the update operator M is a compact (completely continuous) operator on the Hilbert space L 2 ([0, 2π]) -in fact in one dimension it is a Hilbert-Schmidt kernel -and that the finite volume results may be expanded about the infinite volume ones using the Poisson resummation formula. We have chosen to prove our results following the original method of Fredholm, as this seems to be the most direct approach [22].
We start by proving a simple bound on determinants: Proof. If A has any row which is zero then the inequality is trivially satisfied. Otherwise construct a new matrix B whose rows are B i = A i / A i , then the bound is equivalent to | det B| ≤ 1 for all complex matrices whose rows are normalized. If the rows of B are not all orthogonal then there are two rows u and v such that Corollary 1 If all the elements of the matrix A are bounded, |A ij | ≤ C, then | det A| ≤ (C √ k) k .
Theorem 1 (Fredholm) Let ∆ : l 2 L ( ) → l 2 L ( ) be a matrix all of whose elements are bounded by some (L-independent) constant |∆ ij | ≤ C; then det Á + 1 L ∆ ≤ C ′ where C ′ is a constant which is independent of L.
Proof. Let us introduce a parameter λ and then expand in powers of λ to obtain using Hadamard's bound where the ratio test may be applied to show that the series converges for all λ. The desired result follows upon setting λ = 1.

Definition 1
The minor ∆ (ij) of a matrix ∆ is defined by Corollary 2 (Fredholm) Let ∆ be a matrix as in Theorem 1, then where C ′ is an (L-independent) constant.
Proof. If i = j then the matrix Á + 1 L ∆ (ij) is of the form to which Theorem 1 is applicable. If i = j then introducing the parameter λ and expanding in powers of λ as before we obtain where the ratio test may be applied to show that the series converges for all λ. The desired result follows upon setting λ = 1.

Definition 2
The adjoint Adj M of a matrix M is defined to be the transpose of the matrix whose elements are the determinants of the corresponding minors, that is (Adj M) ij ≡ det M (ji) .
Proof. From the algebraic identity (Cramer's rule) and the fact that Adj Á + 1 L ∆ ij = C ij /L for i = j and Adj Á + 1 L ∆ jj = C jj , which follow from Corollary 2, we have From this we find that C jj = D + O(1/L) ⇒ C jj /D = 1 + O(1/L), where we have made use of the fact that D is bounded below by the constant α. The desired result then follows from substituting this relation into Cramer's rule above.
Lemma 2 If the matrix elements ∆ ij referred to above are taken to be the values ∆ 2πi L , 2πj L of a Lebesgue integrable function then we can define the infinite volume limit of the Fredholm determinant, and if this limit is non-zero then there must exist a non-zero constant α which bounds it from below.
Notice that according to our definition A (∅) is a diagonal matrix.
Corollary 3 There is a permutation of the indices p for each matrix A (α) such that it is block diagonal with #α × #α blocks. Each block satisfies the assumptions of Theorem 2, so the inverse may be expressed in block diagonal form where each block is of the form specified in Theorem 2. The product of Kronecker tensors may be split into four disjoint classes, α ∪ β, β − α, α − β, and α ∩ β: It is immediately apparent that only for µ in the last set of indices (α ∩ β) does the sum over q µ have more than one non-vanishing term, and that for µ in the first set of indices (α ∪ β) the expression vanishes unless p µ = r µ . We thus have with |W pr | ≤ C ′ . The result now follows from the observation that #α ∪ β = #α + #β − #α ∩ β. where G (α) its has support on dimensions α and all its elements bounded by a constant.
One application of these results uses the matrices N, O, N D , O D , and R introduced in the body of the paper, and the matrix ∆ defined as ∆ ≡ (Á − N D ) −1 R.
Corollary 4 M has the form given in equation (18).
Proof. Observe that since there is a basis in which N is a strictly triangular matrix it follows that det(Á− N) = 1 (∀L), and thus we may conclude that lim L→∞ det Á + 1 L ∆ = 0. Therefore ∆ satisfies the assumptions of Theorem 3, so Á − 1 L ∆ −1 has the desired form. From equation (29) we have so we may apply Lemma 3 to show that M also has the desired form.
If the infinite volume Markov process has a unique fixed point then from equation (9) we obtain upon averaging over this fixed point distribution of φ (Á− M Q ) φ Q φ,π = P Q , and the uniqueness of φ Q φ,π implies that lim L→∞ det (Á − M * ⊗ M) = 0. As a consequence of Corollary 4 we may write where the matrix elements of A (α) are bounded by an L-independent constant, and thus from Theorem 3 and Lemma 3 we may conclude that a result analogous to Corollary 4 applies to quadratic operators.