{'fraction_non_alphanumeric': 0.044409707095815286, 'fraction_numerical': 0.013632985675114676, 'mean_word_length': 4.174224343675418, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Backpropagating Hybrid Monte Carlo algorithm for fast Lefschetz thimble calculations

The Picard-Lefschetz theory has been attracting much attention as a tool to evaluate a multi-variable integral with a complex weight, which appears in various important problems in theoretical physics. The idea is to deform the integration contour based on Cauchy's theorem using the so-called gradient flow equation. In this paper, we propose a fast Hybrid Monte Carlo algorithm for evaluating the integral, where we"backpropagate"the force of the fictitious Hamilton dynamics on the deformed contour to that on the original contour, thereby reducing the required computational cost by a factor of the system size. Our algorithm can be readily extended to the case in which one integrates over the flow time in order to solve not only the sign problem but also the ergodicity problem that occurs when there are more than one thimbles contributing to the integral. This enables, in particular, efficient identification of all the dominant saddle points and the associated thimbles. We test our algorithm by calculating the real-time evolution of the wave function using the path integral formalism.


Introduction
Many problems in theoretical physics are formulated in terms of a multi-variable integral with some weight, where the number of variables or the system size is typically very large. If the weight is positive semi-definite, we can regard it as a probability distribution and perform Monte Carlo calculation based on the idea of importance sampling. However, it is not straightforward to extend this method to the case with a complex weight. For instance, if we use the absolute value of the complex weight as the probability distribution for sampling and take into account the phase factor by reweighting, huge cancellation occurs among the sampled configurations. For this reason, one needs exponentially large numbers of configurations as the system size increases in order to obtain results with sufficient precision. This is the sign problem, which has been hindering the development of theoretical physics in various branches.
The situation regarding the sign problem has changed drastically in the last decade, however. There is by now a bunch of new techniques that have proven to work beautifully in certain classes of problems. For instance, the tensor renormalization group [1,2,3,4,5] is applicable if one can reformulate the problem in terms of a network of tensors, where the indices of tensors are contracted with those of the adjacent tensors. This method does not rely on the importance sampling, and hence the sign problem is absent from the outset.
Another promising direction is to complexify the integration variables keeping the holomorphicity of the weight and the observables. There are roughly two methods that belong to this category. One of them is the complex Langevin method (CLM) [6,7,8,9,10,11], which can be viewed as a naive extension of the ordinary Langevin method for positive semi-definite weights. In this method, the expectation values of observables can be obtained with configurations generated by solving the complex Langevin equation, which involves the Gaussian noise term and the drift term representing the derivative of the weight. The computational cost is comparable to that of ordinary Monte Carlo methods for positive semi-definite weights. The validity of this method requires, however, that the probability distribution of the drift term and that of the observables should fall off sufficiently fast at large values, which is not always satisfied [10,11]. While the method is definitely worth trying for its simplicity and its capability for large system size (See, for instance, Refs. [12,13,14,15,16,17] for semi-realistic calculations.), one needs to have the luck of satisfying the validity condition in order to explore the most interesting parameter regime.
The other method based on complexification of integration variables uses Cauchy's theorem and deforms the integration contour in such a way that the integrand does not oscillate much in phase. One can then apply standard Monte Carlo methods to perform algorithms solve the Hamilton equation of motion on the deformed contour [24,29] or on the worldvolume [30], which we shall refer to as the deformed manifold collectively in what follows. This requires complicated procedures of determining a constrained motion on the deformed manifold, which is unknown a priori. In addition to the force obtained by taking the derivative of the modulus of the complex weight, one has the normal force characteristic to constrained systems, which has to be determined in such a way that the complexified configuration is constrained on the deformed manifold. This is done by Newton's method, which requires solving the flow equation at each iteration.
Here we consider solving the Hamilton equation of motion on the original contour, which saves us from all the burden associated with treating the constrained motion. The main task is transferred to the calculation of the force, which requires differentiating the modulus of the complex weight evaluated on the deformed manifold, now with respect to the configuration on the original contour 3 . At first sight, one might think that this requires the calculation of the Jacobian matrix associated with the deformation of the integration contour 4 . The crucial point of our proposal is that one can actually avoid this calculation by "backpropagating" the force on the deformed contour to that on the original contour, which reduces the computational cost by a factor of the system size.
The only drawback of our HMC algorithm compared with the existing ones is that the modulus of the Jacobian is not included in sampling and hence should be taken into account by reweighting. Note, however, that the phase of the Jacobian has to be taken into account by reweighting anyway in the GTM in general. The calculation of the Jacobian is time-consuming, but it can be done off-line since it is needed only when one measures the observables for statistically decorrelated configurations obtained by the HMC algorithm. If one uses a crude estimator of the Jacobian [33], the computational cost of the GTM becomes comparable to ordinary Monte Carlo methods just like the CLM. Thus our algorithm opens up a new possibility for applying the GTM to large systems.
The rest of this paper is organized as follows. In Section 2, we briefly review the GTM. In Section 3, we present our HMC algorithm on the original contour. In Section 4, we extend our algorithm to the case in which one integrates over the flow time to avoid the ergodicity problem. In Section 5, we clarify the relationship to the existing HMC algorithms. In Section 6, we test our algorithm by calculating the real-time evolution of the wave function using the path integral formalism. Section 7 is devoted to a summary and discussions. In Appendix A, we provide a simple understanding for the reduction of computational cost by backpropagation. In Appendices B and C, we provide some notes related to Section 5.

Brief review of the GTM
In this section, we briefly review the GTM for a general model given by the partition function In the GTM [25], we deform the integration contour into C N by using the so-called holomorphic gradient flow equation which is solved from σ = 0 to σ = τ with the initial condition z(x, 0) = x ∈ R N . The flowed configurations define an N -dimensional real manifold embedded in C N , which we denote as M τ = {z(x, τ )|x ∈ R N }. According to Cauchy's theorem, the integration contour can be deformed continuously from R N to M τ ⊂ C N without changing the partition function. Then, by noting that z = z(x, τ ) defines a one-to-one map from x ∈ R N to z ∈ M τ , one can rewrite the partition function as and similarly the observable (2.2) as where the Jacobian matrix J(x, τ ) corresponding to the map z(x, τ ) is defined by Taking the derivative with respect to x l on each side of (2.3), one obtains the differential equation for J(x, σ) as which can be solved from σ = 0 to σ = τ with the initial condition J kl (x, 0) = δ kl to obtain the Jacobian matrix J(x, τ ). Here we have defined the Hessian of the action as which plays an important role in this paper. Note that the Hessian is a sparse matrix with only O(N ) nonzero elements if the action is local, meaning that only terms with local couplings among z i 's exist as in the case of local field theories. The virtue of using the holomorphic gradient flow (2.3) in deforming the integration contour can be understood by taking the derivative of the action with respect to the flow time as which implies that the real part of the action S grows monotonically along the flow while keeping the imaginary part constant. Thus, in the σ → ∞ limit, ReS(z(x, σ)) diverges except for such points x = x that satisfy lim σ→∞ ∂ ∂z k S(z(x , σ)) = 0, which are the points that flow into some fixed points z = z defined by ∂S(z ) ∂z k = 0. The deformed contour M τ at τ = ∞ consists of a set of Lefschetz thimbles associated with these fixed points, which correspond to the infinitesimal vicinities of the points x on the original contour. Since ImS(z(x)) is constant on each thimble, the sign problem is maximally reduced. The original proposal was to perform Monte Carlo integration on these thimbles [22,23,24].
However, when there are more than one thimbles, the transition between thimbles does not occur during the Monte Carlo simulation, which leads to the violation of ergodicity. The GTM [25] avoids this problem by choosing a sufficiently small flow time allowing the ImS(z(x)) to fluctuate to some extent, which works as far as the system size is not so large. For a larger system size, one can integrate over the flow time in order to avoid the sign problem and the ergodicity problem at the same time [30].

HMC algorithm with backpropagation
In this section, we discuss our HMC algorithm on the original contour. The crucial point is to calculate the force in the fictitious Hamilton dynamics by using the idea of backpropagation, which reduces the computational cost by a factor of O(N ). Here we consider the case of fixed τ and defer the discussion in the case of integrating τ to Section 4.

HMC algorithm on the original contour
In order to simulate the partition function (2.4), let us consider the partition function where the imaginary part of the action S as well as the Jacobian is omitted. The expectation value (2.5) can be obtained by the standard reweighting formula as where · · · τ represents the expectation value with respect to the partition function (3.1). The most naive way to simulate the partition function (3.1) is the Metropolis algorithm. Starting from some configuration x ∈ R N , one constructs a trial configuration x = x + δx ∈ R N with certain probability distribution for δx. Then one calculates the change of the action and accepts the trial configuration with the probability min(1, e −δS ). Both the numerator and the denominator in (3.2) can be obtained by taking an average over the configurations generated by the above algorithm. While this algorithm is very easy to implement, it is not efficient since the update by δx has to be local in order to keep the acceptance rate reasonably high. The basic idea of the HMC algorithm [34] is to enable a global update δx in the Metropolis algorithm by using the information of the action in its proposal. For that, we introduce auxiliary variables p i corresponding to x i and consider the partition functioñ which is equivalent to (3.1). One can update p i by just generating Gaussian random numbers. After that, using the initial condition for x(0) and p(0) set by the current configuration, one solves the fictitious Hamilton equation where the force F i (s) is defined by for a fixed time s f to obtain x(s f ) and p(s f ), which provide the trial configuration subject to the Metropolis accept/reject procedure with the probability min(1, e −δH ), where If one can solve (3.6) exactly, one obtains δH = 0 due to the Hamiltonian conservation, which implies that the acceptance rate is 100%. In practice, one has to discretize the Hamilton equation, which makes δH = 0. The stepsize ∆s has to be chosen small enough to keep the acceptance rate reasonably high 5 . Note that the detailed balance is satisfied even at finite ∆s thanks to the Metropolis procedure, and hence there is no systematic error due to the discretization of the Hamilton equation. See Section 6.2 for more detail.

Calculating the force by backpropagation
Thus the question boils down to the calculation of (3.7). Let us first note that (3.7) can be rewritten as where J ij (x, τ ) is the Jacobian matrix defined by (2.6) and f i (s) is the gradient of the action on the deformed contour at z(x(s), τ ) corresponding to the point x(s) on the original contour. The calculation of the Jacobian matrix, however, requires us to solve the differential equation (2.7), which involves O(N 2 ) or O(N 3 ) arithmetic operations depending on whether the action S(z) is local or non-local in the sense explained below (2.8). Alternatively, one might think of evaluating the right-hand side of (3.7) by numerical derivative as in Ref. [32]. The required cost, however, is comparable to calculating the Jacobian, and what is worse, there is some systematic error associated with the numerical derivative.
Suppose we were to calculate J ij (x, τ ) v j with v j ∈ R, which appears many times in the HMC algorithm on the deformed manifold [29,30]. In this case, one can avoid the Coming back to our problem, we can actually save a factor of O(N ) as well but in a slightly different way. Let us note that the problem of calculating (3.7) is mathematically equivalent to the calculation in machine learning, where ReS(z(x, τ )) represents the loss function and z i (x, τ ) represents the parameters of the network at a layer labeled with τ . The quantity (3.7) we have to calculate is the derivative of the loss function with respect to x i = z i (x, 0), which represent the parameters at the first layer. The crucial idea here is the use of backpropagation. 6 Let us define the "force" at flow time σ by Applying the chain rule, we obtain This implies that f i (σ) obeys the differential equation which is analogous to (3.11). By solving this differential equation backwards in σ from σ = τ to σ = 0 with the initial condition f i (τ ) = f i , one obtains the force at the first layer Using this, the desired force (3.7) is obtained as Thus the computational cost is reduced by a factor of O(N ). 6 The use of backpropagation in calculating the force in HMC algorithms is actually quite common in lattice gauge theory when the action S is written in terms of smeared field variables, which are obtained by solving some gradient flow equation starting from the original field variables. We would like to thank Akio Tomiya for bringing our attention to their work [35] in this regard.

Discretizing the flow equation
In actual calculations, one has to discretize the holomorphic gradient flow equation (2.3).
In this section, we explain how our idea for calculating the force in the fictitious Hamilton dynamics can be extended to this case. In particular, we show that this can be done in such a way that the discretization causes neither systematic errors nor any decrease in the acceptance rate, which is not the case for the existing HMC algorithms [29,30] as we discuss in Section 5. Let us define σ n = n ε (n = 0, 1, · · · , N τ ) with τ = N τ ε and consider the holomorphic gradient flow (2.3) discretized as with the initial condition z(x, 0) = x ∈ R N . Taking the derivative with respect to x l on each side of (3.18), one obtains By solving this difference equation with the initial condition J kl (x, 0) = δ kl , one obtains exactly the Jacobian matrix (2.6) defined for the solution to the discretized flow equation (3.18). This means that while finite ε affects the deformed contour, thanks to Cauchy's theorem, it does not affect the equations (2.4) and (2.5), and hence does not cause any systematic error. Let us then discuss how the calculation of the force (3.7) should be extended to the case of the discretized flow equation. If one uses (3.9), the Jacobian matrix has to be obtained by solving (3.19), which is time consuming. Instead, we define the "force" at flow time σ n by Applying the chain rule, we obtain By solving this difference equation backwards in σ from σ = τ to σ = 0 with the initial condition f i (τ ) = f i , one obtains the force at the first layer (3.16). Using this, the desired force is obtained as (3.17). The computational cost is reduced by a factor of O(N ) similarly to the case of the continuous flow equation. In Appendix A, we provide a simple understanding for this huge reduction of computational cost.

Integrating over the flow time τ
In order to solve the sign problem, the flow time τ has to be sufficiently large. However, when one has more than one Lefschetz thimbles at τ = ∞, Monte Carlo simulations based on the HMC algorithm suffer from the ergodicity problem. This occurs because the regions that contribute to the integral are either disconnected or separated by singular points, at which the real part of the action diverges. In order to solve the sign problem without suffering from this ergodicity problem, it is useful to integrate over the flow time. This is the idea of the "worldvolume approach" [30], which was put forward in the context of the HMC algorithm on the deformed contour. Here we adapt this idea to our HMC algorithm on the original contour, which not only turns out to be much simpler but also solves a possible problem in the original proposal.

The basic idea
The key observation here is that in the reweighting formula (3.2), the ratio is τ -independent although both the numerator and the denominator on the right-hand side depend on τ . This implies that one can actually integrate over τ and consider where W (τ ) is a real function of τ that can be chosen so that the distribution of τ becomes almost uniform within a certain range [30]. The expectation value of an observable can be calculated as where · · · W represents the expectation value with respect to the partition function (4.1).
We have introduced an arbitrary function W (τ ), which can be chosen to minimize the statistical errors [36].

Extending our HMC algorithm
Let us extend our HMC algorithm to the model (4.1). As we did in the previous section, we introduce auxiliary variables (p i , p τ ) corresponding to (x i , τ ) and consider the partition functionZ W = dτ dp τ dx dp e −H , (4.3) The Hamilton equation for x i (s) and p i (s) are given by (3.6) as before, but here we also have the equations for τ (s) and p τ (s) given by where the force F τ (s) for τ is defined by .   where z k (x, N τ ε) is obtained by solving the difference equation (3.18). Taking the εderivative on each side of (3.18), where σ n (ε) = n ε is treated now as a function of ε, we obtain Solving this difference equation with the initial condition one can obtain the first factor in (4.10). The computational cost for this procedure is comparable to that for calculating the force in the x-direction.
Note that this is an exact calculation of the force for finite ε, and hence there is no systematic error here, either. This is a significant advantage of our algorithm over the existing ones, in which the discretization of the flow equation causes some systematic error as we discuss below.

Relationship to the existing HMC algorithms
As we mentioned in the Introduction, there are HMC algorithms proposed for the GTM in the past. The difference from our HMC algorithm is that the existing ones [29,30] deal with the Hamilton dynamics of the point z(x, τ ) after the flow, whereas we deal with the Hamilton dynamics of the point (x, τ ) before the flow. Since the two points are in one-toone correspondence, we should be able to compare the Hamilton dynamics. In this section, we clarify the relationship between the two types of algorithms from this point of view. Some basic properties of the deformed manifold necessary to understand this section are recapitulated in Appendix B for the readers' convenience.
Here we will discuss the case of integrating τ . For that, in this Section alone, we introduce the notation x µ and p µ (µ = 0, 1, · · · , N ) for the dynamical variables and the corresponding momenta, where we have defined x 0 ≡ τ and p 0 ≡ p τ . Accordingly, we will also use z i (x) instead of z i (x, τ ) for the flowed configuration. The flow equation can be assumed to be either discretized or continuous, which does not matter in our discussion below. The above notations enable us to obtain the results for the case of a fixed flow time τ by simply replacing the index µ = 0, · · · , N by i = 1, · · · , N .
The Hamilton dynamics of the point z i (x) is described by the Hamiltonian where π i is the momentum on the deformed manifold. Therefore, the Hamilton equation on the deformed manifold is given by where N i is the normal force, which is perpendicular to the tangent space of the deformed manifold, and hence does not affect the Hamiltonian conservation. The normal force N i is determined in such a way that the momentum π i (s) resides in the tangent space of the deformed manifold at the point z i (s), which is constrained to be on the deformed manifold. In fact, the Hamilton dynamics on the deformed manifold given above is not equivalent to the Hamilton dynamics on the original contour that we have considered. Instead, it is described by the Hamiltonian on the original contour that takes the form which involves a nontrivial kernel K µν (x) = K νµ (x) in the kinetic term to be specified later. The Hamilton equation therefore takes the form The s-evolution of the flowed configuration z i (x(s)) is then given by Comparing this with (5.2), we obtain the relationship which shows that the momentum π i on the deformed manifold resides in the tangent space as it should. Plugging (5.8) and z i = z i (x) into (5.1), one obtains where we have defined Since (5.9) should be identified with (5.4), we obtain an identity 1 2 K(x) = K(x) g(x) K(x), which implies Next, let us consider the s-evolution of the momentum π i (x) on the deformed manifold. From (5.8), we find where we have defined 14) The first term F i in (5.12) is given by This is nothing but the gradient force − ∂V (z,z) ∂z i projected onto the tangent space. The second term N i in (5.12) can be interpreted as the normal force that is needed to realize the constrained motion. This can be checked by confirming that its projection onto the tangent space vanishes identically; namely, See Appendix B for the details. The appearance of the kernel K = 1 2 g −1 in (5.4) can be understood also from the viewpoint of what we are simulating. The Hamilton dynamics we consider in Section 4.2 is intended to reproduce the partition function (4.1). This can be seen by integrating out the momentum variables in (4.3), which yields (4.1). If we have a nontrivial kernel K = 1 2 g −1 as in (5.4), the integration over the momentum p µ in (5.4) yields an extra factor det g(x) in the partition function. This is nothing but the volume element of the deformed manifold since g(x) defined by (5.10) represents the induced metric on the deformed manifold, which is embedded in C N . Thus our result (5.4) is consistent with the fact that the HMC algorithm using the Hamilton dynamics on the deformed manifold includes the volume element of the deformed manifold in the partition function [29,30].
In the case of fixed τ , the induced metric (5.10) is given by g kl = Re(J † J) kl for the discretized flow equation. In the continuum limit ε → 0 of the flow equation, we have Im(J † J) = 0, which implies that det g(x) = | det J(x)|. (See Appendix C for the details.) Thus, neglecting the finite ε effects, the modulus of the Jacobian is included in sampling, which may have certain advantage, but the price one has to pay seems overwhelming.
In the case of integrating τ , on the other hand, one obtains det g(x) = C(x) | det J(x)| neglecting the finite ε effects, where C(x) ≥ 0 needs to be calculated and taken into account in reweighting [30]. Moreover, if the fixed point of the flow equation that contributes to the integral resides on the original contour, the volume element det g(x) of the worldvolume vanishes at that point, which causes an ergodicity problem. 7 In fact, this problem occurs because the induced metric g µν (x) on the worldvolume becomes singular at that point. Obviously, our HMC algorithm based on the Hamilton dynamics on the original contour is totally free from such problems.
Note also that the HMC algorithms on the deformed manifold in Refs. [29,30] assume Im(J † J) = 0, which is violated upon discretization of the flow equation as we show in Appendix C. There are actually two places in which this assumption is used. One of them is, as already mentioned above, in the statement that the HMC algorithm on the deformed contour includes | det J| in sampling. The other is in constructing the basis vectors of C N orthogonal to the tangent space, which are used in solving the constrained motion on the deformed manifold. Therefore, the HMC algorithms in Refs. [29,30] suffer from systematic errors due to discretization of the flow equation. In contrast, nowhere in our HMC algorithm on the original contour have we used Im(J † J) = 0. Hence, our algorithm is free from systematic errors whatsoever.

