Optimal sampling of dynamical large deviations via matrix product states

The large deviation (LD) statistics of dynamical observables is encoded in the spectral properties of deformed Markov generators. Recent works have shown that tensor network methods are well suited to compute the relevant leading eigenvalues and eigenvectors accurately. However, the efficient generation of the corresponding rare trajectories is a harder task. Here we show how to exploit the MPS approximation of the dominant eigenvector to implement an efficient sampling scheme which closely resembles the optimal (so-called"Doob") dynamics that realises the rare events. We demonstrate our approach on three well-studied lattice models, the Fredrickson-Andersen and East kinetically constrained models (KCMs), and the symmetric simple exclusion process (SSEP). We discuss how to generalise our approach to higher dimensions.


I. INTRODUCTION
The complex behaviour of the non-equilibrium dynamics of stochastic systems can be characterised by studying trajectory ensembles, that is, the set of all possible trajectories alongside the probability that they occur under the evolution defined via a stochastic master operator. This is analogous to standard thermodynamics, where static properties are entirely determined by the equilibrium ensemble of all microstates and the probabilities [1]. Often dynamical behaviour of interest is dominated not by trajectories that are typical under the dynamics, but by "rare events" which are exponentially (in time and in system size) scarce. Studying these rare events is made possible by using the framework of large deviations (LDs) [2][3][4][5][6][7], where in large time limits time-extensive dynamical observables obey a LD principle, and their statistics is encoded in functions which play for dynamics the role that thermodynamic potentials play for statics (see below for definitions).
LD functions can be obtained in principle from a deformation or tilting of the dynamical generator (in the case of continuous-time dynamics) or the Markov matrix (in the case of discrete-time dynamics), through its largest eigenvalue. Obtaining this eigenvalue is not always an easy -or even possible -task, and often one needs to resort to numerical methods. Methods to overcome this difficulty often include techniques based on population dynamics, namely cloning or splitting [8][9][10][11], and importance sampling [12][13][14][15][16] which provide information about the configurations frequently visited by the rare events. Notice that even if one manages to diagonalise the tilted generator (or the Markov matrix), the generation of rare trajectories is non-trivial: while rare trajectories are "generated" by the tilted operator, this is not a proper stochastic operator and these trajectories cannot be directly sampled.
The efficient sampling of rare events can be achieved by searching for another stochastic dynamics which gen-erates trajectories with desirable probabilities that are the same as (or a close approximation to) those of the tilted generator (with any small discrepancy corrected via importance sampling techniques). Methods for doing so currently include optimal control [17,18] and machine learning approaches, where one attempts to "learn" this convenient sampling dynamics [19][20][21]. The optimal choice for a reference dynamics is the so-called generalised Doob dynamics [22][23][24][25][26][27], which generates trajectories with the exact tilting corresponding to the deformed generator. The Doob dynamics thus produces rare trajectories of the original dynamics "on demand". To construct such optimal dynamics, however, requires knowledge of the leading eigenvector of the tilted generator.
Variational tensor network (TN) techniques [28][29][30][31][32][33][34], originally devised as a tool to study quantum many-body systems, are also convenient for studying classical statistical systems [35][36][37][38]. More recently, they have been shown to be useful in the context of LDs in stochastic dynamics [39][40][41][42][43]. In particular, it is often both possible and easy to approximate the leading eigenstate of the tilted generator of a one-dimensional stochastic lattice system using a matrix product state (MPS) ansatz, even those with dynamical (i.e. LD) phase transitions. Recent works have made use of this eigenstate to determine the statistical properties of the dynamics [39][40][41][42][43]. To our knowledge, however, such TN approach has not been exploited yet to sample efficiently rare trajectories. This is what we do in this paper. We present a scheme to use the MPS approximation to the leading eigenvalue of the tilted generator to construct a new dynamics which very closely resembles the optimal Doob dynamics, and we show how we can use this new dynamics to efficiently sample rare events.
The paper is organised as follows. In Sec. II, we review continuous time Markov dynamics and LDs. We also recap how one can apply an MPS ansatz to study KCMs. In Sec. III, we define the Doob dynamics and introduce a scheme to approximate it with a reference dynamics, constructed using an MPS approximation to the leading eigenstate of the tilted generator. In Sec. IV we present the numerical results from our method applied to the three models. We show how our approach can effectively be used to accurately measure the statistics of time-extensive observables. We provide an outlook on possible generalisations and our conclusions in Sec. V.

