Preconditioned flow as a solution to the hierarchical growth problem in the generalized Lefschetz thimble method
Jun Nishimura Katsuta Sakai Atis Yosprakob





 


  KEK-TH-2618

 









  Preconditioned flow as a solution to the hierarchical growth problem
 in the generalized Lefschetz thimble method

  
  
  
  















         Jun Nishimura^1,2)[ E-mail address : jnishi@post.kek.jp],

         Katsuta Sakai^1,3)[ E-mail address : sakai.las@tmd.ac.jp] and
         Atis Yosprakob^2,4)[ E-mail address : ayosp@phys.sc.niigata-u.ac.jp]








^1)KEK Theory Center,
Institute of Particle and Nuclear Studies,

High Energy Accelerator Research Organization,

  1-1 Oho, Tsukuba, Ibaraki 305-0801, Japan

^2)Graduate Institute for Advanced Studies, SOKENDAI,

1-1 Oho, Tsukuba, Ibaraki 305-0801, Japan 


^3)College of Liberal Arts and Sciences, 
Tokyo Medical and Dental University,

2-8-30 Kounodai, Ichikawa, Chiba 272-0827, Japan


^4)Department of Physics, Niigata University,

8050 Igarashi 2-no-cho, Nishi-ku, Niigata-shi, Niigata 950-2181, Japan







  The generalized Lefschetz thimble method is a promising approach that
attempts to solve the sign problem
in Monte Carlo methods

by deforming the integration contour

using the flow equation.


Here we point out 

a general problem that occurs

due to the property of the flow equation,
which extends a region on the original contour
exponentially to a region on the deformed contour. 


Since the growth rate for each eigenmode is governed 
by the singular values of the Hessian of the action,

a huge hierarchy in the singular value spectrum,
which typically appears for large systems,
leads to various technical problems in numerical simulations.



We solve this hierarchical growth problem by preconditioning the flow


so that the growth rate becomes identical for every eigenmode.






As an example, we show that the preconditioned flow enables us
to investigate the real-time quantum evolution of an 
anharmonic oscillator with the system size 
that can hardly be achieved by using the original flow.




 



























§ INTRODUCTION



In theoretical physics, we are often faced with


a system which cannot be solved analytically.

In that case, we have to either perform perturbative calculations
by taking some limits in the parameter space, 
or make some plausible simplifications so that
the system becomes analytically tractable.

Such approximations, however,
quite often obscure the physics we are most interested in.




In this regard,
numerical simulation is of particular importance

since it provides us with a powerful tool to investigate various systems
nonperturbatively from first principles without such approximations.


Here we are concerned with numerical simulation of multi-variable integrals 
with some weight, where the number of variables N
representing the system size is large. 
In the case with a positive-definite weight
such as the Boltzmann weight that

appears
in classical statistical mechanics,

the numerical integration

can be realized by
Monte Carlo (MC) simulation
based on the importance sampling,
where we interpret the weight as the probability distribution.
However, we often encounter


a weight w=e^-S with a complex action S, 
which cannot be interpreted as the probability distribution.
A naive prescription is the so-called reweighting,
where we interpret the absolute value of the weight as the probability
distribution,
and take into account the phase factor when we take the average of the observables.
This does not really work for large systems since the phase factor

oscillates violently depending on the generated configurations, which
leads to huge cancellations among them.
Thus,

in order to obtain the correct expectation value,


one needs a huge number of configurations,
which typically grows exponentially with the system size N. 
This is the notorious sign problem, which has been hindering

nonperturbative studies of various important systems,
such as finite density QCD, theories with a θ-term,
strongly correlated electron systems,

and the Lorentzian path integrals for the real-time quantum evolution.

In the last decade, various new approaches have been developed to overcome 
the sign problem.

In particular, within the framework of MC methods,

two promising approaches,

the complex Langevin method <cit.>
and the Lefschetz thimble
method <cit.>
have been studied intensively[
As a promising approach without MC simulation, 
the tensor renormalization group has been developed <cit.>. 
The basic idea is to rewrite the integral
as a network of tensors and to perform a 
coarse-graining procedure iteratively using the singular value decomposition.

Since there is no need to interpret the weight as the probability
distribution,

the method is free from the sign problem from the outset.].

A common feature of these two approaches is that
one complexifies 
the variables and extends the weight and the observables
as holomorphic functions of the complexified variables. 
In the complex Langevin method,


one generates complexified
configurations using the Langevin equation
with the drift term given by the gradient of

the action.


While the cost of this method is O(N) as in typical MC simulation,

it is applicable only to 
a limited class of systems due to
the conditions for
justification <cit.>. 




On the other hand, the Lefschetz thimble method

is based on the Picard-Lefschetz theory,


in which one deforms the integration contour in the complex space

to a set of steepest descent paths
emanating from some saddle points of the action.
These paths are called the Lefschetz thimbles (or thimbles in short).
The phase of the complex weight is constant along each thimble,
and the sign problem is solved
as far as the phase factor coming from the
integration measure can be taken into account by reweighting <cit.>. 




In practice, the

contour deformation
can be realized
by solving the so-called anti-holomorphic gradient flow equation 
for the variables <cit.>,
where the amount of flow (the flow time)
plays the role of the deformation parameter.

While the Lefschetz thimbles are obtained in the long flow time limit,
one can also use a deformed contour obtained at
finite flow time of the order of log N,
which is expected to be long enough to solve the sign problem
for the system size N.

This is called the generalized thimble method (GTM),
which has a big advantage over the

original Lefschetz thimble method
in that it does not require prior knowledge of the Lefschetz thimbles.
(See Ref. <cit.> for a review.)





When there are more than one thimbles
that contribute to the integral,
the GTM


suffers from

an ergodicity problem
for a large flow time

due to the infinite potential barriers between different thimbles.

For a relatively small system,


one can just choose 
the flow time
long enough to solve the sign problem but not too long
so that the ergodicity problem can be avoided <cit.>.








The range of the flow time that can avoid both problems,
however, shrinks
as the system size increases and eventually vanishes.











In order to solve both problems even in that case,
one can integrate
over the flow time[This proposal is a significant
improvement over
the related ones <cit.>
based on tempering with respect to the flow time,
which requires the calculation of the Jacobian when one swaps the replicas.],
which amounts to treating the flow time as an extra dynamical variable
in the simulation <cit.>.








As an efficient algorithm for numerical simulations in general,
the Hybrid Monte Carlo (HMC) algorithm is widely used.

In this algorithm, one uses a fictitious Hamilton dynamics
to update the configuration.
When one applies this algorithm to the GTM,
there are actually two approaches depending on whether one considers
the Hamilton dynamics on the deformed contour
or on the original contour.
In the first approach,
one has to solve the Hamilton dynamics of
a constrained system in order to make sure that
the configuration remains on the deformed
contour <cit.>.
This requires some complicated calculations since the deformed
contour is given only implicitly by solving the

gradient flow equation.



A big advantage of this approach, however, is that the probability distribution
of the generated configurations

automatically 
includes the modulus of the Jacobian associated with the flow of variables,
which is not the case

in the second approach
as we explain below.



In the second approach,
one can just solve the Hamilton dynamics of
an unconstrained system on the real axis, which makes this part of
the algorithm much simpler than in the first approach.
The force term of the Hamilton equation has to be calculated
by taking the derivative of the action on the deformed contour
with respect to the variables on the real axis, 



which is actually possible without large computational cost
if one uses the idea of backpropagation <cit.>.
A drawback of this second approach, however, is that
the probability distribution of the generated configurations
does not include 
the modulus of the Jacobian associated with the flow of variables,

and therefore it has to be included by reweighting
together with the phase of the Jacobian.





In fact, the modulus may fluctuate considerably
depending on the generated configurations.
In that case, only a few configurations that have a large modulus of the Jacobian
dominate the ensemble, and one cannot increase the statistics efficiently.
This is the so-called overlap problem that
can occur in general when one uses reweighting.
In what follows,
we call these
two approaches the on-thimble approach[This is a bit of
  abuse of the word “thimble” since here it actually represents
  a deformed contour obtained by the
  
  gradient flow
  at finite τ.]