Practical applications
In this section, we consider an important improvement of our algorithm suggested also from our discussion in Section 5. This makes the discretization of the Hamilton equation slightly more nontrivial as we describe below. We test our algorithm by applying it to the time evolution of the wave function in the path integral formalism.

Introducing a mass parameter in the HMC
As we have seen in Section 5, the Hamilton dynamics on the deformed manifold is equivalent to introducing a nontrivial kernel K(x) = 1 2 g −1 (x) in (5.4) in the Hamilton dynamics on the original contour 8 . In a similar spirit, for practical applications, we find it very important to generalize our HMC algorithm by introducing a τ -dependent mass parameter m(τ ) as Z W = dτ dp τ dx dp e −H , (6.1) Note first that integration over the momentum gives a factor m(τ ) N/2 in the partition function. In the case of fixed τ , this is just a constant factor, while in the case of integrating τ , it can be absorbed into the definition of W (τ ) in (4.1). This also suggests that m(τ ) should be chosen to be proportional to the typical value of | det J(x, τ )| 2/N for various x with fixed τ . The need for this improvement can be understood from the fact that the flow equation typically maps a small region on the original contour to an exponentially large region on the deformed contour at large τ , where the scale factor is given by | det J(x, τ )| 1/N . Introducing the mass in the HMC changes the effective stepsize in solving the Hamilton equation on the original contour by the factor of 1/ m(τ ). Therefore, the above choice of m(τ ) enables a random walk on the deformed manifold with almost uniform discretization.
From the Hamiltonian (6.2), the Hamilton equation is obtained as where the forces F i (s) and F τ (s) are defined by (3.7) and (4.7), respectively, as before.