II. LARGE DEVIATIONS AND MATRIX PRODUCT STATES
In this section we introduce continuous-time Markov dynamics, giving specific examples in the context of kinetically constrained models (KCMs) and exclusion processes. We then also review the framework of large deviations (LDs) and how variational matrix product states (MPS) can be used to determine the LD statistics.

A. Continuous time Markov dynamics for KCMs and exclusion processes
We consider stochastic Markov dynamics which evolves continuously in time. Suppose we have some system with the set of configurations {x 1 , x 2 , . . . , x M } where M is the size of the configuration space. The probability that the system is in some configuration x at the time t is encoded in the probability vector |P (t) = x P (x, t) |x which evolves under the stochastic master equation Here the generator of the dynamics W is given by where w x→x are the transition rates from configuration x to x and R x = x =x w x→x is the escape rate from x. The largest eigenvalue of the generator is zero, with left eigenvector the flat state −| = x x|, and right eigenvector the steady state |ss = x P (x) |x , which describes the probability of finding any configuration at equilibrium. If our system obeys detailed balance, then we are guaranteed that any initial state will eventually relax to some equilibrium state given enough time. Here we assume this to be the case.
We will focus on two broad areas of 1D constrained systems. The first is KCMs (for reviews see [6,44,45]), for which configuration changes are governed by a kinetic constraint which is explicitly encoded in the generator. For concreteness, we focus on the 1D spin facilitation Fredrickson-Andersen (FA) [55] and East [46] models. Both models are defined on a 1D lattice of N binary variables (spins) n j = 0, 1 for j = 1, . . . , N , and configuration changes are only allowed via single-spin flips. The Markovian generators for both models are given by where σ ± i are the Pauli raising/lowering operators acting on site i and c ∈ (0, 0.5] controls the rates at which spins flip, given they satisfy the kinetic constraints where the first only allows a transition if the spin attempting to flip has a neighbouring excitation, and the second only if the neighbouring spin to the left is excited. (For the FA model the constraint is sometimes defined as the projector n i−1 + n i+1 − n i−1 n i+1 , but in practice it makes little difference with the definition above.) The second area we consider are exclusion processes [49,50] -particles hopping around sites on a lattice, with a hardcore exclusion such that we can have at most one particle per site. We focus on the 1D symmetric simple exclusion process (SSEP), adopting the lattice notation we used for KCMs, where now n j = 1(0) implies the site is occupied (empty). In the SSEP, a particle can hop left or right to its neighbouring sites, both with the same rate (γ = 1/2) if the neighbouring site is not already occupied. The generator for the dynamics is where σ a i are the Pauli operators acting on site i. For the entirety of this paper, we will assume open boundary conditions (OBC), which will later reduce the computational cost of tensor network contractions. This formally means that we set n 0 = n N +1 = 0. Furthermore, we impose certain restrictions on the state space. For the FA model, we simply exclude the disconnected zero state n i = 0, ∀i. On the other-hand, we set n 1 = 1 for the East model which ensures the state space remains fully connected on each dynamical site i > 1. Finally, we restrict SSEP such that the total number of particles N p = i n i is fixed, with particle density n p = N p /N which will be assumed to be n p = 1/2.

