Efficient inference in the transverse field Ising model

In this paper we introduce an approximate method to solve the quantum cavity equations for transverse field Ising models. The method relies on a projective approximation of the exact cavity distributions of imaginary time trajectories (paths). A key feature, novel in the context of similar algorithms, is the explicit separation of the classical and quantum parts of the distributions. Numerical simulations show accurate results in comparison with the sampled solution of the cavity equations, the exact diagonalization of the Hamiltonian (when possible) and other approximate inference methods in the literature. The computational complexity of this new algorithm scales linearly with the connectivity of the underlying lattice, enabling the study of highly connected networks, as the ones often encountered in quantum machine learning problems.


I. INTRODUCTION
Many relevant systems in statistical mechanics, combinatorial optimization and quantum information theory can be mapped to realizations of the so-called transverse field Ising model (TFIM) [1,2]. Notable examples include Kitaev 1D chains [3,4], many NP-complete problems such as 3SAT [5] and the foundations of the whole field of adiabatic quantum computation [6][7][8]. Given its numerous applications, the study of the properties of this lattice spin model has attracted considerable attention in the last decades. Thanks to that, much have been learned about the nature of quantum phase transitions, especially for the 1D case. However, despite all efforts, doing inference (i.e. computing observables) remains very hard in the general case. The problem stems, as usual, from the exponential increase in size of the Hilbert space with the number of variables. There are two potential ways around the dimensionality curse. One is using a quantum computer to directly implement the required Hamiltonian and let Nature deal with the Hilbert space. Unfortunately, the development of a general purpose quantum computer seems to be at least some decades in the future and near term devices are challenged by decoherence effects and thermal fluctuations. The other option is to make use of suitable approximations and make computations in a classical device.
The quantum cavity method (QC), introduced in [9,10], enables a reduction of the complexity of the inference tasks for a TFIM. It combines an imaginary-time path integral expansion of the density matrix with the cavity formalism, well known in the study of classical disordered systems [11]. As in the conventional cavity approach, QC establishes that all the information of the system is encoded in the set of single variable cavity distributions for interaction networks lacking short loops. This is not true for general topologies but QC remains a useful mean field approximation in that scenario.
The set of coupled equations that determine the cavity distributions can be, at least in principle, solved by an iterative procedure, formally equivalent to the message passing algorithm known as Belief Propagation (BP) [12]. In practice, the implementation cost is similar to a population dynamics algorithm [13], with the important difference that in the present case the sampled variables are trajectories and not real values. There are problems for which the inference calculations must be repeated many times, e.g. for the learning loop of a quantum Boltzmann machine [14,15]. Therefore, it is highly desirable to develop time-and memory-efficient (albeit approximated) methods to solve the QC equations. This is the main purpose of the present work. By focusing of the structure of the QC equations and the properties of the cavity distributions we will be able to derive a fast, precise and flexible solution to the inference problem for the TFIM. To put our analysis in context, we will consider the predictions of other previously developed inference methods. Moreover, we will show how these approximations fit into the QC formalism.
Following the introduction of QC, approximated solution schemes were devised to further reduce its computational burden. Already in [10], the static approximation was proposed as a cheaper way to obtain qualitatively correct results. The static numerical predictions, however, are far from ideal. Of special interest to us is the projective procedure employed in [16], by which the information in the (cavity) distribution is compressed to just a couple effective field parameters. Remarkably, the projected cavity method (PCM), developed therein, is quite precise in the numerical estimation of the observables and the critical temperature of sparse ferromagnetic lattices. Two drawbacks of PCM are the disadvantageous scaling (exponential) of the running time with the system connectivity and that it is specially fine tuned for longitudinal observables, somewhat at the expense of less transversal precision. Other approximations considered in [16] are based on more qualitative mean field arguments. This is the case of the so-called naive mean field and the cavity mean field (CMF) methods. The CMF proposal is simpler and faster than PCM but, as we will see, it does not reduce to the cavity solution in the classical limit of no tranverse field. This should be a minimal requirement to all approximate solutions of the quantum cavity equations that aim at being useful for a wide range of values of the external parameters.
An inference method applicable to the TFIM is the quantum Cluster Variation Method and the related algorithm called quantum Belief Propagation (qBP) [17,18]. At first glance, the standard derivation based on the minimization of a variational Bethe free energy functional is somewhat removed from the QC formalism. However, we can show that it is possible to establish a connection and consider qBP as an approximate solution scheme of the QC equations.
The rest of the paper is organized as follows. In section II, we give a precise definition of the TFIM and review the quantum cavity approach. We recast the exact cavity equations in a form that is suitable for subsequent approximate treatments. In section III we will introduce the main result of this work, an inference algorithm that we call quantum cavity mean field (qCMF), and discuss the connections with other inference approximations. Numerical simulations in different scenarios and topologies are discussed in section IV. Finally, in section V we provide an overview of the results, some concluding remarks and perspectives for future developments.