Discretizing the Hamilton equation in the HMC
In actual calculations, the Hamilton equation in the HMC algorithm has to be discretized. In the case of fixed τ , the Hamiltonian takes the canonical form H = 1 2 p 2 + V (x) so that one can use the standard leap-frog discretization. This respects the reversibility and the preservation of the phase space volume [34], which guarantees the detailed balance. In the case of integrating τ , one should note that the first term of the Hamiltonian (6.2) mixes the momentum variables p i with τ , which is now treated as one of the coordinate variables. For this reason, we have to generalize the leap-frog discretization slightly. Let us note first that the Hamilton equation takes the form dp τ ds = A(x, τ, p) , dx ds = B(τ, p) , dτ ds = C(p τ ) , dp ds = D(x, τ ) . (6.7) We discretize this equation as A(x(s n ), τ (s n ), p(s n )) , (6.8) x(s n+1/2 ) = x(s n ) + ∆s 2 B(τ (s n ), p(s n )) , (6.9) τ (s n+1/2 ) = τ (s n ) + ∆s 2 C(p τ (s n+1/2 )) , (6.10) p(s n+1 ) = p(s n ) + ∆s D(x(s n+1/2 ), τ (s n+1/2 )) , (6.11) τ (s n+1 ) = τ (s n+1/2 ) + ∆s 2 C(p τ (s n+1/2 )) , (6.12) x(s n+1 ) = x(s n+1/2 ) + ∆s 2 B(τ (s n+1 ), p(s n+1 )) , (6.13) p τ (s n+1 ) = p τ (s n+1/2 ) + ∆s 2 A(x(s n+1 ), τ (s n+1 ), p(s n+1 )) , (6.14) where s ν = ν∆s with ν being an integer or a half integer. We repeat the above procedure for n = 0, 1, · · · , (N s − 1), where s f = N s ∆s represents the total fictitious time for the Hamilton evolution. One can prove the reversibility and the preservation of the phase space volume as in the standard leap-frog discretization. Note that the procedure (6.14) for n = k and the procedure (6.8) for n = k + 1 can be combined into one step as for k = 0, 1, · · · , (N s − 2). In actual simulation, it is useful to specify a finite range of τ over which we integrate in (4.1). If we realize this by choosing W (τ ) in (4.1) to be large outside the range of τ , the force in the Hamilton dynamics becomes very large when τ gets out of the range. Instead, we introduce walls at both ends of the range. At the discretized level, when τ gets out of the range at (6.10), we go back to (6.8), flip the sign of p τ (s n ) there and continue. One can easily prove that introducing this procedure does not violate the reversibility of the Hamilton dynamics.