B. Trajectories and large deviations
Consider some general trajectory ω t = {x 0 → x t1 → · · · → x t K } where the system moves into the configura-tion x ti at time t i and has the total time t > t K . The dynamical activityK [3][4][5][6]56] is a trajectory observable which measures the number of configuration changes for a given trajectory. The probability of observing some activity K can then be calculated as the sum over all trajectories with K configuration changes, and the probability they occur, where π(ω t ) is the probability of observing ω t . For large times, this obeys the large deviation (LD) principle [2][3][4][5] where ϕ(K/t) is called the LD rate function and plays the role of entropy density for trajectories. Alternatively, one can consider the moment generating function (MGF) [2] Z t (s) = which contains equivalent information to Eq. (7) and can be considered the partition function. From Eq. (8), we see that the weighting of each trajectory is the probability that the trajectory occurs, exponentially re-weighted by its dynamical activity. The MGF also obeys a LD principle, where θ(s) is the scaled cumulant generating function (SCGF), whose derivatives evaluated at s = 0 give the cumulants of K scaled by time. The SCGF plays the role of the thermodynamical free energy of trajectories and is related to the LD rate function by a Legendre transform θ(s) = − min k (sk + ϕ(k)) [2]. The MGF Eq. (8) can be expressed as where |in is some initial probability vector and W s is a new operator which we name the tilted generator, and is a deformed version of Eq.
(2) where we tilt with respect to the dynamical observable of interest [2][3][4][5]. For the case of the dynamical activity [3][4][5], we simply tilt the off-diagonals of W with the same factor to obtain The largest eigenvalue of W s is the SCGF θ(s), with associated left and right eigenvectors l s | and |r s . Since l s | in general is not the flat state, W s is not a proper stochastic generator for s = 0 [3][4][5]. If one could exactly diagonalise Eq. (11) to find its leading eigenvalue and eigenvectors, then they would entirely unravel the LD statistics. We now briefly recap how this can be achieved using numerical TN techniques [40][41][42][43].