II. QUANTUM CAVITY SOLUTION OF THE TFIM
In this section we present a summary of the quantum cavity formalism as applied to the transverse field Ising model on a diluted graph. We will introduce the basic Hamiltonian of the system, the corresponding density matrix and sketch the mapping that turns the quantum statistical problem into a classical probability distribution. We skip the detailed derivation and refer the interested reader to the original papers [9,10].
The TFIM is defined by the following Hamiltonian: The transverse fields B i introduce non-commuting terms to an otherwise classical (diagonal) Hamiltonian. In the context of the cavity method, the interaction topology specified by the symmetric matrix J ij corresponds to a diluted, tree-like network. Realizations of such geometries include Cayley trees as well as random regular and Erdos-Renyi graphs. One can also apply the cavity method in other situations, such as in finite dimensional regular lattices. In that case, instead of producing an exact solution, it is interpreted as an structured mean field approximation. The quantum cavity method is based on a Suzuki-Trotter expansion of the density matrixρ = 1 Z exp −βĤ [19]. This procedure reveals a mapping between the quantum system defined by (1) and an equivalent classical model with an extra dimension (the so-called imaginary time coordinate). The classical equivalent has the same interaction topology as the original, with trajectories (paths) replacing the quantum spins in the nodes of the network structure. The statistics of the new system will be given by the joint probability distribution ρ(σ) of all trajectories σ = ( σ 1 , . . . , σ N ). Each path, denoted σ i , is a piecewise constant function taking the values ±1 alternatively in the imaginary time direction, see In equation (3) we have defined m( σ i ) = 1 Ns τ σ τ i and σ i · σ j = 1 Ns τ σ τ i σ τ j , where N s is the number of Suzuki-Trotter slices used in the expansion. The τ index is the imaginary time coordinate, taking values in the range [1, . . . , N s ]. We will, however, normally work in the continuous time limit N s → ∞. Then τ becomes a continuous parameter and the previous sums become definite integrals in the interval [0, 1]. Notice that the energy function E z (σ) does not contain interactions between different Trotter slices. Moreover, the external fields h i act homogeneously in the imaginary time direction. Different slices are coupled only by the w( σ i , B i ) factors. The form of these w factors imply that in the classical regime, B i → 0, the trajectories with non-zero probability are only those which are constant in the imaginary time direction. It is also easy to see in the definition of w( σ i , B i ) that this interaction represents a one-dimensional ferromagnetic chain in the imaginary time direction.
The quantum-classical mapping is completed once we specify the correspondence of quantum operators representing physical observables with functions of the spin trajectories:F → F (σ). The quantum average F = Tr[ρF ] becomes a classical average σ F (σ)ρ(σ). For example, to compute m z i = σ z i one has to average m( σ i ) over the marginal ρ i ( σ i ): See Appendix (A) for the definitions of other observables of interest such as the m x i magnetization and pairwise correlations. The important point to note here is that these quantities typically rely on the computation of marginals of the joint density ρ(σ), that is, on the distributions of small subsets of trajectories. The inference of such marginals is computationally unfeasible, much like in the original quantum system, where the exponentially growing size of the Hilbert space forbids any attempts of computing partial traces of the density matrix. At this point is where the cavity method comes into action.
The cavity approach exploits the structure of the interaction network to simplify the computation of local marginals such as ρ i ( σ i ) and ρ ij ( σ i , σ j ). The functional form of the resulting local distributions include a Boltzmann-like term with the interactions appearing in the original Hamiltonian. The effect of the neighbors is included by extra multiplicative factors. The single spin (trajectory) distribution has the form: here ∂i represents the set of spins that are neighbors of spin i. For interacting pairs we have a similar structure: The µ and M factors above are called cavity distributions. Any given µ i→j ( σ i ) is interpreted as the probability distribution of the variable σ i in a system where the interaction with site j is removed. Contrastingly, although M k→i ( σ i ) is positive for all σ i and can be normalized to a probability distribution, its physical meaning is not so clear in this quantum scenario. A reader familiar with the classical Belief Propagation will identify (6) and (7) as the belief distributions and the M and µ functions will correspond to messages in a factor graph [12].
The cavity distributions µ and M must take values such that the consistency of (6) and (7) holds. The probability function of any pair must be consistent with each single trajectory distribution: The constraints introduced by (8) give the relation between the µ and M distributions: An exact solution of the system of equations given by (9) and (10) is technically unfeasible in the general case, due to the characteristics of the sample space. One possibility is to implement a numerical solution as described in [10]. The basic idea is to represent the µ distributions with sampled estimates, that is, each distribution becomes a weighted collection of trajectories.
In the next section we will describe alternative approaches. We will explore approximations where the complexity of the distributions is reduced with the aid of mean field arguments. Only a reduced number of parameters will be needed to encapsule the statistical properties of the cavity distributions.
As a final note in this summary, it is worth mentioning that it is not strictly necessary to use both µ and M distributions. One can perfectly choose one type at convenience and write everything in terms of that. In fact, previous contributions privilege the use of µ distributions, probably due to the clearer physical interpretation. We think, however, that adopting a combined treatment provides further insights and flexibility regarding where and how approximations can be done. For completeness, we write down the self-consistent equations that are obtained after eliminating the M distributions: and the ones resulting from eliminating µ:

III. APPROXIMATE SOLUTIONS OF THE CAVITY EQUATIONS
Approximate solutions of the cavity equations are usually of a projective type. This means that one looks for a solution within a tractable subspace (e.g. of factorized distributions). The main problem with this approach is that equations (9) and (10) do not preserve the simplified structure in general. For example, if a factorized (in the imaginary time direction) µ distribution is plugged in the RHS of (10), the resulting M will not be factorized in the same way. One needs to project said M back to the easier subspace. Different projections lead to approximations with varying performance. Below we mention some possibilities already explored in the literature, as this will provide the necessary context for our own contribution.
The static approximation [10] essentially amounts to introduce the ansatz µ i→j ( σ i ) = µ i→j (m( σ i )). The probability is projected over surfaces of constant m( σ i ). In other words, all trajectories with the same m( σ i ) get the same probability. Interestingly, the ansatz is correct in the classical limit, where trajectories are constant in the imaginary time direction. However, in the quantum regime this mean field treatment of the intrisically one dimensional structure of µ i→j ( σ i ) [due to the w factor in it, see equation (11)] leads to inaccurate numerical results.
A significantly more precise treatment is the projected cavity mapping (PCM) of [16,20]. The idea in this case is to use a projection of the form One can put this ansatz on the RHS of (11) for the distributions µ k→i ( σ k ) and compute h i→j on the LHS such that both parts of the equation give the same m z i magnetization. The computation in the RHS is equivalent to determining the m z i magnetization of spin i with the cavity Hamiltonian:Ĥ In the LHS, the equivalent Hamiltonian is just −B iσ (13) reveals the fundamental disadvantage of PCM, namely, the scaling of the computation time with the number of neighbors. This is a serious issue if one is interested in using the PCM to compute the approximate statistics of a model with high connectivity.
In the same paper [16], the authors propose what they call a cavity mean field (CMF) approximation to compute the single site statistics. It is a less refined approximation where the cavity Hamiltonian is defined as the single particle operator:Ĥ In this case the projected cavity field h i→j is simply the expression inside brackets in equation (14). The quantity σ z k k\i is computed in the cavity construction corresponding to spin k, that is, using the HamiltonianĤ k\i . This gives a recursive relation for the cavity fields: The computational cost of this approximation scales much better (linearly) with the connectivity of the spins. The simplification comes at a cost: it does not provide the correct classical limit when B → 0. In effect, if we take this limit, we get the iterative equations: whereas the correct classical relations for the cavity fields are [13]: Also, it is not clear from the original derivation what is the relation between the CMF ansatz and the formalism of trajectory distributions. To make this connection clearer we can go back to equation (10) and approximate the average over σ k by making: This can be interpreted as substituting the trajectory σ k in the exponent by its average value m k→i = m( σ k ) µ k→i , computed with the cavity distribution µ k→i . After using this ansatz in equation (9), the effective field h i→j is identical to the CMF result. Although the CMF can be obtained without making any contact with the trajectory formalism, in doing so we can readily find pointers for improvements. This is an attractive goal due to the convenient scaling properties of this algorithm. We note in passing that the naive mean field approximation of [16] can also be obtained from the exact cavity equations if we make the approximation instead of the cavity one µ k→i ( σ k ). The fixed point equations for the magnetizations m i turn out to be the same in both formalisms. Another inference method we want to discuss is a quantum version of the Belief Propagation algorithm (qBP) [17,18]. The usual derivation follows a minimization of the Bethe free energy of the system in a restricted subspace of compatible local density matricesρ i andρ ij . While searching for effective projection strategies, we have found that it is also possible to obtain the qBP algorithm from the full quantum cavity equations. To see this, one can use the M distributions and project to the exponential form M k→i ( σ k ) ∝ exp [βu k→i m( σ k )]. This time, however, the value of u k→i will not be fixed by projecting in (10) or (12). Instead, just as in the variational derivation of qBP, the effective field u k→i is fixed by forcing the consistency relation: One can show that the effective fields u k→i obey the same recurrence equations as in qBP. This quantum extension of BP has an important structural difference with the classical one. For classical BP, in the computation of u k→i only the set of effective fields {u l→k } is needed (with l running over the neighbours of k except for spin i). This means that a change in u k→i would never have an effect on itself if the network is a tree. For classical trees, information flows always in loopless paths. In the quantum case, however, the situation is quite different. Equation (19) uses also information from neighbors of spin i to determine u k→i . It is not hard to see that in this case there will be feedback loops for the flow of information. Numerical simulations discussed in section IV show that this induced loopiness has consequences for the behaviour of the algorithm.