Time evolution of the wave function by the path integral
In this section, we test our algorithm by applying it to calculations in quantum mechanics using the path integral formalism (See, for instance, Refs. [37,38,39] for calculations of the real-time correlator.). The fundamental object in this formalism is the transition amplitude where the integral is taken over all the paths x(t) (t i ≤ t ≤ t f ) with the constraints x(t i ) = x i and x(t f ) = x f . The time evolution of the wave function is given by where the path integral is taken now with the constraint x(t f ) = x f only.
Here we calculate the time-evolved wave function (6.17) numerically for the action where we use m = 1 and the quartic potential V (x) = 1 4! x 4 . We assume, for simplicity, a Gaussian form for the initial wave function where we set γ = 1 and α = 1. Let us first discretize the time t as t n = (n − 1) and define x n = x(t n ), where n = 1, · · · , (N + 1). We also define t i = t 1 , x i = x 1 and t f = t N +1 , x f = x N +1 . Thus the time-evolved wave function Ψ(x f , t f ) can be obtained by evaluating the partition function with the dynamical variables x i (i = 1, 2, · · · , N ), where the action S(x; x f ) is given by We consider the time evolution for t f = 2, which is discretized by N = 9. The flow equation is discretized with N τ = 10 steps. In the HMC algorithm, we use N s = 15 steps with the total fictitious time s f = 0.25. We have confirmed that the computational cost for making one step in the HMC algorithm scales linearly with N .
In Fig. 1 (Left), we plot 2 N log | det J(z(x, τ ))| obtained for each configuration against τ . We can fit the data points to a linear function c τ with c = 11.795, which is consistent with the fact that the typical scale factor between the original and deformed contours grows exponentially with the flow time τ . According to our discussion in Section 6.1, we therefore use m(τ ) = e c τ in the Hamiltonian (6.2). In reweighting (4.2), we choose W (τ ) = 1 2 N c τ in order to cancel the reweighting factor | det J(x, τ )| ∼ m(τ ) N/2 . In Fig. 1 (Right), we plot the phase distribution of the reweighting factor (4.2) for configurations obtained within various regions of τ . We find that the peak of the distribution becomes sharper as one goes to larger τ , which confirms that the sign problem becomes weaker as τ increases.
Instead of calculating the time-evolved wave function (6.17) directly, we calculate its log derivative by where the expectation value is taken with respect to the partition function (6.20). The results for N = 9 are shown in Fig. 2, which is in good agreement with the results obtained by diagonalizing the Hamiltonian that corresponds to a continuous time evolution (N = ∞).