and the on-axis approach, respectively. 


In this paper,
we point out a general problem that occurs when the GTM
is applied to large systems.




In the GTM, the deformation of the integration contour
by the

gradient flow equation is really the key to
solve the sign problem.
The crucial feature of the gradient flow equation
that makes this possible
is that it maps a small region on the real axis
to an exponentially large region on the deformed contour as the flow time
becomes longer.


The growth rate of this exponential behavior, however,
depends on the eigenmode, and it is governed by the singular values of
the Hessian of the action.
As the system size increases, the singular value spectrum typically
exhibits a huge hierarchy.

As a consequence,
if we choose the flow time long enough
to solve the sign problem associated with the eigenmodes
corresponding to small singular values,
the eigenmodes corresponding to large singular values
tend to diverge and easily get out of control[This
  problem
  should not be confused with
  
  the blow-up problem of the flow equation
  discussed in Ref. <cit.>, which occurs
  when the integrand
  
  becomes zero at some point in the complexified configuration space.
  
  In fact,
  
  we consider that the blow-up problem does not occur
  if one uses the preconditioned flow
  that we propose in this paper
  as we discuss 
  in Appendix <ref>.].
  








We show that this hierarchical growth problem of the GTM
can be solved by preconditioning the

gradient flow equation.

For that, we first point out that
we are free to

introduce a Hermitian positive-definite kernel on the
right-hand side of the gradient flow equation without spoiling
its property that is necessary to solve the sign problem.
One can actually use this freedom to make

the growth rate

identical for
every eigenmode,
which implies that
the hierarchical growth problem can be solved by the preconditioned flow.
Moreover, when applied to the on-axis approach,
the preconditioned flow suppresses the fluctuation of the Jacobian,
which can otherwise cause the overlap problem as we mentioned above.


Actual implementation of this preconditioned flow can be done

by using the techniques developed for Rational Hybrid Monte Carlo
algorithm, which is widely used in numerical simulation.
In order to demonstrate
how the preconditioning works,


we apply the GTM

to the real-time quantum evolution of an anharmonic oscillator[See also
  Refs. <cit.>
  for calculations of the real-time quantum evolution using the GTM
  based on the Schwinger-Keldysh formalism, which actually simplifies the
  simulation <cit.>.].

In particular, we show that
the preconditioning enables us to

simulate a large system
that cannot be achieved otherwise.
We also show that
the computational cost
for generating configurations

grows only linearly with the system
size N
for a local system if we implement 

the preconditioning in
the flow appropriately.
Preliminary discussions on the precondition flow
are presented in our previous paper <cit.>
and a proceedings article <cit.>,
where


we establish a new picture of
quantum tunneling in the real-time path integral using the GTM.















The rest of this paper is organized as follows. 
In Section <ref>,
we discuss the hierarchical growth problem of the original flow 


due to the hierarchy in the singular value spectrum of the Hessian. 
In particular, we show that the problem becomes severer when the system
size increases using a simple example of a harmonic oscillator.
In Section <ref>, we explain
our proposal of the preconditioned flow equation,
and discuss how it can be implemented in practice.






In Section <ref>,
we demonstrate how the preconditioning solves the hierarchical growth problem
of the flow as well as the overlap problem in the on-axis approach
by applying the GTM to
the real-time quantum evolution of an anharmonic oscillator.










Section <ref> is devoted to a summary and discussions.

In Appendix <ref>, we
discuss the absence of the so-called

blow-up problem

in
the preconditioned flow.
In Appendix <ref>, we show that
the preconditioned flow actually changes the thimble
and yet it solves the sign problem using a simple example. 
In Appendix <ref>, we provide a brief review
on the application of the HMC to the GTM based on the on-axis approach.
In Appendix <ref>,
we discuss how the idea of backpropagation for calculating the force
term in the HMC algorithm works with the preconditioned flow.
In Appendix <ref>, we present the

parameters
of the GTM
used in our simulations.














§ THE HIERARCHICAL GROWTH PROBLEM IN THE GTM


  



In this section, we discuss the crucial properties

of the

gradient flow equation
that make it possible to solve the sign problem.
In particular, we show
that it maps a region on the real axis
to a region on the deformed contour,
which becomes exponentially large with the flow time.



Moreover the growth rate for each eigenmode is governed
by the singular value spectrum of the Hessian of the action.
When

the system size N is large,
the growth rate tends to exhibit a large hierarchy,
which is

the hierarchical growth problem.





 §.§ Hierarchical growth property of the flow equation



Let us consider a partition function

with a complex action S(x)∈ℂ defined by


  Z=∫_ℝ^Nd^Nx e^-S(x)  .
  
  

In the GTM,
one deforms the integration contour by using
the

gradient flow equation


  d z_j(σ)/d σ
  =∂S(z(σ))/∂z_j  .


By solving this equation
from σ=0 to σ=τ
with the initial condition z_j(0)=x_j,
we obtain a map x ↦ z(x,τ) ≡ z(τ).



Then the deformed contour is defined by
Σ_τ  = {  z(x,τ) | x ∈ℝ^N  },
where τ is the so-called flow time.



Due to Cauchy's theorem, the partition function

can be rewritten as

  Z   =∫_Σ_τd^Nz e^-S(z)   .
  


The crucial

property of the flow equation is


  d S(z(σ))/d σ
  = ∂S(z(σ))/∂z_j
   d z_j(σ)/d σ
  =   |∂S(z(σ))/∂z_j|^2
  > 0  ,

 
which implies that
Re S(z(σ)) increases monotonically
with increasing σ,
while Im S(z(σ)) remains constant.
In the infinite flow time limit τ=∞,

the deformed contour Σ_∞
is composed of

a set of Lefschetz thimbles.
Each thimble can be obtained by the flow from some saddle
point (defined by S(z)/ z_j = 0),
and hence Re S(z(σ)) increases monotonically
keeping Im S(z) constant as one flows away from the
saddle point.
This implies that the sign problem is solved on each thimble
as far as

the phase factor coming from the integration measure d^N z
in (<ref>)
can be treated by reweighting <cit.>. 





The key mechanism for solving the sign problem here
is that there is some point P on ℝ^N which flows
into a saddle point in the τ→∞ limit,
and
the thimble associated with that saddle point
is obtained by mapping an
infinitesimal vicinity of the point P using the flow.


This property of the flow, together with the fact that
Im S(z) is kept constant along the flow,
makes it possible to solve the sign problem.





In fact, the sign problem can be solve even at finite τ
as far as τ

is large enough to suppress the fluctuation of Im S(z).
















Here we discuss in more detail
how a small region on the original contour ℝ^N
is extended

in the deformed contour.
For that, let us consider two points x and x+δ x
on the real axes, which differ only by some infinitesimal displacement δ x.

The displacement
ζ(σ) ≡ z(x+δ x,σ)-z(x,σ)
of the two points after some flow time σ

is given by the equation





  d ζ_j(σ)/d σ
  = H_jk(z(σ))   ζ_k(σ) 
  


with ζ_j(0)=δ x_j, 
where H is the Hessian of the action defined by

    H_ij(z)
         = ∂^2 S(z)/∂ z_i∂ z_j .



In order to solve (<ref>),
let us rewrite it as

    d/d σ[   ζ(σ); ζ̅ (σ) ]   =
     H(σ)
    [  ζ(σ); ζ̅(σ) ] ,

where
we use the vector notation
ζ = (ζ_1 , ⋯ , ζ_N)^⊤,
ζ̅ = (ζ̅_1 , ⋯ , ζ̅_N)^⊤
and define the 2N × 2N Hermitian[Note that the Hessian
  is symmetric H^⊤ = H and hence H̅=H^†.]
matrix H(σ) as

    H(σ)    = 
      [         H(z(σ)); H(z(σ))         ] .


Thus the solution to
the differential equation (<ref>)

can be written down formally as

    [   ζ(τ); ζ̅ (τ) ]   =
     Pexp( ∫_0 ^τ  d σ  H(σ) )
    [ δ x; δ x ] ,

where
δ x = (δ x_1 , ⋯ , δ x_N)^⊤


and Pexp represents the path-ordered exponential,
which ensures that H(σ) with smaller σ
comes on the right after Taylor expansion.