A. The quantum cavity mean field algorithm
In order to account for the inaccuracy of CMF in the small B regime, we should look more closely to the behaviour of the classical and quantum parts of the trajectory distributions. Numerical studies for the static approximation [10] already suggest that the constant trajectories with σ(τ ) = 1(−1) ∀τ bear a significant fraction of the probability weight for small and even moderate transverse field B values. Separating the classical and quantum parts of the µ distributions before making any projective approximation is therefore an idea worth exploring. Using this insight we introduce below a quantum Cavity Mean Field (qCMF) algorithm.
Let I c [ σ, s] be the indicator function that is different from zero only for the classical trajectories σ ≡ s, with s = ±1: I c [ σ, s] = 1 σ(τ ) = s ∀τ 0 otherwise (20) and let I q [ σ] be the complement of I c . We can then separate: where the first and second term correspond to quantum and classical trajectories, respectively. For convenience we will denote c ki (s) = µ k→i ( σ k ≡ s) the probability weight of the constant trajectories. It is clear that q ki = 1− s=±1 c ki (s) equals the trace of the quantum part of µ k→i . With this definition we can renormalize the distribution restricted to the quantum trajectories by writing: and obtain finally: The key idea to outperform CMF is to use (23) back in (10) and project only the quantum part: The cavity magnetization m q k→i is computed using only the quantum part µ q k→i . The functional form of (24) is more involved than just the single exponential term used in CMF but, since M k→i ( σ i ) is a function of m( σ i ) and not of all the details of the trajectory σ i , the problem is still tractable. With the projected M k→i ( σ i ), computed for all k ∈ ∂i \ j, we can find the parameters c ij (s) and m q i→j for µ i→j ( σ i ) and generate a new M i→j ( σ j ) with an expression similar to (24) [see Appendix (B)]. This iteration is continued until a stationary value is reached for all the parameters. Algorithm (1) provides a summary of the process. It is not difficult to prove analitically that the algorithm we are proposing here has the correct classical limit. Also, it is important to stress that the scaling of the computational cost of computing the µ distributions is linear in the number of spins in the cavity construction.