Summary and discussions
In this paper, we have proposed a new HMC algorithm for fast Lefschetz thimble calculations. Unlike the existing HMC algorithms, which solve the Hamilton dynamics on the deformed contour, we solve the Hamilton dynamics on the original contour using the action evaluated at the point obtained by solving the holomorphic gradient flow equation. The crucial point of our proposal is to calculate the force in the Hamilton dynamics by backpropagating the force on the deformed contour to that on the original contour, analogously to the well-known idea that plays a central role in machine learning. The computational cost in the present calculation is reduced by a factor of the system size.
The possible ergodicity problem can be avoided by integrating over the flow time [30], which can be naturally implemented in our algorithm. We discussed, in particular, that the flow equation can be discretized even in this case without causing systematic errors unlike the existing HMC algorithms. We have also discussed the relationship to the existing HMC algorithms from the viewpoint of the Hamilton dynamics, which revealed the existence of a nontrivial kernel in the kinetic term. This kernel actually reproduces the volume element of the deformed manifold upon integrating the auxiliary momentum variables. Inspired by this observation, we introduced a flow-time dependent mass parameter in the Hamilton dynamics, which turns out to be important in practical applications.
Our algorithm is particularly useful in identifying the saddle points and the associated thimbles that contribute to the integral in question. This may shed light, for instance, on quantum tunneling phenomena from the viewpoint of the real-time path integral [40,41,42,43,44]. On the other hand, when one calculates the expectation values, one needs to calculate the Jacobian that appears in the reweighting factor. This is time-consuming, but it can be done off-line only when one measures the observables. If one uses a crude estimator of the Jacobian whose cost is only proportional to the system size [33], the computational cost becomes comparable to ordinary Monte Carlo simulation for positive semi-definite weights. We hope to apply our algorithm to a large system in that way, although it would be certainly nice if there is a better way to treat the Jacobian.
Last but not the least, the "original contour" in our algorithm does not have to be the real axis but can be deformed in the complex plane as far as one does not pass through singularities. This points to the possibility of combining the GTM with various proposals for path optimization [45,46,47,48,49,50,51,52], which may be useful in reducing the flow time required to solve the sign problem, and hence in reducing the effects of the Jacobian in reweighting.
The solution to the differential equation (A.1) can be written down formally as where 1 N is the N × N unit matrix, and P exp represents the path-ordered exponential, which ensures that H(σ) with smaller σ comes on the right after Taylor expansion. On the other hand, F j defined by (3.9) can be written as A naive method would be to obtain the Jacobian by (A.3) and to obtain F j by (A.4), which corresponds to taking the products in (A.5) from the right. Let us consider how the backpropagation calculates the same quantity (A.5). First, the "force" f i (σ) at flow time σ defined by (3.12) can be written as This quantity satisfies the differential equation which is equivalent to (3.15). By solving this backwards in σ from σ = τ to σ = 0 with the initial condition f j (τ ) = f j , we obtain Using this in (A.5), we obtain Thus, the backpropagation amounts to reversing the order of multiplications in (A.5), which replaces the matrix-matrix products by the matrix-vector products, thereby reducing the computational cost by a factor of O(N ).
The point just mentioned becomes clearer in the case of discretized flow time discussed in Section 3.3. Corresponding to (A.1), we can rewrite (3.19) in the form where the 2N × 2N matrix H(σ) is defined by (A.2) as before. The formal solution corresponding to (A.3) can be written as where the product is taken in such a way that a factor with smaller k comes on the right. Plugging (A.11) in (A.4), we obtain Corresponding to (A.6), let us define which represents the force propagated backwards in σ. Note that this quantity satisfies the difference equation which is equivalent to (3.22). By solving this difference equation with the initial condition f j (τ ) = f j , we obtain the desired force as (3.17). Thus, in the discretized version, the backpropagation simply amounts to taking the products in (A.12) from the left, which reduces the cost by a factor of O(N ) compared with taking the products from the right.