C. Variational matrix product states
A matrix product state (MPS) is an ansatz for describing vector states of many-body systems [28-30, 57, 58], where each subsystem k has its own rank-3 tensor A k with the dimensions d×D×D. The allowed entanglement within the state is controlled by the bond dimension D [32]. It is often convenient to represent tensor networks in a diagrammatic form using shapes to represent tensorial objects, and (connecting) lines to represent contractions over tensors. For example, the corresponding diagram for an MPS is where each circle corresponds to one of the tensors A k . Similarly, one can also attempt to write some operator O as a matrix product operator (MPO) [33,34,[59][60][61][62].
Operators which act locally on the sub-systems, such as Eqs. (3)(4)(5), can be efficiently described as a MPO. That is to say we can represent them exactly in MPO form with only a small constant bond dimension. The diagrammatic representation for MPOs iŝ MPS allow for the easy and efficient implementation of the widely used density matrix renormalization group (DMRG) method [63,64], an algorithm designed to iteratively minimize the energy of a state E Ψ with respect to some HamiltonianĤ. In the language of MPS [30], we start with some guess at some fixed bond dimension, and sweep through each tensor applying local optimizations with all other tensors fixed. This is done until we reach convergence, which is usually when the change in energy of the state per sweep is small. At the end of the routine, one can efficiently calculate the variance of the state with respect to the Hamiltonian where · Ψ = Ψ| · |Ψ denotes an expectation value. We check to see if it has fallen below some desired value, ; if not, we run the algorithm with an increased bond dimension, where we typically use the state from the previous run as an initial guess. For more details on the workings of variational MPS (vMPS) algorithms, see the reviews [30,65]. Many recent works have shown that vMPS algorithms are very effective for studying the LD statistics of classically constrained systems which obey detailed balance [40,41,43]. In particular, if we write the tilted generator in a way such that it is Hermitian then the state we are searching for is the ground state. This guarantees each update is an improvement upon the last. For dynamics obeying detailed balance, the activity-tilted generator can be brought to a Hermitian form using a similarity transformation that is independent of s [5], For the case of the East/FA models [5], the diagonal operator Q is given by and for the SSEP by Q SSEP = I. The new Hamiltonian H s has the ground state |ψ s with energy −θ(s). The ground state is related to the left and right eigenvectors of W s in the following way [40], where l s (x) = l s |x and r s (x) = x|r s .

III. DOOB TRANSFORMATION AND OPTIMAL SAMPLING
We now define the so-called generalised Doob transformation [19, 22-25, 27, 66], and show how one can use our MPS solution to Eq. (16) to construct a reference dynamics which closely resembles the true Doob dynamics. We then present a method to optimally sample the rare events of our toy models using these new dynamics.

A. Generalised Doob dynamics
The goal is to find a proper stochastic generator which generates trajectories with the same probabilities as those in the tilted dynamics W s , cf. Eq. (11). This can be achieved using the (long-time) generalised Doob transformation [19, 22-25, 27, 66], defined as where L = diag( l s |) is the left eigenvector l s | as a diagonal matrix. It is easy to check that Eq. (21) is annihilated by the flat state −|, which means that W Doob s is a stochastic operator. Its stationary state is The generator W Doob s can also be expressed as a sum of its diagonal and off-diagonal elements Thus our new dynamics has the transition rates and escape ratesw respectively. That is to say the transition rates are reweighted by e −s and by some ratio l s (x )/l s (x) which depends on the structure of the configurations, and the escape rate is shifted by θ(s).
We now consider some general time-dependent observ-ableÂ, and ask what is the expectation value in the tilted dynamics, (26) One can now apply importance sampling to arrive at whereπ(ω t ) is the probability of observing ω t in the dynamics generated by W Doob s and · Doob denotes an expectation value with respect to trajectories with probabilities from the Doob dynamics. At a first glance, it might look that we have not gained much from expressing the expectation of A using the Doob generator W Doob s . However, if one calculates the ratio of probabilities in Eq. (27) then the power of this expression becomes apparent.
Let us first consider the original dynamics described by Eq. (2). If we have some system in configuration x, then the probability it flips to some other state x at the time ∆t is It then follows that the trajectory ω t occurs with probability where we have also accounted for the fact that the system must remain in the same state after the final flip for the remainder of the time and the probability of the initial configuration P (x 0 ) (where we assume it is in the steady state). The probability of the trajectory under the Doob dynamics has a similar form, with the substitutions where all but the endpoint factors of l s (x) cancel out telescopically. The ratio of probabilities then goes as where we have usedP ( And so it follows that one can exactly sample the expectation value of a trajectory observable in the tilted ensemble defined by the non-stochastic tilted generator, by sampling it directly from trajectories generated by the stochastic Doob dynamics Eq. (21), up to factors at the endpoints of each trajectory (which become negligible in the long time limit ifÂ is time-extensive). We note that Eq. (32) can also be derived by means of linear algebra, using Eq. (21) and the ratio P (x 0 )/P (x 0 ) (see Ref. [67] for details).

B. Reference Dynamics
While the above shows how to optimally sample if one has access to the Doob generator, which is obtained from the exact minimisation of the tilted generator, we now consider how to approximate it efficiently.
Suppose we have an MPS approximation |ψ ref We then construct the generator of the so-called reference dynamics, which goes as Eq. (2) with the transition rates and escape rates given by respectively. Note that here we have not used Eq. (25) for the escape rates, as these reference dynamics only act as an approximation to the Doob dynamics, and thus would not give a true stochastic dynamics. In appendix A, we show the steady-state solution to the reference dynamics is given by where dt∆R is the time integral of the difference of escape rates between the reference dynamics and the orig- We can estimate a sampling error when using Eq. (37) in the following way [66]. First, let us assume the effects of the time-edge factors is negligible (as they are not exponential in time) and try to sample the quantity where R(ω α ) = e sK(ω α )+ dt∆R(ω α ) is the umbrella which compensates for change in sampling dynamics and we estimate for a fixed number of samples, N sp . The variance of Eq. (38) gives a way to quantify the sampling error, In appendix B we show The last approximation holds for δE small enough (tδE 2 1). In Eq. (40), δE 2 is the calculated variance on our MPS approximation of the leading eigenvector, see Eq. (15).

C. Simulating trajectories
We are now in a position to efficiently simulate trajectories from our reference dynamics. The sampling of trajectories from a classical generator is usually achieved using a continuous time Monte Carlo (CTMC, otherwise known as the BKL algorithm) [68]. Given that our system is in some configuration x at time t , we need to calculate the next jump in the trajectory. That is, we need to decide the next configuration the system will move into, and the time it does so. Calculating this can be split into five separate steps: 1. Find each configuration x the system can move into from x.
2. Calculate the transition rates w x→x for each x .
3. Calculate the escape rate R x as the sum of all transition rates.
4. Randomly choose one x , each with the probability w x→x /R x 5. Randomly choose the jump time ∆t with probability P (∆t) = R x e −Rx∆t .
By starting at a configuration sampled from equilibrium (which in the case of the reference dynamics can be efficiently done using the MPS |ψ ref s [69,70]), or otherwise, one can simply repeat this procedure until some total time t has elapsed.
We can use this method for our reference dynamics, where the only step that needs slight adjustment is the second. While one must still calculate the transition rates of the original dynamics in the usual way, we must also calculate the left vector components l ref s (x) and l ref s (x ). Let us assume the former is carried over from the previous jump in the algorithm. Then all one needs to do is calculate each l ref s (x ). We start by noting that any configuration x can be written in MPS form with bond dimension 1, and then we can simply calculate the left component as a MPS-MPS contraction The transition rates for the reference dynamics are then calculated using Eq. (34), and the method proceeds as before. The total computational cost for calculating each l ref s (x) is O(D 2 N ), and thus the total cost of each Monte Carlo (MC) step is O(D 2 N N F ), where N F is the total number of configurations x for a given step.
Let us now consider our KCMs where we have singlespin flip dynamics. We first note that the number of possible configuration changes from x is bounded by the number of sites, that is, 1 ≤ N F ≤ N . Using the method described above, the computational cost for each step is at worst, quadratic in the system size. However, by realising that the tensor network contractions l ref s |x and l ref s |x are identical apart from just one tensor (corresponding to the spin which would flip), we can reduce the computational cost by recycling partial contractions from the edges. We first need to identify the first and last sites on the lattice which are able to flip, which we label i l and i r respectively. In a similar fashion to variational algorithms, we then contract from the left edge of tensor network l ref s |x up to i r −1, and saving each tensor block along the way.

→ →
We do the same but from the right and up to i l + 1. This initialization of partial contractions has a one-time cost of Calculating each l ref s (x ) at site j is then easy. We just contract our remaining tensors at site j with the previously saved left and right blocks, This is done for each possible site which can flip, and thus entails a computational cost O(D 2 N F ). Once a choice is made for which site to flip, which we will label i, we must update the blocks of partial contractions up to (the now possibly different) i l and i r . Note that this time we do not have to start from the edges of the MPS, but just from site i as the previous partial contractions that come before do not change. The total cost of updating the partial contractions is The total computation cost for each MC step is the sum of the cost for calculating each l ref s (x) and updating the partial blocks after a choice is made, Consequentially, the cost of each MC step is reduced to one which is at most linear in system size.

A. Approximating the Doob Dynamics
We put to the test the general method presented above by approximating the Doob dynamics of each model defined in Sec. II. We show that the Doob dynamics is well estimated using the MPS reference dynamics, and can even be well approximated with truncated MPS.
Each of the three models is known to exhibit a trajectory phase transition (when tilted against the activity) for long times and in the thermodynamic limit N → ∞, manifested in the SCGF θ(s) at s = 0 with a discontinuous drop in the dynamical activityK(s) = −θ (s)/N [3,5,43,52,54]. We call the dynamical phase for s < 0 the active phase, and that for s > 0, the inactive phase. One is able to do a detailed investigation of this firstorder phase transition by considering the finite-size scaling of the model [40,43,53,71,72]. We can estimate a critical point s c (N ) 0 by finding the peak of the dynamical susceptibility χ(s) = θ (s), which shows a drastic change in a small region around the transition point.
We start by taking the usual approach of approximat-ing the ground states |ψ s using vMPS. That is, we run the algorithm allowing the bond dimension to increase until the variance of the energy (with respect to the Hamiltonian) falls sufficiently, cf. Eq. (15). The resultant MPS is then used to construct the reference dynamics, which approximates the Doob dynamics to a high accuracy, as explained in the previous section. Note that because the vMPS tries to keep entanglement as low as possible, for s > s c (N ) the approximated ground state exhibits localisation at just one edge of the system [40]. While for the East case this corresponds to the structure of the ground state in the sector with fixed occupation 1 in the leftmost site, the FA and SSEP models have reflection symmetry, spontaneously broken for s > 0 and large N . Thus, in order to maintain the symmetry in the latter two cases, we construct an MPS which is a superposition of the result from vMPS and its spatially reflected state to obtain our dynamics in the inactive phase. 1. Direct sampling with the reference dynamics but without re-weighting We first check that the CTMC algorithm with our MPS reference dynamics gives the expected results. We do this without using the trajectory re-weighting, cf. Eq. (37). This amounts to only considering infinite-time dynamics, and assuming that our approximation is actually exact. Despite this strong assumption, we find that it produces excellent results as shown in Fig. 1. The expected dynamical activity (per unit site and time, dashed lines) can be calculated as a TN contraction over our MPS and MPO,k The same quantity can be calculated on a trajectory level (symbols) by counting the total number of configuration changes, K and taking its time (and spatial) average, where t is the run time for each trajectory. We show results for each model, for a range of system sizes of N ∈ [20, 400]. The expected and measured results have excellent agreement. This simplified algorithm struggles most around the transition point, s c (N ), due to the required large bond dimension (see Refs. [40,43]). We also show representative trajectories for the active (s = −1) and inactive phases (s = 1). Each model excellently demonstrates the difference in dynamics between the two phases. The active phase displays very rapid changes with structures that allow for unconstrained dynamics. For the FA and East models this means having a large number of excitations, while SSEP requires particles to be spaced apart. Conversely, this inactive phase has just few configuration changes with highly constrained dynamics. This means minimizing the number of excitations for the FA and East resulting in the dynamics responsible for the so-called "space-time bubble" in local regions of space [5,47,73], while for SSEP we restrict the activity by clustering the particles [54,74]. To our knowledge, direct dynamical sampling of trajectories for these systems sizes and values of s = 0 is unprecedented for these three models.

Reference dynamics with truncated bond dimensions
While in the extreme active/inactive limits we can achieve a good MPS description with just a bond dimension of O(10), one may need a bond dimension of O(100) for the more difficult regions such as around s = 0 [40,43]. One reason for the necessity of this high bond dimension could be that the state has longer-ranged spatial correlations. Another could be that when one runs the vMPS, we run it against some constraint in the state space. For the FA model, this is the weak constraint that restricts to the connected component of all configuration but the one with n i = 0 for all i. For the SSEP, we have the stronger constraint that we are within the state space with fixed N p particles.
The goal is to look for a state with a smaller bond dimension than we currently have which still contains all the necessary interactions, but if necessary, discards the information which enforces the constraint. Then, by starting our CTMC algorithm in a state which satisfies the constraint, we will automatically enforce it for the rest of the trajectory, as the dynamics keeps the system in the constrained subspace.
Approximating a TN by another one with a small bond dimension is known as truncation. For MPS as we use, this can be achieved via a singular value decomposition across each bond, where only the largest D < D singular values are kept. We show this truncation in Fig. 2(a) for SSEP (as this typically requires the largest bond dimension), where we run the vMPS to at least (but higher if required) D = 50 to find the state |ψ D , and then truncate to |ψ D with the bond dimension D < D. We measure the truncation error ε = 1 − | ψ D |ψ D | 2 between the two states, where we assume both are normalised. We find that when far from the critical point, we can describe the original state to a high accuracy with bond dimensions as small as D ∼ 20. Conversely, we cannot attain the same level of accuracy for s ∼ s c , where the state exhibits larger amounts of entanglement.
There are multiple reasons that one may want to find a state with a truncated bond dimension. The first is that our Monte Carlo algorithm scales quadratically with the bond dimension -this could hinder the convergence of time-dependent observables at large times, which can require a large sample size to be determined with sufficient accuracy. For such situations, reducing the scaling of the algorithm would be desired. Another reasoning could be that we want to investigate a system which requires a higher complexity of TN, such as 2D system with projected entangled pair states (PEPS) [65,75]. Not only would the scaling of our CTMC algorithm increase, but so would the scaling of the variational algorithm used to find the reference dynamics. In this case, one may not be able to reach a bond dimension large enough to give a desirable variance.
We show the measured dynamical activity for SSEP and the East model (symbols), with a reference dynamics constructed from states with a truncated bond dimension in Fig. 2(b, c), and compare to the expected result from the non-truncated MPS (dashed line). Sur-prinsingly, we find that for the most part, even for bond dimensions as small as D = 2, we can accurately reproduce the correct dynamical activity for each of the models. As expected, the truncation struggles mostly around the transition point. Nevertheless, we can achieve good results for the FA (not shown) and East with a truncated bond dimension of D = 4, and D = 10 for SSEP.
The calculations done thus far have been with a reference dynamics constructed using a truncated bond dimension without any trajectory re-weighting. In principle, Eq. (37) is exact and thus allows for further improvements by using the umbrella g(ω) = e −tθ(s)+ dt∆R(ω) .
We implement this re-weighting via transition path sampling (TPS) with the shifting method, see Refs. [12,66] for further details. Figure 2(d) shows the results of this umbrella sampling for the FA model with a s value close to the critical point, s c . It is here the discrepancy is the largest, and we can do a more detailed analysis by looking at a larger range of bond dimensions. We see a significant improvement when using the re-weighting factor Eq. (48). It might be that we could see further improvements with more TPS iterations.
The main point to take from this is that we are able to achieve accurate results for the dynamical activity (the observable we are tilting) and some local observables with a relatively small bond dimension. This of course comes at a cost however, as when we truncate we discard some of the information that accounts for the long-ranged spatial correlations. For the case of SSEP, even though apparently we are discarding a large amount of information when truncating (cf. Fig. 2(a)), it seems that we keep the relevant information needed to reproduce the correct dynamics, but at the cost of not maintaining the conservation law. We note however that it is possible to explicitly implement the conservation laws in the MPS [76], but it is not clear how this will affect the quality of the reference dynamics in the CTMC algorithm.

B. Sampling rare events of finite times
For the previous results, we disregard finite-time effects by considering our sampled trajectories to be a "slice" of a larger infinite-time trajectory. We now look to incorporate these effects back into our sampling by considering the full re-weighting factor Note that previously, for a large bond dimension, the part of Eq. (49) which accounts for the difference in escape rate had a negligible effect, and could be ignored. This is not always the case here, as the umbrella sampling at the time-edges of the trajectory causes the system to visit configurations which are atypical in the Doob dynamics, and not well described by our MPS approximation. As a proof of principle, we start by comparing results from TPS with the original dynamics against TPS with the reference dynamics for a small system size N = 40, and a variety of times, as is shown for the FA model at s = ±0.1 in Fig. 3(a). For small times, both show excellent agreement. For large times however, the normal dynamics struggles to correctly account for the expected activity shown by the dotted lines, a result of the exponential time-dependence in Eq. (26) (asK is time extensive). While sampling with our reference dynamics reduces the exponential cost in time, the time-edges still suffer from an exponential sampling difficulty in the system size. This is most noticeable for the inactive phase, where each model exhibits an exponential localization [40,43,77] at the spatial edge(s) of the system. This causes the l ref s (x) values to exponentially vary. Nevertheless, it is still a significant improvement on the previously exponential cost in space, time and s.
The average occupations n i s (at site i) for s = 0.1 and t = 100 is shown in Fig. 3(b), while Fig. 3(c) shows the average excitation density n s = N −1 i n i s for s = −0.1 and t = 10. It is here the time-edge effects become obvious; we start at a state which lies somewhere between the expected s = 0 (dashed line) dynamics and the expected long-time dynamics, which depends on the whole spectrum of W s , as well as the total trajectory time. The system quickly evolves and resembles the Doob dynamics. Note that at the end of trajectory, it is again described by the original probability vector, as is expected due to the time-symmetry in Eq. (37).
Finally, Fig. 3(d) shows the average dynamical activity as a function of s and time, t. We show the expected activity in the infinite time limit t = ∞ as a black dashed line, and the measured activity for finite times as symbols. Notice that as time decreases, the drop in activity becomes less sharp. Furthermore, the transition from the active to inactive phase happens at decreasing s. While the methods presented here could allow for a detailed investigation into the temporal scaling of the critical point [52,53,71,72], doing so for desirable system sizes would be at a large computational cost. We hope to investigate this more extensively using time-evolution methods (see e.g. [30,65,78]).

V. CONCLUSIONS
We have expanded on previous applications of TNs to classical constrained models [40][41][42][43], using the MPS approximation of the leading eigenstates of a tilted stochastic generator in 1D to construct a reference dynamics which well approximates the exact Doob dynamics. This allows us to (nearly) optimally sample the rare events of 1D constrained systems with just a polynomial cost in both space and time -rather than the exponential cost of most sampling methods. We have demonstrated here the efficiency of this approach by generating tilted trajectory ensembles for the FA and East KCMs and the symmetric simple exclusion process. Our simulations are for sizes and times unprecedented for such large deviation studies.
Furthermore, our results show that it is possible to obtain an accurate dynamics away from the dynamical transitions of the models we studied with a truncated bond dimension, which enables close to optimal sampling simulations at little cost. Further extensions of our work includes generalising our methods to higher dimensions, for example by using two-dimensional variational tensor network techniques, such as PEPS [65,75] to approximate the leading eigenvectors of 2D classical generators, as is done in Ref. [42]. From the associated leading eigenvectors, as we have shown here, we can in turn construct a reference dynamics which is nearly optimal for sampling rare trajectories. While PEPS algorithms do not currently allow for bond dimensions comparable to vMPS, they remain a fruitful area of research which is constantly being improved on [79][80][81][82][83][84][85][86][87]. Recent works [88] have shown the effectiveness of using recurrent neural networks (RNN) to approximate the leading eigenstates of tilted generators in two dimensions. The methods presented here can be generalized to RNN to allow for the efficient sampling of 2D rare events.
Another area that deserves exploration is to apply similar TN methods to systems which do not obey detailed balance, and for which their generators cannot be brought to a Hermitian form. While this would damper the effectiveness of variational algorithms, approaches based on time evolution may offer a promising solution (see e.g. [30,65,78]). Such approaches could also offer further insights into intermediate time rare events, where both usual sampling methods and large deviation approaches fall short. We hope to report on such studies in the near future.