IV. RESULTS AND DISCUSSION
In this section we present a numerical comparison of the quality of the different approximated treatments of the quantum cavity method. We compare the static approximation (SA), the projected cavity method (PCM) and its efficient approximation (CMF), our proposed quantum cavity mean field (qCMF) and the quantum belief propagation (qBP).
Our qCMF can be understood an upgrade of the cavity mean field of [16], which in turn is a simplification of the PCM. Therefore, we want to clarify how our proposal compares with PCM and establish whether it is able to recapture the desired accuracy while preserving a linear scaling of the computational complexity. To this end, we define a simple star-like lattice with a spin i in the center and N − 1 spins connected to it as in figure (2). We then compare the m z magnetizations predicted by the three approximations with the exact value for different connectivities. We use an external field h z = 0.1 to break the z-symmetry of the Hamiltonian and B = 0.7, which is an intermediate value between We study the quality of these approximations for different system sizes N ∈ [3,12] using as a reference the computation of the exact density matrix, see figure (3).
cavity distribution of a border spin, µ j→i ( σ j ), is trivial and does not require any projection. In the absence of i the border spin is isolated. As a consequence, the PCM statistics of the central spin correspond trivially with the exact calculation. It is more informative to look at one of the spins j in the border and compare its observables when µ i→j ( σ i ) (the cavity distribution of the central spin) is computed via PCM or qCMF. Figure (3a) shows that both algorithms accurately estimate the m z magnetization for all system sizes. Calculations also show that qCMF gives better results than PCM for the transverse magnetization m x , especially for high connectivities. This last finding is somewhat unexpected because PCM makes use of an exact diagonalization of the Hilbert subspace of spin i and its neighbors. A plausible explanation is that PCM is fine-tuned to approximate m z through moment matching of the real and projected µ distributions, while m x remains uncontrolled. We also display the results for the CMF algorithm (the one without the classical correction) which, as expected, gives poorer estimates than its counterparts. From this analysis we can conclude that a careful treatment of the balance between classical and quantum trajectories in the distributions µ i→j ( σ i ) is enough to obtain a precise estimate of m z while reducing the computational cost from exponential to linear in the number of neighbors. Moreover, we observe that our approximation (qCMF) is able to outperform the more costly procedure of exact diagonalization for the transverse magnetization.
We now proceed to a more thorough benchmarking on another typical toy model: an infinite tree (or random regular graph) of fixed connectivity and homogeneous ferromagnetic interactions. The connectivity of the graph is c = 3. We show in figure (4) the B vs. T heatmap of deviations of the total magnetization |∆ m| = | m approx − m ref | from the exact values for the SA, PCM, qCMF and qBP approximations. By comparing the total magnetization we can get a better idea of the overall performance and make an assessment that includes the transverse magnetization as well. As a reference value m ref we use the prediction of the cavity method in continuous time computed as in [10]. The first thing we notice is that all approximations suffer in the vicinity of the quantum phase transition, located at B c ≈ 2.232. Close to this region, the best results are obtained by the qBP and PCM. This is in line with the fact that these two approximations are designed to be accurate for m z , the order parameter of this transition. Considering the overall performance, the best results correspond to PCM and the worst to the static approximation. In the regions away from the transition, we observe that the qCMF gives precise estimates, being in this regime faster than PCM and more accurate than qBP. The main source of error for qBP is in the estimate of the transverse magnetization m x (not shown).
The above results suggest that the qBP and qCMF algorithms are reasonable candidates for efficient approximate inference. The PCM is substantially better only near the phase transition. We now test both qBP and qCMF in finite networks to gain more insight on the effect of the system size and short loops. We consider two finite lattices of 12 spins each. One is a tree, where the connectivity of internal nodes was kept at c = 3 for as many nodes as possible, and the other is a random regular graph with c = 3. We show in figure (5) the average longitudinal magnetization m z as a function of the transverse field B. In the plots shown, the temperature of the systems is low (T = 0.1) and the longitudinal external field h z = 0.001 is very small, although not negligible. The exact calculation (dashed line) shows that for both topologies there is a nonzero m z magnetization. This polarization is small due to the finite system size. For many more spins the external h z would be enough to break the symmetry more evidently. Actually, in the thermodynamic limit, the state of the system corresponds to a ferromagnetic phase when B < B c = 2.232.
For the loopy random regular graph (RRG, right) both qCMF and qBP converge to the (incorrect) infinite system ferromagnetic solution. This is a consequence of the structure of the RRG. For local algorithms the topology is indistinguishable from the infinite case and the short loops makes them converge to the selfconsistent solution of the infinite lattice. For the finite tree (left) we have a different outcome. In this case there are no real loops and the qCMF succeedes in converging to the correct results. The qBP fails in this case, converging to a magnetized solution. This behavior is connected to the induced loopiness of qBP, explained in the previous section. The tree and the RRG are both loopy in what concerns to qBP. Another complication due to the loopy behaviour of qBP in a tree is that one cannot exploit this structure to speedup the convergence. With qCMF (as in the classical BP algorithm) one can schedule the computation of the distributions so that only one pass over the lattice is required. For qBP one has to perform a much slower fixed point iteration. These qualitative differences between the tree and the random regular graph are lost for larger values of h z , for which the symmetry of the low temperature states is broken. In that regime both qBP and qCMF approximate well the exact value of the observables.
Summarizing, amongst the methods used to compute the approximate solution of the quantum cavity equations, qBP and the newly proposed qCMF stand out for two fundamental properties. First, both scale linearly with the connectivity of the network, which allows to consider them for more densely connected systems and second, both provide accurate estimates of observables for an ample range of physical conditions and different topologies. For transverse observables, our qCMF is more accurate. The qBP algorithm by construction cannot take full advantage of factorized tree structures.