B Some basic properties of the deformed manifold
In this appendix, we briefly review some basic properties of the deformed manifold [30], which is necessary in understanding Section 5. We also provide a proof for the statement (5.17).
The deformed manifold consists of points z i (x), which can be obtained by solving either the continuous or discretized version of the holomorphic gradient flow equation for the flow time τ = x 0 with the initial point given by x i . Therefore, the deformed manifold is parametrized by x µ (µ = 0, 1, · · · , N ), and it is embedded in C N . As we mentioned in Section 5, our discussion applies to the case of fixed τ as well by replacing the index µ = 0, 1, · · · , N with i = 1, · · · , N .
The tangent space at each point on the deformed manifold is a real linear space spanned by the basis vectors and any element of it can be represented as Since the deformed manifold is embedded in C N , one can define the inner product of two tangent vectors u i and v i as Expanding the tangent vector v i as (B.2) and similarly for u i as we can rewrite the inner product u, v in terms of the coefficients b µ and c µ as where g µν is defined by which gives the induced metric (5.10) on the deformed manifold. Let us consider a vector w ∈ C N and project it onto the tangent space. For that, let us we decompose w as where a µ ∈ R is specified later, and y i is a vector orthogonal to the tangent space; i.e., Re E i µ y i = 0 for all µ . Plugging this in (B.7), the projection of w onto the tangent space is given by This confirms that (5.16) is indeed the projection of the gradient force − ∂V (z,z) ∂z i onto the tangent space.