In order to understand the behavior of the solution (<ref>),
let us diagonalize the Hermitian matrix H(σ).
For that, we





consider the singular value decomposition (SVD)
of the Hessian H_ij(z(σ)) given as[This is known
  as the Takagi decomposition,
  which is the SVD for a complex symmetric matrix.]

    H(z(σ)) = U^⊤ (σ)  Λ(σ)   U(σ)  ,

where U is a unitary matrix and
Λ= diag(λ_1 , ⋯ , λ_N)
is a real diagonal matrix with positive[Here we assume that
  there are no zero singular values for simplicity. However, our idea can
  be easily extended to the case in which this assumption does not hold.
See footnote <ref>.
]
entries λ_j > 0.

Using Λ(σ) and U(σ), we can diagonalize H(σ) as

    H(σ)
         =  U^†(σ)
      [   Λ(σ)       ;        - Λ(σ) ] U(σ)  , 
      
    
    
    
    
       U(σ)
         =
      1/√(2)[    1_N    1_N; -  1_N    1_N ][ U(σ)     ;      U(σ) ] .

This implies that

the displacement grows exponentially with σ
for a small region of σ, in which H_jk(z(σ)) can be regarded
as constant.
Moreover,




the growth rate

is governed by the 
singular values of the Hessian H(z), and it depends on the
eigenmode ζ_ E^(i) =  Re(U_ijζ_j) ∼ e^λ_i σ.
When the singular value spectrum has a large hierarchy,


the growth rate becomes very different for different eigenmodes.





Since the flow time τ is common to all the eigenmodes,
if we choose τ in such a way that the sign problem is solved
for the eigenmodes ζ_ E^(i)
with small λ_i,


the eigenmodes ζ_ E^(i)
with large λ_i tend to diverge.

 
This problem is severer in the on-axis approach
since one has to update the configuration on the real axis.
In order to avoid the hierarchical growth problem, one has to choose a very small
step size in solving the Hamilton dynamics of the
HMC

algorithm,
which makes the simulation extremely slow.
In the case of the on-thimble approach, one updates the
configuration on the deformed contour directly, but the problem

may occur when one solves the flow equation to
ensure that the updated configuration is still on the deformed contour.
(See Section <ref> for further discussions.)








 §.§ Example of a harmonic oscillator

 

In this section, we show that the hierarchical growth problem discussed in the previous
section indeed occurs in a simple example of a harmonic oscillator.
This example provides a clear understanding of the problem since the
gradient flow equation can be solved explicitly.



The path integral that describes

the real-time quantum evolution of a harmonic oscillator
is given by

(<ref>)
with the action





  S   =-i∑_j=0^Nϵ{  1/2
  (x_j+1-x_j/ϵ)^2
 - 1/2   m^2   x_j+1^2+x_j^2/2}  ,



where we fix x_0 = x_N+1=0 for simplicity and
treat x_1,⋯,x_N as variables.

The path integral (<ref>) then
represents
the transition amplitude from the origin to the origin
in time T=ϵ (N+1) up to some known normalization factor. 

In this example,
the Hessian H(z) of the action defined by
(<ref>) is a constant matrix,
which does not depend on the complexified configuration z_j.

Furthermore,

it is not only symmetric but also pure imaginary.
Therefore, we can diagonalize it as


    H = -i    O^⊤ Λ̃  O  ,

where O∈ SO(N)

and Λ̃ is a real diagonal matrix
Λ̃ = diag(λ̃_1, ⋯ , λ̃_N)
with


  λ̃_j =
  ϵ{ 4/ϵ^2sin^2πj/2(N+1)-m^2}  .
  

Note that λ̃_j can be negative in general,
and the singular values[The singular value decomposition (<ref>)
  of H= U^⊤ Λ  U is given by U=V   O,
  where V is a diagonal unitary matrix with the entries given by
  either e^π/4i or e^-π/4i depending on
  whether the corresponding eigenvalue λ̃_j is
  negative or positive, respectively.]
are given by λ_j = |λ̃_j|.













Using the Hessian H and its diagonalized form (<ref>),
the action can be written as

    S= 1/2x_jH_jkx_k
       = -  i/2∑_j=1^N λ̃_j (y_j)^2   ,

where we have defined the variables y_j=O_jk x_k.
Thus the integral (<ref>)

is rewritten as


  Z   =
  
  ∫_ℝ^Nd^Ny exp(i/2 ∑_j=1^N λ̃_j   y_j^2)  .




The flow equation for the complexified variables z_j(y,τ) is given by

  ∂z_j(y,σ)/∂σ
     =i λ̃_jz_j(y,σ)  ,


z_j(y,0)   =y_j  ,

whose solution

can be readily obtained as

  z_j(y,σ)={cosh(λ̃_jσ)
  + isinh(λ̃_jσ)}   y_j  .
  

Note that the gradient flow magnifies the displacement δ y_j on the real axis
by a factor of

    R_j(σ) ≡√(cosh^2(λ̃_jσ) + sinh^2(λ̃_jσ)  )
      = √(cosh(2λ̃_jσ) )∼
      e^λ_j σ ,
 
which grows exponentially with σ, and the growth rate is given by
the singular value λ_j = |λ̃_j|
as suggested by the general discussion in the previous section.


Using (<ref>),

the integral on the deformed contour is written as 

  Z   =∫_Σ_τd^Nz exp(i/2
  ∑_j=1^N λ̃_j   z_j^2)


   =  ∫_ℝ^Nd^Ny  det   J    
exp(- 1/2
∑_j=1^N λ̃_j {sinh(2λ̃_jτ)
- i}   y_j^2)  .


Here
the Jacobi matrix J associated with the deformation of the contour is given by

  J_jk   =∂z_j(y,τ)/∂y_k
  
  = { cosh(λ̃_jτ)+ isinh(λ̃_jτ) }   δ_jk  ,
  

which is independent of y and hence can be factored out of the integral.


From (<ref>),
one can see how the sign problem is solved by increasing τ. 
The real part of the argument of the exponential function implies
that the region of y_j that contribute to the integral is |y_j| ≲Δ_j,
where

    Δ_j =
      
      {λ̃_j sinh(2λ̃_jτ) }^-1/2 .


In that region, the phase factor of the integrand is close to unity as far as

    1/2  |λ̃_j| (Δ_j)^2 ≪π .
 
Combining (<ref>) and (<ref>),
the sign problem associated with the eigenmode j is solved if

    τ≫1/2 λ_j arcsinh( 1/2π)  .
 
Thus, in order to solve the sign problem of the whole system,
one has to satisfy (<ref>)
for the smallest singular value λ_ min.
In that case, the magnification factor (<ref>)
for the largest singular value λ_ max becomes

    R_ max(τ)    ∼exp( c  λ_ max/λ_ min)  ,
      
           c    ≫1/2 arcsinh( 1/2π) ∼ 0.08  .

When the singular value has a large hierarchy 
λ_ max/λ_ min≫ 1,
R_ max(τ) can easily diverge,
which is nothing but the hierarchical growth problem.

The ratio η(H) ≡λ_ max/λ_ min
represents

the condition number of the matrix H.
Let us evaluate this quantity

explicitly in the case of a harmonic oscillator
using (<ref>).
For mT < π, all the eigenvalues become positive λ̃_j>0
for sufficiently large N. In that case, η(H) can be easily evaluated as



η(H)   =4(N+1)^2sin^2Nπ/2(N+1)-(mT)^2/4(N+1)^2sin^2π/2(N+1)-(mT)^2



    ∼
4/π^2-(mT)^2   N^2  ,

in the N→∞ limit with fixed m^2 and T,
which corresponds to the continuum limit.

In Fig. <ref>,
we plot η(H) against N for m^2=1 and T=2.
By plugging these values in (<ref>), one finds that
R_ max is as large as 2.4× 10^9

even for N = 20

used in our numerical simulation in
Section <ref>.


























 








Let us also note that the 
Jacobian that appears in (<ref>)

has a modulus that can be written
in terms of the magnification factor
R_j(σ) defined in (<ref>) as

  | det   J |
  
     =∏_j=1^N R_j(τ)
  

  
  ∼e^Λτ  ,



where Λ = ∑_j=1^N λ_j grows quadratically
with N in the continuum limit.
In the present example of a harmonic oscillator,
the Jacobian is a constant that does not depend on the variables y_j,

which simply factors out of the integral in (<ref>)
as we already mentioned.
This is not the case

in a general system, however.


The Jacobian with
a huge modulus depending on the configuration
can cause an overlap problem in the on-axis approach as we discussed
in the Introduction.





In the example discussed above,

the origin of the large
η(H) is the ratio of the momenta in the UV and IR
regions as one can see from (<ref>).



This is closely related to the recent discussion


that the flow equation in the continuum theory
does not allow a well-behaved solution due to the high frequency
modes <cit.>. 








In order to solve this problem,

it was proposed to
modify the real part of the action that appears
on the right-hand side of the flow equation

in the vicinity of the saddle points.





As we have seen in Section <ref>,
a similar problem occurs
generally in a discretized theory with
a large number

N of variables,




which

typically exhibits
a large hierarchy in the singular value spectrum of the Hessian.





We solve this problem

by

modifying the flow equation so that the growth rate of each mode becomes
identical,
while keeping the crucial property of the flow intact.























§ A SOLUTION TO THE HIERARCHICAL GROWTH PROBLEM




Our basic strategy to solve the hierarchical growth problem
is to normalize the growth rates at each step of the flow.
For instance, in the example of a harmonic oscillator discussed in the
previous section,
we can replace (<ref>) by

  ∂z_j(y,σ)/∂σ
     =i   λ̃_j/|λ̃_j|  
  z_j(y,σ)  ,


without spoiling the crucial property (<ref>)
of the flow equation.
Thus the hierarchical growth problem can be completely solved.
What we aim to do here is to generalize this prescription to an
arbitrary system.



 §.§ Preconditioned gradient flow equation

  
Let us first note that the crucial property (<ref>)
of the gradient flow equation holds for
more general flows described by[Recently, it has been found that
  introducing a kernel in the complex Langevin equation is
  useful in stabilizing simulations in the CLM for real-time quantum
  evolution <cit.>.]




d z_j/d σ
   = 𝒜_jk(z,z̅)   ∂S/∂z_k  ,


with 𝒜(z,z̅) being an arbitrary positive-definite hermitian matrix,
which is not necessarily a holomorphic function of z.
This can be proved easily as


d S(z(σ))/d σ
=∂S/∂z_j𝒜_jk(z,z̅)
∂S/∂z_k>0  .

Note that 𝒜 changes not only the deformed contour Σ_τ
but also the thimbles Σ_∞.

However, since the crucial property of the flow is kept intact,
Im S is constant on each thimble, which implies that
the sign problem can be solved in the same way.
We discuss this point
in Appendix <ref> with an explicit example.









For the generalized flow (<ref>),
we find that the displacement
ζ(σ) ≡ z(x+δ x,σ)-z(x,σ)
of the two points after some flow time σ

is given by the equation

    d/d σζ_i(σ)
    
      =𝒜_ikH_kl(z(σ)) ζ_l(σ)
      + (
      ∂𝒜_il/∂ z_kζ_k(σ)
      + ∂𝒜_il/∂z̅_kζ_k(σ))
        ∂ S(z(σ))/∂ z_l ,





instead of

the original one (<ref>).
Let us here assume that the first term is
dominant[This assumption is valid when z(σ) is
  close to a saddle point, for instance. Otherwise, it should be 
simply regarded as a working hypothesis.]
in (<ref>).
Then

the solution to (<ref>)
can be written formally as (<ref>)
with H(σ) replaced by[Note that
  H(σ) is not Hermitian in general
  unlike H(σ).]

    H(σ)    = 
      [               A H(z(σ)); A̅  H(z(σ))             ] .



Recalling

the singular decomposition (<ref>) of H(z(σ)),









we can choose

    𝒜 = U(σ)^† Λ^-1(σ)  U(σ)  ,

so that
the problematic hierarchy of singular values 

in H(z(σ)) is completely eliminated as

    A̅  H(z(σ))
       = U^⊤ (σ)   U(σ)  .













In this case, the 2N × 2N matrix H(σ)
is again Hermitian and it has N eigenvalues of 1 and -1, respectively;
namely the growth rate of each eigenmode becomes identical.
(See (<ref>).)
From this point of view, (<ref>) can be regarded
as the optimal choice for the “preconditioner” 𝒜 
in the generalized

flow equation (<ref>).



In order to implement this idea in practice,

let us note 
that (<ref>) can be written
as[
  When H has zero singular values, the preconditioner 𝒜
  has to be regularized, for instance,
  as 𝒜=(H̅(z̅)H(z)+ε)^-1/2
  with a small positive ε.
  This is needed
  
  when one treats a system with symmetries. 
  While the deformed contour changes with ε, the integral
  remains unaltered.
  
] 

    𝒜(z(σ),z(σ)) 
      
         =
      { H^†(z(σ)) H(z(σ)) }^-1/2
      = {H(z(σ)) H(z(σ)) }^-1/2 .

Here we use the rational
approximation










    x^-1/2   ≈
    a_0  + ∑_q=1^Q a_q/ x + b_q  ,

which can be made accurate for a wide range of x>0

with the real positive parameters a_q and b_q
generated by the Remez algorithm.




Thus we obtain

    𝒜(z,z̅)
           ≈ a_0  1_N
        + ∑_q=1^Q a_q  {H(z) H(z)+b_q 1_N }^-1 .

The matrix inverse (H̅H+b_q1_N)^-1 does not have to
be calculated explicitly since it only appears in the algorithm 
as a matrix that acts on some

vector, which allows us to
use an iterative method for solving a linear equation
such as the conjugate gradient (CG) method.
The factor of Q in the computational cost
can be avoided by the use of a multi-mass CG
solver <cit.>.
These techniques are well known 
in the so-called Rational HMC
  algorithm <cit.>,

which is widely used in

  QCD with dynamical strange quarks <cit.> and
  supersymmetric theories such as the BFSS and IKKT matrix 
models (See Refs. <cit.>, 
for example.). 




 §.§ Calculation of the Jacobian



In the GTM, one has to calculate
the Jacobian associated with the change of variables
defined by the gradient flow.
In the on-thimble approach, only the phase factor of the Jacobian
has to be reweighted,

while in the on-axis approach,
not only the phase factor but also the modulus has to be reweighted.

The calculation of the Jacobian
gets modified when one introduces the preconditioner
to the flow equation.
In this subsection, we discuss how this can be done efficiently.


Let us note first that
the flow of the Jacobi matrix is given
by[The flow equation (<ref>)
  for the displacement
  
  can be obtained
from (<ref>) by using
ζ_j(σ) = J_jk(σ)  δ x_k.
] 

    d/d σJ_ij(σ)
    
      =𝒜_ikH_kl(z(σ))  J_lj(σ)
      + (
      ∂𝒜_il/∂ z_k J_kj(σ)
      + ∂𝒜_il/∂z̅_kJ_kj(σ))
        ∂ S(z(σ))/∂ z_l ,

for the preconditioned flow equation (<ref>).
(The corresponding flow for the original flow equation (<ref>)
can be retrieved by setting 𝒜 to an identity matrix.)
Using the expression (<ref>),

the derivative of 𝒜
in (<ref>)

can be calculated as


    ∂𝒜/∂ z_k   =-∑_q=1^Qa_q (H̅H+b_q1_N)^-1H̅∂ H/∂ z_k
      (H̅H+b_q1_N)^-1 ,
    
      
    ∂𝒜/∂z̅_k   =-∑_q=1^Qa_q(H̅H+b_q1_N)^-1∂ H/∂ z_k
      H(H̅H+b_q1_N)^-1 .


Note that
the matrix inverse (H̅H+b_q1_N)^-1 appears twice
in the above expressions.
In order to use

the idea of the multi-mass solver,
which enables us to avoid repeating the CG procedure Q times
(See the explanation below
Eq. (<ref>).),

we have to apply the matrix inverse to some different vectors separately.
For that, we regard the matrix equation
(<ref>)
as a set of vector equations for i=1 , ⋯ , N.
We can calculate the right-hand side of the vector equation for
a particular i=I by multiplying a unit vector e_i^(I)=δ_iI
to the right-hand side of (<ref>).
Thus
the matrix inverse
(H̅H+b_q 1)^-1 that appears
in
(<ref>)


can be applied to

the two vectors e_i^(I) and ∂ _l S
separately.










In fact, the expressions (<ref>) have to be used


also when we compute the force term in the HMC algorithm
in the on-axis approach. (See Appendix <ref>
for a brief review.)
In that case, too, we can apply the idea of the multi-mass solver

thanks to the structure of calculations in the backpropagation
as we discuss in Appendix <ref>.


The calculation of the Jacobi matrix that we discussed above

requires O(N^2) computational cost for a local action
since it involves multiplications of a sparse matrix 
to a dense matrix[For a non-local action,

  the cost increases
  by a factor of O(N) since the sparseness is lost.].

On the other hand,
the calculation of its determinant
requires O(N^3) computational cost,


which is not affected by the preconditioning at all.
These calculations appear only in the reweighting procedure,
which can be done off-line after generating configurations,
and the procedure can be
  parallelized trivially without any overhead due to
  communications.






On the other hand,
all the calculations

for generating configurations

require
only O(N) computational cost for a local action
even after implementing the preconditioner in the gradient flow equation
since it just involves multiplications of a sparse matrix to a
vector.










































In passing,
let us also mention that

there is actually
a cheaper preconditioner,
which can be obtained by considering only the free part of the action
when we calculate the Hessian H to be used in (<ref>).
Since the preconditioner 𝒜(z,z̅) in this case becomes
a constant matrix independent of the configuration,
we only have to compute 𝒜
and diagonalize it once and for all.
The rational approximation (<ref>)
is not needed any more.
This cheaper preconditioner works efficiently in the weak coupling regime,
and the computational cost can be reduced typically by

an order of magnitude compared to the full preconditioner.
























§ DEMONSTRATION OF THE PRECONDITIONED FLOW





In this section, we demonstrate
how the preconditioned flow equation works
in the application of GTM to the real-time quantum evolution
of an anharmonic oscillator.



Here we use the on-axis approach for the HMC algorithm,
where the force term is calculated by
the backpropagation <cit.>.

In particular,
we show that
the long autocorrelation in the generated
configurations

is drastically reduced by preconditioning the flow.
We also show that the preconditioning
drastically reduces the modulus of the Jacobian
and its fluctuations,


thereby solving the overlap problem in the on-axis approach.









 §.§ Simulation setup



The system we deal with in this section

is the same as
(<ref>) except that
we consider an anharmonic potential

    V(x)    =
    1/4! λ  x^4  .



Here we introduce the initial wave function

    ψ_ i (x) ∝exp{-1/4γ 
      (x-x_ i)^2} ,

and calculate the time-evolved wave function ψ_ f(x)
after some time T,
which is given by the path integral (<ref>)
up to some known normalization factor
using the action

S(x;x_f,x_i)=
   -i∑_j=1^Nϵ{1/2 (
x_j+1-x_j/ϵ )^2
-  V(x_j+1) + V(x_j)/2}


+1/4   γ  (x_1-x_i)^2 ,


where T=Nϵ and x_N+1≡ x_ f.



Similar calculations have been

done

using the original flow
in Ref. <cit.>
although N had to be restricted to small values such as N=9.


The main point here is that the preconditioned flow enables us
to increase N without
any problems.















In order to avoid the ergodicity problem

concerning

multiple thimbles, we integrate over the flow time <cit.>
as described in Ref. <cit.>

in the on-axis approach.

The Hamiltonian to be used in the fictitious Hamilton dynamics
of the HMC algorithm is

given by



  H=1/2m^2(τ)∑_j=1^Np_j^2
   + 1/2m̃^2   p_τ^2
  + Re  S(z(x,τ)) +  W(τ)  ,


where p_j and p_τ are the conjugate momenta corresponding
to x_j and the flow time τ, respectively.
The function m(τ) 
and the parameter m̃
are introduced for optimization as we explain shortly.
The total time s_ f and the step size Δ s
for the fictitious Hamilton dynamics are the parameters of the HMC algorithm
that can be optimized in a standard manner.
In all our simulations,
we choose s_ f=1 and Δ s=0.05,
and the number of time steps for solving the flow equation is set to
N_τ=10.

The parameter m̃>1 in (<ref>) is introduced
only for simulations with the original flow

in order to avoid the force in the τ-direction getting too large depending on x,
which causes the drop of the acceptance rate.
We have increased m̃ from unity until the acceptance rate becomes
reasonably high.




On the other hand, the function m(τ) 
in (<ref>)
is chosen to be
the typical largest singular value of the Jacobi matrix J(x,τ)
based on

the discussion in Section <ref>





so that the eigenmode

with the largest growth rate does not diverge.[
  In Ref. <cit.>,
  
  m(τ) was chosen to be

  the typical value of |det J(x,τ)|^1/N,




which
corresponds to taking the geometric average of the growth rate.



However, the effective stepsize in this case 
becomes too large for the eigenmode with the largest growth rate,
which causes divergence during simulations at large N. 



] 



Then the hierarchical growth problem manifests itself as

long autocorrelation due to the modes with small growth rate.

The function W(τ)
in (<ref>) is determined so that the distribution
of τ becomes flat within the region τ∈ [τ_ min,τ_ max],
where τ_ max has to be chosen to be
large enough to solve the sign problem
and τ_ min has to be chosen to be
small enough to solve the ergodicity problem.
The actual form of the function W(τ) as well as that of m(τ)
is determined iteratively by improving them step by step
using the results obtained by the previous simulation.





















Let us discuss how we





determine
the optimal region [τ_min, τ_max]
of τ to be used in the simulation,
which actually depends on whether we use the preconditioned flow or not.


















 















In Fig. <ref>,
we show the distribution
of the phase of the integrand in (<ref>)
including the phase coming from the Jacobian
for the original flow (Left) and the preconditioned flow (Right)




in the case of
N=6, T=2, λ=1, x_ i=1, γ=1 and x_ f=0.













The sign problem is solved in
the region of τ that exhibits some non-uniformity
in the distribution of the phase.

However, in order to avoid the ergodicity problem,
one should also sample the region of τ
in which the distribution of the phase is almost uniform.
Based on these criteria,
we choose the range of τ to be

[0.02, 0.2] for the
original flow,
and 

[0.4, 0.8] for
the preconditioned flow in this case.


See Appendix <ref> for

the choice of the parameters

m̃, τ_ min, τ_ max, m(τ) and W(τ)
determined in the way described above for each case.






 §.§ Reduction of the autocorrelation
 



Let us first show that the autocorrelation is drastically
reduced by preconditioning the flow equation.

We define the autocorrelation
in the generated configurations x = (x_1 , ⋯ , x_N )


as follows.
Let x^(k) (k = 1, ⋯ , n)
be the k-th configuration.
We denote the average of the n configurations as
x = 1/n∑_k=1^n  x^(k).
Then the autocorrelation

can be defined as

  
  C(k)
  ≡1/v 1/n-k∑_m=1^n-k
( x^(m)-  x )
·(x^(m+k)- x )  ,





where v is the variance defined by

  v =   x ·x 
  -   x ·x     .




  












In Fig. <ref>,
we show the autocorrelation
of x=(x_j) obtained by simulations
with
the same parameters as in Fig. <ref>.
The region of τ is chosen to be
[0.02, 0.2] and [0.4, 0.8]
for the original flow and the preconditioned flow, respectively.





We find that the autocorrelation is reduced by a factor of 10
for the preconditioned flow, which suggests
that the hierarchical growth problem is avoided.




























 §.§ Reduction of the modulus of the Jacobian





Next we discuss the Jacobian associated with the change of variables
by the gradient flow, which can typically have a large modulus
as we mentioned below (<ref>) in Section <ref>.

This, in particular, causes the overlap problem in the on-axis approach.
We will see below that the preconditioned flow reduces
the modulus of the Jacobian and its fluctuation drastically,
thereby solving this problem as well as the hierarchical growth problem.

In

Fig. <ref> (Left),
we show the scattered plot of log_10|det J|

obtained by simulations
for various N (4 ≤ N ≤ 20)

with T=2, λ=1, x_ i=1, γ=1 and x_ f=0.


We have fixed the range of τ as





[0.02, 0.2] and [0.4, 0.8] 
for the cases with and without preconditioning, respectively. 
The blue crosses and the orange circles
represent the results for the original flow
and the preconditioned flow, respectively.



















Our results show that 
the preconditioning reduces the values of
|det J| drastically as expected.




Furthermore, from Fig. <ref> (Right),
we find that the preconditioning suppresses
the fluctuation of |det J| drastically,
which suggests that it can

also solve the overlap problem in the on-axis approach.



Without preconditioning,
it was not even possible to

perform simulations for N>16

since the typical magnitude of the HMC force in the τ-direction
changes too much with x, and we were not able to control
it by the x-independent potential

W(τ).
This may be viewed as

another manifestation of the hierarchical growth problem.



Such behaviors were
not seen in simulations with

the preconditioned flow.









 §.§ Result for the time-evolved wave function





Finally we show that the preconditioned flow enables us
to obtain

results for
the time-evolved wave function ψ_ f(x_ f)
even at large N and at strong coupling (large λ).


As has been done in Ref. <cit.>,
we calculate









  ∂/∂x_f   log(ψ_f(x_f))
  

  
  =- ⟨∂S/∂x_f ⟩
  
  =⟨i(x_f-x_N/ϵ
   - ϵ/2   V'(x_f ))⟩ ,


which is directly accessible by calculating
the expectation value on the right-hand side.
In Fig. <ref>,
we show our results for

N=20, T=2, λ=30, x_ i=0.3 and γ=4.
For each x_ f, we have determined the range of τ
as presented in Appendix <ref>.
Our results are in good agreement with the results obtained by
solving the Schrödinger equation with Hamiltonian diagonalization.

















In Ref. <cit.>, it was difficult to obtain reliable results
as one goes

beyond N=9
because of the overlap problem
even at λ=1.









Here we are able to obtain results for N=20












even at strong coupling λ=30 

without any problems.



















§ SUMMARY AND DISCUSSIONS




In this paper we have pointed out a problem in the gradient
flow equation, which is used
in the GTM to deform the integration contour into
the complex plane. The property of the flow
that plays an important role in solving the sign problem
is that it maps a small region on the
real axis to a region on the deformed contour
which becomes exponentially large with the flow time.
The problem is that the growth rate for each mode
typically has a huge hierarchy


when the system size becomes large.
If one chooses the flow time to be large enough to solve the sign problem
associated with the slowly growing modes,
the fast growing modes tend to diverge.

In order to solve this hierarchical growth problem,
we have proposed to modify the flow equation

by the preconditioner, which makes the growth rate equal
without spoiling the crucial properties of the flow equation.
This preconditioner can be implemented practically in the GTM
with the standard techniques used in the Rational HMC algorithm.

We applied this method to the 

real-time quantum evolution

of an anharmonic oscillator
using the on-axis approach with backpropagation for calculating the
HMC force <cit.>.


In the on-axis approach, the hierarchical growth problem manifests itself in the long
auto-correlation time since one cannot choose the parameters
in the HMC algorithm for each mode separately.
Our results indeed show that the auto-correlation time is reduced drastically.
Moreover, the modulus of the Jacobian, which has to be taken into account
by reweighting, is also reduced drastically.
This solves the overlap problem in reweighting, which is caused
by large fluctuations of the modulus of the Jacobian.

Thus we were able to simulate the system with the size N
that was not accessible

without preconditioning the flow <cit.>.
Note also that
this is achieved with a very strong coupling λ=30 in the
quartic potential (<ref>), which makes the
flow equation highly non-linear.


As a side remark,

the preconditioned flow

solves

yet another problem of
the original gradient flow, which occurs

when
the integrand has zeroes <cit.>.
In that case, there are points on the real axis that flow
into the zeroes within finite flow time,
and the right-hand side of the flow equation

blows up
as one approaches the zeroes.

A practical solution proposed in Ref. <cit.>
is to


add some regulator to the flow equation,
which modifies the flow only in the vicinity of the zeroes.
We consider that

this problem does not occur for the preconditioned flow

since the preconditioner cancels the divergence of the right-hand side
that occurs as one approaches the zeroes.
See Appendix <ref> for 

discussions on this issue.








While we focused on the on-axis approach in this paper,
the hierarchical growth problem of the flow
may also 

affect the on-thimble approach.




There, we have to solve a set of


equations iteratively by Newton's method
in order to make sure that the updated configurations
are constrained

on the deformed contour. 

This procedure appears at every step of the fictitious time evolution
in the HMC algorithm.

The coefficient matrix of these equations
involves the Jacobian, which has a huge condition number for the original flow.


Therefore, the solution of the equations may
suffer from large numerical errors and/or slow convergence
when the system size becomes large.

The preconditioned flow may be useful in solving such problems.
(See Ref. <cit.> for a new proposal related to this issue.)























As we have shown in this paper,
the preconditioned flow enables us to apply
the GTM to much larger systems than those accessible with the original flow.

In particular, it gives us an access to the continuum limit

and to the strong coupling regime,


which enabled us to establish a new picture of quantum tunneling
in the real-time path integral formalism <cit.>.



We expect that there are many other applications awaiting us to explore.




 §.§ Acknowledgements

We would like to thank Yuhma Asano, Genki Fujisawa,
Masafumi Fukuma and Nobuyuki Matsumoto


for valuable discussions.
The computations were carried out on


the PC clusters in KEK Computing Research Center
and KEK Theory Center.
K.S. was supported by the Grant-in-Aid for JSPS Research Fellow,
No. 20J00079.
A.Y. is supported by a Grant-in-Aid for Transformative Research Areas
“The Natural Laws of Extreme Universe—A New Paradigm for Spacetime and
Matter from Quantum Information” (KAKENHI Grant No. JP21H05191) from
JSPS of Japan.










§ THE ABSENCE OF THE BLOW-UP PROBLEM




Here we discuss
a problem that occurs when the integrand of the partition function has zeroes.


In such a case, the gradient flow can reach the zeroes within finite flow time
and the right-hand side of the flow equation blows up.
This is

the blow-up problem of the GTM, which was discussed in Ref. <cit.>.
Let us first emphasize that
this problem occurs
although the partition function is totally well defined.

The deformed contour one obtains

at sufficiently long flow time
is nothing but a set of contours connected at the

zeroes[See Ref. <cit.>
for recent discussions on this case.] .
It is therefore just a technical problem of how
to deal with the divergence that occurs while solving the flow equation.













Here we argue that this problem

is naturally avoided in the preconditioned flow
since
it

slows down the flow as one approaches the zeroes.







To simplify the argument, we consider a single-variable case N=1
in (<ref>) with real S(x).
In this case, the original flow equation (<ref>) reads


  d x(σ)/d σ

     = S'(x(σ))   ,


whereas the preconditioned flow (<ref>) with
(<ref>) reads


  d x(σ)/d σ
     =
S'(x(σ))/|S”(x(σ))|  ,
  


where ' represents the derivative with respect to x.



For instance, let us consider the case with

    S(x)    = 1/2(n+1)  x^2(n+1) ,      n=1,2,⋯ ,


which gives e^-S(x)→ 0 for x →±∞.

The saddle point is x=0, and there are flows towards the
singularities at x = ±∞.
Solving the original flow equation (<ref>), one obtains

    x(σ)    =
    {[   { x(0)^-2n - 2n σ}^-1/2n                x(0) > 0  ,; - { x(0)^-2n - 2n σ}^-1/2n                x(0) < 0  , ].
 
which reaches the singularities within finite flow time
σ= 1/2nx(0)^-2n.
However, for the preconditioned flow (<ref>),
one obtains

    x(σ)    = x(0)   e^σ/(2n+1) ,
 
which reaches the singularities only in the σ→∞ limit.




Next we consider the case with[Strictly speaking, the partition function
  is not finite
  
  in this case.
  In order to make it finite, one can add a term like

1/2ϵ x^2 in the action, which clearly does not affect the
flow near the singularity x=0.]

    S(x)    = - log (x^2n)   ,      n=1,2,⋯ ,


which gives e^-S(x)→ 0 for x → 0.

The saddle point is x=±∞, and there are flows towards the
singularity at x = 0.
Solving the original flow equation (<ref>), one obtains

    x(σ)    =
    {[    √(  x(0)^2 - 4n  σ)            x(0) > 0  ,; -  √(  x(0)^2 - 4n  σ)            x(0) < 0  , ].
 
which reaches the singularity within finite flow time σ= 1/4nx(0)^2.
However, for the preconditioned flow (<ref>),
one obtains

    x(σ)    = x(0)   e^-σ ,
 
which reaches the singularity only in the σ→∞ limit.


Thus the preconditioning (<ref>)
makes the singularities unreachable within finite flow time.
While we have discussed a single variable case

with a real action
for simplicity,
this property of the preconditioned flow is considered to be quite general.

For instance, the quantum mechanical system with
the anharmonic potential (<ref>)
suffers from a blow-up problem for sufficiently long flow time
if one uses the original flow equation
similarly to the example (<ref>).
We consider that the preconditioned flow naturally avoids this problem.

A more careful study on this issue shall be left for future investigations.





§ MODIFICATION OF THE THIMBLES BY PRECONDITIONING




The preconditioning of the gradient flow equation changes the
deformed contour.
In fact, it also modifies the Lefschetz thimbles,
which appear as the deformed contour in the long flow time
τ→∞ limit.
On the other hand, the preconditioning does not alter
the saddle points, which are defined by S(z)/ z=0
independently of the flow.

In this section, we discuss these points
using a simple example. 










Let us consider

an integral (<ref>)
with N=2 variables, where the action is given by 

S(z)=1/2(z_1^2-iz_2^2)  . 

The saddle point is given by z_1 = z_2 = 0.
The original flow equations read

  d z_1/d σ
  =z̅_1  ,     ∂z_2/∂σ=iz̅_2  ,
  

which is linear.

Using

the real-variable notation 
z_1=x_1+iy_1, z_2=x_2+iy_2,
the flow equation (<ref>)
can be rewritten as



  d/d σ
  ([ x_1; y_1; x_2; y_2 ])
    =
    M
    
    
    
    
    
    
([ x_1; y_1; x_2; y_2 ])  ,     M =
      ([  1         ;    -1      ;           1;        1    ])  .

The Lefschetz thimble Σ_∞ is spanned by
the eigenvectors of the matrix M corresponding to positive
(degenerate) eigenvalue 1, which are given by







v_1=([ 1; 0; 0; 0 ])  ,     v_2=([ 0; 0; 1; 1 ])  .


Next we consider the generalized flow (<ref>)
with

a preconditioner[Note that
  
  the optimal choice is actually 𝒜=  1 in the present case.
  The aim here is to demonstrate that nontrivial 𝒜 modifies
  the thimble.]



  𝒜
  =([ 5 1; 1 5 ])  .



The preconditioned flow equations are given by

  ∂z_1/∂σ=5z̅_1+iz̅_2  ,
      ∂z_2/∂σ=z̅_1+5iz̅_2  .
  

Using the real-variable notation,
the preconditioned flow equation (<ref>)

becomes

  d/d σ
  ([ x_1; y_1; x_2; y_2 ])
    =
    M
    
    
    
    
    
    
([ x_1; y_1; x_2; y_2 ])  ,     M =
      ([  5        1;    -5  1   ;  1        5;    -1  5    ])  .

The Lefschetz thimble Σ_∞ is spanned by
the eigenvectors of the matrix M
corresponding to positive eigenvalues 4√(2) and 3√(2),
which are given, respectively, by



ṽ_1=([ 5 + 4√(2); 5 - 3√(2); 1 + 5√(2);         7 ])  ,
    ṽ_2=([ 5 + 3√(2); 5 - 4√(2); 1 - 5√(2);        -7 ])  .

Thus we find that the thimble
Σ_∞ is different from the original one Σ_∞.

Note, however, that both thimbles are
embedded in a real 3-dimensional hypersurface
defined by
{(x_1 y_1 x_2 y_2)
 | Im S=2x_1y_1-x_2^2+y_2^2=0}∈ℝ^4,
and therefore the sign problem is solved in both cases. 








§ APPLICATION OF THE HMC USING THE ON-AXIS APPROACH




In this section, we briefly review how one can apply

the HMC algorithm to the GTM using the on-axis approach.
The key idea that makes this feasible is the backpropagation
in calculating the force term in the fictitious Hamilton dynamics
of the HMC algorithm <cit.>.
Here we fix the flow time for simplicity, but
it is straightforward to implement
the integration over the flow time,
which is needed in avoiding the ergodicity problem.
(See Ref. <cit.> for the details.)







Let us consider the general partition function given by
(<ref>).
The expectation value of the observable is defined by

  ⟨𝒪⟩   =1/Z
  ∫_ℝ^N dx e^-S(x)𝒪(x)  .
  

In the GTM, this integral is evaluated by

    ⟨𝒪⟩   = 1/Z∫_Σ_τ dz e^-S(z)𝒪(z)  ,
    

where τ is the flow time and
the integration contour Σ_τ is defined by solving
the

gradient flow equation (<ref>).








When one applies the HMC algorithm to evaluate this integral,
one has to define a fictitious Hamilton dynamics to update
the configuration on the deformed contour.
In the on-thimble approach,
one defines the Hamilton dynamics on Σ_τ,
which is a classical mechanics of a constrained system.

Since the deformed contour Σ_τ is not given explicitly
but has to be determined by solving the flow equation,
complicated procedures are necessary in solving the Hamilton equation
in such a way that the configuration is always on Σ_τ.






In the on-axis approach, on the other hand, one rewrites the integral
(<ref>) as

  ⟨𝒪⟩   =
  1/Z∫_ℝ^N dx  
  det  J(x,τ)  e^-S(z(x))𝒪(z(x))  ,

and defines the Hamilton dynamics on the real axis ℝ^N.







Here 







J(x,τ) is the Jacobi matrix
associated with the change of variables
x↦ z(x,τ),

which obeys the following flow equation

∂/∂σJ_kl(x,σ)=H_km(z(x,σ))J_ml(x,σ)  ,

where H_km is the Hessian defined by (<ref>).









The fictitious Hamilton dynamics is defined by
the Hamiltonian

H=1/2  ∑_i=1^N p_i^2 + Re S(z(x,τ))  ,


where p_i

are the momentum variables conjugate to
the coordinate variables x_i.
Note that there is no constraint that complicates
the Hamilton dynamics unlike the on-thimble approach.


In order to obtain
the expectation value of the observables, one uses the reweighting

  ⟨𝒪⟩=⟨𝒪(z(x,τ))   
    det  J(x,τ)   e^-i  Im   S(z(x,τ))⟩_HMC/⟨det  J(x,τ)   e^-i   Im   S(z(x,τ))⟩_HMC      ,


where the expectation value ⟨ · ⟩_ HMC
represents an ensemble average over the configurations
generated by the HMC algorithm with the Hamiltonian (<ref>).




Note that the force term in the Hamilton equation is given by

  F_j   = - ∂ ReS(x,τ)/∂x_j
  =-Re(f_i(x,τ) J_ij(x,τ))  ,


f_i(x,τ)
   =.∂S(z(x,τ))/∂z_i|_z=z(x,τ)  .


If one naively uses this formula to calculate the force term,
the Jacobi matrix J_ij(x,τ) that appears in (<ref>)


requires the computational cost of O(N^2) for a local system
since it involves matrix-matrix products.


In fact, this can be completely avoided by
using the backpropagation <cit.> as we discuss below.









Let us first note that 
z(x,τ) and hence S(z(x,τ)) can be regarded
as functions of z(x,σ) and z̅(x,σ)
(0≤σ≤τ)
since z(x,σ+δσ)≃
z(x,σ)+∂ S(z(x,σ)).
We can therefore define the force at each σ as

  F_j (σ)=∂S(z(x,τ))/∂z_j(x,σ)  ,
      F̅_j (σ)
=∂S(z(x,τ))/∂z̅_j(x,σ)  ,

which obeys the relation

F_j(σ-δσ)
   =F_i(σ)
∂z_i(x,σ)/∂z_j(x, σ-δσ)
+F̅_i(σ)
∂z̅_i(x,σ)/∂z_j(x, σ-δσ)  .



Thus we can calculate F_j =F_j(0) from f_j = F_j(τ)
by solving (<ref>) backward in σ.

Note that
this is the same procedure as the backpropagation used in machine learning.
Since this procedure
involves the matrix-vector multiplication only,
one can calculate the force term
with the computational cost of

O(N) for a local system.












§ BACKPROPAGATION WITH THE PRECONDITIONER




As we have seen in the previous section,

backpropagation (<ref>)
is the key idea that
makes the calculation of the force term (<ref>) 
in the HMC algorithm feasible in the on-axis approach.
In this section,
we discuss how the backpropagation works with the preconditioned
gradient flow in a way compatible with the idea of the multi-mass solver.




In the case of the preconditioned flow (<ref>),
the explicit form of (<ref>) is given by

  F_j(σ-δσ)
     =F_j(σ) + δσ  G_j(σ)   , 

  G_j(σ)
     = 
  F_i(σ)
  (∂_j𝒜_ik)∂_k S+F̅_i(σ)
  [𝒜̅_ikH_kj+(∂_j𝒜̅_ik)
    ∂_kS]  .
  









By using the rational approximation 
(<ref>) for 𝒜,
we can rewrite G_i(σ)

in (<ref>) as


G_j(σ)
=    -∑_q=1^Qa_q(F^⊤(σ)

(H̅H+b_q1)^-1
H̅)_i(∂_jH_ik)(

(H̅H+b_q1)^-1
∂S)_k


+(F̅^⊤𝒜̅H)_j

   -∑_q=1^Qa_q(F̅^⊤(σ)

(H̅H+b_q1)^-1
)_i(∂_jH_ik)(H̅

(H̅H+b_q1)^-1
∂S)_k  . 

As in the calculation of the Jacobian using (<ref>),
the matrix inverse
(H̅H+b_q 1)^-1 in the first line

can be applied to

the two vectors F_j(σ) and ∂ _ j S
separately, and similarly for the matrix inverse in the second line.
Therefore we can avoid the computational cost of O(Q) by 

using the idea of the multi-mass solver.


As we discussed in Section <ref>,

in order to avoid the ergodicity problem in the GTM,
we have to integrate

over the flow time,
which implies that we have to treat τ

as a dynamical variable in the simulation.
In the HMC algorithm, we therefore have to calculate the force term
for τ as described in Section 4 of Ref. <cit.>.
Below we discuss how this can be done efficiently
even in the presence of the preconditioner in the gradient flow.


Since the flow time is discretized as τ = N_τδτ
when we solve the gradient flow equation,

the derivative with respect to τ used in defining the force
F_τ(x,τ) for τ is replaced by the derivative with respect
to the spacing δτ.
The force F_τ(x,τ) for τ is defined by




  F_τ(x,τ)   =
  Re(∂_jS(x,τ=N_τδτ)
    ż_j(x,τ=N_τδτ) )  , 


where ż_j(x,τ=N_τδτ)
is defined by the difference equation

        (n+1)   ż_j(x,(n+1)  δτ) 

    =      n     ż_j(z,n  δτ) + 𝒜_jk
  ∂_kS

      + n   δτ[𝒜_jk   H̅_ki
      ż̅̇_i(x,nδτ) 
    +(∂_i𝒜_jk)   
    ∂_k S   ż_i(x,nδτ)
    +(∂̅_i𝒜_jk)   ∂_k S   
    ż̅̇_i(x,nδτ)] 
  

with the condition ż_j(z,0)=0.
Note that
the second term on the right-hand side

is

an inhomogeneous term, which does not appear
in the flow equation (<ref>)
for the Jacobi matrix.

Nevertheless, we

can calculate (<ref>) by backpropagation as follows.

Using a 2N-dimensional complex vector 
ζ^(n) =(ż(x,nδτ) , ż̅̇(x,nδτ))^⊤,
(<ref>) can be written formally as





ζ^(n+1)=ℬ^(n)ζ^(n)+β^(n)  ,

where ℬ^(n) and β^(n) are a 2N×2N matrix
and a 2N-dimensional vector, respectively, both depending on ζ^(n).
Then, 



the force (<ref>) for the flow time can be expanded as

2F_τ   =V^(N_τ)⊤ ζ^(N_τ)

   =V^(N_τ)⊤(ℬ^(N_τ-1)ζ^(N_τ-1)+β^(N_τ-1))

   =V^(N_τ)⊤ℬ^(N_τ-1)(ℬ^(N_τ-2)ζ^(N_τ-2)+β^(N_τ-2))

+V^(N_τ)⊤β^(N_τ-1)







   = V^(N_τ)⊤ ∑_n=1^N_τ
    ( ∏_m=0^n-1
    ℬ^(m))  β^(N_τ-n)  ,
    

where
we have introduced a 2N-dimensional
vector V^(N_τ)=(v^(N_τ) , v̅^(N_τ))^⊤
with v_j^(N_τ)=∂_jS(x,N_τδτ)
and used the initial condition ζ^(0)=0.
The product ∏_m=0^n-1 in the last line
should be understood as the time-ordered product, 

in which

ℬ^(m) with smaller m appears on the right.


Here we define a vector
V^(n)⊤≡ V^(N_τ)⊤ℬ^(N_τ-1)⋯ℬ^(n),
which can be obtained recursively by using the relation

  V^(n-1)⊤    = V^(n)⊤ ℬ^(n-1)   .
  

Then we can obtain the force (<ref>) by

  2F_τ   = ∑_j=0^N_τ-1 V^(j+1)⊤ β^(j)  .
  
  















  















More explicitly, the algorithm to calculate the force for τ
can be derived as follows.
Let us represent V^(n) as
V^(n)=(v^(n) , v̅^(n))^⊤.






From (<ref>), we obtain 

  v_j^(n)   =n/n+1
  {v_i^(n+1)[δ_ij+δτ(∂_j𝒜_ik)
    ∂_k S]+v̅_i^(n+1)
  δτ[𝒜̅_ikH_kj
    +(∂_j 𝒜_ik)∂_kS]} ,


v_j^(N_τ)   =∂_jS(x,N_τδτ)  .



Then, the force (<ref>)
can be obtained as

  F_τ(x,τ)   =
  Re(f_τ^(0))   ,

where f_τ^(0) is given by solving

f_τ^(n)   =f_τ^(n+1)+1/n+1v_j^(n+1)𝒜_jk
∂_kS  ,


f_τ^(N_τ)   =0  .






From (<ref>),

we find that the computational cost
of O(Q) can be avoided by the multi-mass solver here as well.


















§ THE PARAMETERS OF THE GTM USED IN OUR SIMULATION




In this appendix, we present the parameters chosen for the GTM
such as
the range [τ_ min, τ_ max] of the flow time
as well as m̃, m(τ) and W(τ)
that appear in (<ref>).




In Table <ref>,
we present the choice of

τ_ min, τ_ max, m̃
  and m(τ), which is parameterized as

    m(τ)   =exp(∑_j=0^2a_jτ^j )  .

The symbols O and P in the “flow” column
stand for the original and preconditioned flows, respectively[For
the preconditioned flow with N=20 and x_ f=0, 
the region of τ is chosen to be
[0.4,0.8] in Fig. <ref> for the sake of comparison
with other values of N,
whereas it is chosen to be [0.6,1] in Fig. <ref>,
which turns out to be more optimal.]. 


















In Table <ref>,
we present

the choice of W(τ), which is parameterized as

    W(τ)   =exp(-15τ)+∑_j=1^6b_jτ^j  .

The first term is

introduced to avoid the dominance of τ=0 that occurs otherwise.
One can see that the functions 
m(τ) and W(τ) for the preconditioned flow 
do not have strong dependence on N and x_ f,
which implies that we do not need to fine-tune these functions.

    
    



























JHEP