V. CONCLUSIONS
In this work we have explored different approximate inference methods for the transverse field Ising model and proposed a new ansatz for the solution of the quantun cavity equations. We have shown that the quantum cavity formalism is a starting point for a wide range of approximations, including those in the preceding literature (qBP, PCM, naive mean field, CMF, static). Our proposal for a new algorithm, the quantum Cavity Mean Field, proved to be efficient and accurate in most regimes, with the largest errors in the vicinity of a phase transition. There are several advantages in the use of this method: it is consistent in the classical limit of zero transverse field, the computational cost scales linearly with the connectivity of the lattice, gives better estimates of transversal observables and preserves the local tree-like structure of the network (no induced loopiness). The method is based on a separation of the classical The reference value m ref is computed from a numerical solution of the exact cavity equations [10]. The static approximation gives the worst results whereas the PCM is the most precise overall. Far from the quantum phase transition (T = 0, B c ≈ 2.232), our qCMF is the more efficient alternative to make a precise inference. The PCM and qBP, by design focused on longitudinal m z estimation, are more suited for analysing the critical region where this quantity is the order parameter.
and quantum parts of the cavity distributions. We can speculate that this separation could be useful in the design of sampling strategies for quantum MC simulations given that classical paths are very cheap to sample.
Looking forward to applications, we expect qCMF to be a useful approximation for highly connected lattices, like the ones found in learning problems, e.g. for quantum Boltzmann machines. Also, qCMF can contribute to theoretical approaches to disordered quantum systems. In this context, the advantage with respect to PCM and qBP is that a fixed point intermediate step, inconvenient when averaging over quenched disorder, is completely avoided. An extension of the algorithm presented here to Hamiltonians with higher order interactions in the Z direction should be straighforward. In this case, the distintion between µ and M distributions is even more meaningful and the most convenient graphical representation of the system structure is given by a factor graph [12]. An interesting direction for future work would be to analyze other members of the stoquastic Hamiltonians family, for example, those includinĝ σ x iσ x j interaction terms. This is a harder task since the ST expansion of the corresponding density matrix cannot be decomposed in independent factors as for the TFIM.
FIG. 5: Comparing qBP and qCMF average m z magnetization for a ferromagnetic finite tree (left) and a random regular graph (right) as a function of the transverse field for low temperature. The simulation shows that for a tree, qBP is unable to find the correct solution and gets caught by the effects of induced loopiness. In contrast, qCMF works well in the tree, accounting for the balance between all spins up and down solutions. In the random regular graph, with physical loops, both algorithms fail to provide an accurate result. An increase of the longitudinal field h z would break the symmetry of the low temperature solution and the difference between the exact and approximate results would diminish.
Appendix A: Observables in the path integral formulation ObservablesF that are diagonal in the computational basis of the eigenvectors of |σ z have a very simple form in the path integral formalism [10]. If the operator expression isF = f (σ z 1 , . . . ,σ z N ) then the corresponding observable F (σ) is: f (σ τ 1 , . . . , σ τ N ) (A1) and the average value F = Tr[ρF ] is just F = σ F (σ)ρ(σ). The index τ labels the Suzuki-Trotter slices. These formulas have obvious generalizations in the N s → ∞ limit. Often needed observables are the spin magnetization and the correlations. For the Z magnetization we haveF ≡σ z i . Then, according to (A1): σ τ i and marginalizing the joint probability distribution ρ(σ) we get: In the same way, for the nearest neighbor correlation in the Z direction we obtain: where σ i · σ j = 1 N s Ns τ =1 σ τ i σ τ j , as defined also in the main text. Transversal observables are not so straightforward. The single spin m x i = σ x i magnetization is related to the average of the number of flips n f ( σ i ) in the trajectory σ i : Likewise, the transversal correlation c xx ij corresponds to the average number of simultaneous jumps in trajectories σ i and σ j , denoted n ij f ( σ i , σ j ): with Z i being the normalization constant of the distribution ρ i ( σ i ). To compute m x i we use the relation m x i = 1 β ∂ ln Zi ∂(Bi) to get: