2005 Taipei Summer Institute

1
Simulation Algorithms for Lattice QCD

• A D Kennedy
• School of Physics, The University of Edinburgh

2
Contents

• Monte Carlo methods
• Functional Integrals and QFT
• Central Limit Theorem
• Importance Sampling
• Markov Chains
• Convergence of Markov Chains
• Autocorrelations
• Hybrid Monte Carlo
• MDMC
• Partial Momentum Refreshment
• Symplectic Integrators
• Dynamical fermions
• Grassmann algebras
• Reversibility
• PHMC
• RHMC

• Chiral Fermions
• On-shell chiral symmetry
• Neubergers Operator
• Into Five Dimensions
• Kernel
• Schur Complement
• Constraint
• Approximation
• tanh
• ?????????
• Representation
• Continued Fraction
• Partial Fraction
• Cayley Transform
• Chiral Symmetry Breaking
• Numerical Studies
• Conclusions

• Approximation Theory
• Polynomial approximation
• Weierstraß theorem
• ?????????? polynomials
• ???????s theorem
• ??????? polynomials
• ??????? approximation
• Elliptic functions
• Liouvilles theorem
• Weierstraß elliptic functions
• Expansion of elliptic functions
• Addition formulæ
• Jacobian elliptic functions
• Modular transformations
• ?????????s formula
• Arithmetico-Geometric mean
• Conclusions
• Bibliography
• Hasenbusch Acceleration

3
Functional Integrals and QFT

• The Expectation value of an operator ? is
  defined non-perturbatively by the Functional
  Integral

• The action is S (?)

• Defined in Euclidean space-time
• Lattice regularisation

• Continuum limit lattice spacing a ? 0
• Thermodynamic limit physical volume V ? ?

4
Monte Carlo methods I

• Monte Carlo integration is based on the
  identification of probabilities with measures
• There are much better methods of carrying out low
  dimensional quadrature
• All other methods become hopelessly expensive for
  large dimensions
• In lattice QFT there is one integration per
  degree of freedom
• We are approximating an infinite dimensional
  functional integral

5
Monte Carlo methods II
6
Central Limit Theorem I
7
Central Limit Theorem II

• Note that this is an asymptotic expansion

8
Central Limit Theorem III

• Distribution of the average of N samples

• Connected generating function

9
Central Limit Theorem IV
10
Central Limit Theorem V
11
Importance Sampling I
12
Importance Sampling II
13
Importance Sampling III
14
Importance Sampling IV
1 Construct cheap approximation to sin(x)
2 Calculate relative area within each rectangle
3 Choose a random number uniformly
4 Select corresponding rectangle
5 Choose another random number uniformly
6 Select corresponding x value within rectangle
7 Compute sin(x)
15
Importance Sampling V

• With 100 rectangles we have V 16.02328561

• But we can do better!

• With 100 rectangles we have V 0.011642808

16
Markov Chains I

• State space ?

• (Ergodic) stochastic transitions P ? ? ?

• Deterministic evolution of probability
  distribution P Q ? Q

17
Convergence of Markov Chains I

• The sequence Q, PQ, P²Q, P³Q, is Cauchy

18
Convergence of Markov Chains II
19
Convergence of Markov Chains III
20
Markov Chains II

• Use Markov chains to sample from Q
• Suppose we can construct an ergodic Markov
  process P which has distribution Q as its fixed
  point
• Start with an arbitrary state (field
  configuration)
• Iterate the Markov process until it has converged
  (thermalized)
• Thereafter, successive configurations will be
  distributed according to Q
• But in general they will be correlated
• To construct P we only need relative
  probabilities of states
• Do not know the normalisation of Q
• Cannot use Markov chains to compute integrals
  directly
• We can compute ratios of integrals

21
Markov Chains III

• Integrate w.r.t. y to obtain fixed point
  condition
• Sufficient but not necessary for fixed point

• Sufficient but not necessary for detailed
  balance

22
Markov Chains IV

• Composition of Markov steps
• Let P1 and P2 be two Markov steps which have the
  desired fixed point distribution
• They need not be ergodic
• Then the composition of the two steps P2?P1 will
  also have the desired fixed point
• And it may be ergodic
• This trivially generalises to any (fixed) number
  of steps
• For the case where P1 is not ergodic but (P1 ) n
  is the terminology weakly and strongly ergodic
  are sometimes used

23
Markov Chains V

• This result justifies sweeping through a
  lattice performing single site updates
• Each individual single site update has the
  desired fixed point because it satisfies detailed
  balance
• The entire sweep therefore has the desired fixed
  point, and is ergodic
• But the entire sweep does not satisfy detailed
  balance
• Of course it would satisfy detailed balance if
  the sites were updated in a random order
• But this is not necessary
• And it is undesirable because it puts too much
  randomness into the system

24
Markov Chains VI

• Coupling from the Past
• Propp and Wilson (1996)
• Use fixed set of random numbers
• Flypaper principle If states coalesce they stay
  together forever
• Eventually, all states coalesce to some state ?
  with probability one
• Any state from t -? will coalesce to ?
• ? is a sample from the fixed point distribution

?
25
Autocorrelations I

• Exponential autocorrelations
• The unique fixed point of an ergodic Markov
  process corresponds to the unique eigenvector
  with eigenvalue 1
• All its other eigenvalues must lie within the
  unit circle

26
Autocorrelations II

• Integrated autocorrelations
• Consider the autocorrelation of some operator O

27
Autocorrelations III
28
Hybrid Monte Carlo I

• In order to carry out Monte Carlo computations
  including the effects of dynamical fermions we
  would like to find an algorithm which
• Update the fields globally
• Because single link updates are not cheap if the
  action is not local
• Take large steps through configuration space
• Because small-step methods carry out a random
  walk which leads to critical slowing down with a
  dynamical critical exponent z2

• Does not introduce any systematic errors

29
Hybrid Monte Carlo II

• A useful class of algorithms with these
  properties is the (Generalised) Hybrid Monte
  Carlo (HMC) method
• Introduce a fictitious momentum p corresponding
  to each dynamical degree of freedom q
• Find a Markov chain with fixed point ?
  exp-H(q,p) where H is the fictitious
  Hamiltonian ½p2 S(q)
• The action S of the underlying QFT plays the
  rôle of the potential in the fictitious
  classical mechanical system
• This gives the evolution of the system in a fifth
  dimension, fictitious or computer time
• This generates the desired distribution exp-S(q)
  if we ignore the momenta p (i.e., the marginal
  distribution)

30
Hybrid Monte Carlo III
Hybrid Monte Carlo III

• The HMC Markov chain alternates two Markov steps
• Molecular Dynamics Monte Carlo (MDMC)
• (Partial) Momentum Refreshment
• Both have the desired fixed point
• Together they are ergodic

31
MDMC I

• If we could integrate Hamiltons equations
  exactly we could follow a trajectory of constant
  fictitious energy
• This corresponds to a set of equiprobable
  fictitious phase space configurations
• Liouvilles theorem tells us that this also
  preserves the functional integral measure dp dq
  as required
• Any approximate integration scheme which is
  reversible and area preserving may be used to
  suggest configurations to a Metropolis
  accept/reject test
• With acceptance probability min1,exp(-?H)

32
MDMC II

• We build the MDMC step out of three parts

• A Metropolis accept/reject step

• with y being a uniformly distributed random
  number in 0,1

33
Partial Momentum Refreshment

• The Gaussian distribution of p is invariant under
  F
• The extra momentum flip F ensures that for small
  ? the momenta are reversed after a rejection
  rather than after an acceptance
• For ? ? / 2 all momentum flips are irrelevant

34
Hybrid Monte Carlo IV

• Special cases
• The usual Hybrid Monte Carlo (HMC) algorithm is
  the special case where ? ? / 2
• ? 0 corresponds to an exact version of the
  Molecular Dynamics (MD) or Microcanonical
  algorithm (which is in general non-ergodic)
• the L2MC algorithm of Horowitz corresponds to
  choosing arbitrary ? but MDMC trajectories of a
  single leapfrog integration step (? ??). This
  method is also called Kramers algorithm.

35
Hybrid Monte Carlo V
Hybrid Monte Carlo V

• Further special cases
• The Langevin Monte Carlo algorithm corresponds to
  choosing ? ? / 2 and MDMC trajectories of a
  single leapfrog integration step (? ??).
• The Hybrid and Langevin algorithms are
  approximations where the Metropolis step is
  omitted
• The Local Hybrid Monte Carlo (LHMC) or
  Overrelaxation algorithm corresponds to updating
  a subset of the degrees of freedom (typically
  those living on one site or link) at a time

36
Symplectic Integrators I

• Baker-Campbell-Hausdorff (BCH) formula

37
Symplectic Integrators II

• In order to construct reversible integrators we
  use symmetric symplectic integrators

38
Symplectic Integrators III

• The basic idea of such a symplectic integrator is
  to write the time evolution operator as

39
Symplectic Integrators IV
40
Symplectic Integrators V

• From the BCH formula we find that the PQP
  symmetric symplectic integrator is given by

• In addition to conserving energy to O (??² ) such
  symmetric symplectic integrators are manifestly
  area preserving and reversible

41
Symplectic Integrators VI

• For each symplectic integrator there exists a
  Hamiltonian H which is exactly conserved
• For the PQP integrator we have

42
Symplectic Integrators VII

• Substituting in the exact forms for the
  operations P and Q we obtain the vector field

43
Symplectic Integrators VIII

• Note that H cannot be written as the sum of a
  p-dependent kinetic term and a q-dependent
  potential term

44
Dynamical fermions I

• Fermion fields are Grassmann valued
• Required to satisfy the spin-statistics theorem
• Even classical Fermion fields obey
  anticommutation relations
• Grassmann algebras behave like negative
  dimensional manifolds

45
Grassmann algebras I

• Grassmann algebras
• Linear space spanned by generators ?1,?2,?3,
  with coefficients a, b, in some field (usually
  the real or complex numbers)
• Algebra structure defined by nilpotency condition
  for elements of the linear space ?² 0
• There are many elements of the algebra which are
  not in the linear space (e.g., ?1?2)
• Nilpotency implies anticommutativity ?? ?? 0
• 0 ?² ?² (? ? )² ?² ?? ?? ?² ??
  ?? 0
• Anticommutativity implies nilpotency, 2?² 0
• Unless the coefficient field has characteristic
  2, i.e., 2 0

46
Grassmann algebras II

• Grassmann algebras have a natural grading
  corresponding to the number of generators in a
  given product
• deg(1) 0, deg(?i ) 1, deg(?i?j ) 2, ...
• All elements of the algebra can be decomposed
  into a sum of terms of definite grading

• A natural antiderivation is defined on a
  Grassmann algebra
• Linearity d(a? b?) a d? b d?
• Anti-Leibniz rule d(??) (d?)? P(?)(d?)

47
Grassmann algebras III

• Definite integration on a Grassmann algebra is
  defined to be the same as derivation

• There is no reason for this function to be
  positive even for real coefficients

48
Grassmann algebras IV

• Gaussian integrals over Grassmann manifolds

• Where Pf(a) is the Pfaffian
• Pf(a)² det(a)

• Despite the notation, this is a purely algebraic
  identity
• It does not require the matrix a gt 0, unlike its
  bosonic analogue

49
Dynamical fermions II

• Pseudofermions
• Direct simulation of Grassmann fields is not
  feasible
• The problem is not that of manipulating
  anticommuting values in a computer

• The overall sign of the exponent is unimportant

50
Dynamical fermions III

• Any operator ? can be expressed solely in terms
  of the bosonic fields

• E.g., the fermion propagator is

51
Dynamical fermions IV

• The determinant is extensive in the lattice
  volume, thus again we get poor importance sampling

• The fermion kernel must be positive definite (all
  its eigenvalues must have positive real parts)
  otherwise the bosonic integral will not converge
• The new bosonic fields are called pseudofermions

52
Dynamical fermions V
53
Dynamical fermions VI

• It is not necessary to carry out the inversions
  required for the equations of motion exactly
• There is a trade-off between the cost of
  computing the force and the acceptance rate of
  the Metropolis MDMC step
• The inversions required to compute the
  pseudofermion action for the accept/reject step
  does need to be computed exactly, however
• We usually take exactly to by synonymous with
  to machine precision

54
Reversibility I

• Are HMC trajectories reversible and area
  preserving in practice?
• The only fundamental source of irreversibility is
  the rounding error caused by using finite
  precision floating point arithmetic
• For fermionic systems we can also introduce
  irreversibility by choosing the starting vector
  for the iterative linear equation solver
  time-asymmetrically
• We do this if we to use a Chronological Inverter,
  which takes (some extrapolation of) the previous
  solution as the starting vector
• Floating point arithmetic is not associative
• It is more natural to store compact variables as
  scaled integers (fixed point)
• Saves memory
• Does not solve the precision problem

55
Reversibility II

• Data for SU(3) gauge theory and QCD with heavy
  quarks show that rounding errors are amplified
  exponentially
• The underlying continuous time equations of
  motion are chaotic
• ??????? exponents characterise the divergence of
  nearby trajectories
• The instability in the integrator occurs when ?H
  1
• Zero acceptance rate anyhow

56
Reversibility III

• In QCD the ??????? exponents appear to scale
  with ? as the system approaches the continuum
  limit ? ? ?
• ?? constant
• This can be interpreted as saying that the
  ??????? exponent characterises the chaotic nature
  of the continuum classical equations of motion,
  and is not a lattice artefact
• Therefore we should not have to worry about
  reversibility breaking down as we approach the
  continuum limit
• Caveat data is only for small lattices, and is
  not conclusive

57
Reversibility IV

• Data for QCD with lighter dynamical quarks
• Instability occurs close to region in ?? where
  acceptance rate is near one
• May be explained as a few modes becoming
  unstable because of large fermionic force
• Integrator goes unstable if too poor an
  approximation to the fermionic force is used

58
PHMC

• Polynomial Hybrid Monte Carlo algorithm instead
  of using ??????? polynomials in the multiboson
  algorithm, Frezzotti Jansen and deForcrand
  suggested using them directly in HMC
• Polynomial approximation to 1/x are typically of
  order 40 to 100 at present
• Numerical stability problems with high order
  polynomial evaluation
• Polynomials must be factored
• Correct ordering of the roots is important
• Frezzotti Jansen claim there are advantages
  from using reweighting

59
RHMC

• Rational Hybrid Monte Carlo algorithm the idea
  is similar to PHMC, but uses rational
  approximations instead of polynomial ones
• Much lower orders required for a given accuracy
• Evaluation simplified by using partial fraction
  expansion and multiple mass linear equation
  solvers
• 1/x is already a rational function, so RHMC
  reduces to HMC in this case
• Can be made exact using noisy accept/reject step

60
Chiral Fermions
61
On-shell chiral symmetry I

• It is possible to have chiral symmetry on the
  lattice without doublers if we only insist that
  the symmetry holds on shell

62
On-shell chiral symmetry II
63
Neubergers Operator I

• We can find a solution of the Ginsparg-Wilson
  relation as follows

64
Into Five Dimensions

• H Neuberger hep-lat/9806025
• A Boriçi hep-lat/9909057, hep-lat/9912040,
  hep-lat/0402035
• A Boriçi, A D Kennedy, B Pendleton, U Wenger
  hep-lat/0110070
• R Edwards U Heller hep-lat/0005002
• ? ? ? (T-W Chiu) hep-lat/0209153,
  hep-lat/0211032, hep-lat/0303008
• R C Brower, H Neff, K Orginos hep-lat/0409118
• Hernandez, Jansen, Lüscher hep-lat/9808010

65
Neubergers Operator II

• Is DN local?
• It is not ultralocal (Hernandez, Jansen, Lüscher)
• It is local iff DW has a gap
• DW has a gap if the gauge fields are smooth
  enough
• It seems reasonable that good approximations to
  DN will be local if DN is local and vice versa
• Otherwise DWF with n5 ? 8 may not be local

66
Neubergers Operator III

• Four dimensional space of algorithms

• Constraint (5D, 4D)

• Representation (CF, PF, CTDWF)

67
Kernel
68
Schur Complement

• It may be block diagonalised by an LDU
  factorisation (Gaussian elimination)

69
Constraint I

• So, what can we do with the Neuberger operator
  represented as a Schur complement?

• Consider the five-dimensional system of linear
  equations

70
Constraint II

• and we are left with just det Dn,n det DN

71
Approximation tanh

• Pandey, Kenney, Laub Higham Neuberger
• For even n (analogous formulæ for odd n)

?j
72
Approximation ?????????
?j
73
Approximation Errors

• The fermion sgn problem
• Approximation over 10-2 lt x lt 1
• Rational functions of degree (7,8)

74
Representation Continued Fraction I
75
Representation Continued Fraction II
76
Representation Partial Fraction I

• Consider a five-dimensional matrix of the form
  (Neuberger Narayanan)

77
Representation Partial Fraction II

• Compute its LDU decomposition

78
Representation Partial Fraction III
79
Representation Cayley Transform I

• Compute its LDU decomposition

• Neither L nor U depend on C

80
Representation Cayley Transform II

• In Minkowski space a Cayley transform maps
  between Hermitian (Hamiltonian) and unitary
  (transfer) matrices

81
Representation Cayley Transform III
µP
µP-
P-
P
82
Representation Cayley Transform IV
83
Representation Cayley Transform V

• With some simple rescaling

84
Representation Cayley Transform VI

• It therefore appears to have exact off-shell
  chiral symmetry
• But this violates the Nielsen-Ninomiya theorem
• q.v., Pelissetto for non-local version
• Renormalisation induces unwanted ghost doublers,
  so we cannot use DDW for dynamical (internal)
  propagators
• We must use DN in the quantum action instead
• We can us DDW for valence (external)
  propagators, and thus use off-shell (continuum)
  chiral symmetry to manipulate matrix elements

85
Chiral Symmetry Breaking

• G is the quark propagator

• mres is just one moment of ?L

86
Numerical Studies

• Matched p mass for Wilson and Möbius kernels
• All operators are even-odd preconditioned
• Did not project eigenvectors of HW

87
Comparison of Representation
Configuration 806, single precision
88
Matching mp between HS and HW
89
Computing mres using ?L
90
mres per Configuration
e
91
Cost versus mres
92
Conclusions I

• Relatively good
• Zolotarev Continued Fraction
• Rescaled Shamir DWF via Möbius (tanh)
• Relatively poor (so far)
• Standard Shamir DWF
• Zolotarev DWF (Chiu)
• Can its condition number be improved?
• Still to do
• Projection of small eigenvalues
• HMC
• 5 dimensional versus 4 dimensional dynamics
• Hasenbusch acceleration
• 5 dimensional multishift?
• Possible advantage of 4 dimensional nested Krylov
  solvers
• Tunnelling between different topological sectors
• Algorithmic or physical problem (at µ0)
• Reflection/refraction

93
Approximation Theory
We review the theory of optimal polynomial and
rational ??????? approximations, and the theory
of elliptic functions leading to ?????????s
formula for the sign function over the range Rz
e z1. We show how Gauß arithmetico-geometr
ic mean allows us to compute the ?????????
coefficients on the fly as a function of e. This
allows us to calculate sgn(H) quickly and
accurately for a Hermitian matrix H whose
spectrum lies in R.
94
Polynomial approximation

• What is the best polynomial approximation p(x) to
  a continuous function f(x) for x in 0,1 ?

95
Weierstraß theorem
96
?????????? polynomials
97
???????s theorem

• The error p(x) - f(x) reaches its maximum at
  exactly d2 points on the unit interval

98
???????s theorem Necessity

• Suppose p-f has less than d2 extrema of equal
  magnitude

• Then at most d1 maxima exceed some magnitude

• This defines a gap

Lagrange was born in Turin!

• And whose magnitude is smaller than the gap

• The polynomial pq is then a better approximation
  than p to f

99
???????s theorem Sufficiency

• Therefore p - p must have d1 zeros on the unit
  interval

• Thus p p 0 as it is a polynomial of degree d

100
??????? polynomials

• Convergence is often exponential in d

• The notation is an old transliteration of ???????
  !

101
??????? rational functions

• ???????s theorem is easily extended to rational
  approximations
• Rational functions with nearly equal degree
  numerator and denominator are usually best
• Convergence is still often exponential
• Rational functions usually give a much better
  approximation

102
??????? rationals Example

• This is accurate to within almost 0.1 over the
  range 0.003,1

• Using a partial fraction expansion of such
  rational functions allows us to use a multishift
  linear equation solver, thus reducing the cost
  significantly.

• This appears to be numerically stable.

103
Polynomials versus rationals

• This has L2 error of O(1/n)

• Optimal L8 approximation cannot be too much
  better (or it would lead to a better L2
  approximation)

104
Elliptic functions

• Elliptic function are doubly periodic analytic
  functions
• Jacobi showed that if f is not a constant it can
  have at most two primitive periods

105
Liouvilles theorem
106
Liouvilles theorem Corollary I
107
Liouvilles theorem Corollary II

• It follows that there are no non-constant
  holomorphic elliptic functions
• If f is an analytic elliptic function with no
  poles then f(z)-a could have no zeros either

108
Liouvilles theorem Corollary III

• Where n and n and the number of times f(z) winds
  around the origin as z is taken along a straight
  line from 0 to ? or ? respectively

109
Weierstraß elliptic functions I

• A simple way to construct a doubly periodic
  analytic function is as the convergent sum

• Note that Q is an odd function

110
Weierstraß elliptic functions II
The name of the function is ?, and the function
is called the Weierstraß function, but what is
its name called? (with apologies to the White
Knight)
The only singularity is a double pole at the
origin and its periodic images
111
Weierstraß elliptic functions III
112
Weierstraß elliptic functions IV
Simple pole at the origin with unit residue
113
Weierstraß elliptic functions V

• That was so much fun well do it again.

No poles, and simple zero at the origin
114
Expansion of elliptic functions I

• Let f be an elliptic function with periods ? and
  ? whose poles and zeros are at aj and ßj
  respectively

115
Expansion of elliptic functions II

• By Liouvilles theorem f/g is an elliptic
  function with no poles, and is thus a constant C

116
Addition formulæ I

• Similarly, any elliptic function f can be
  expanded in partial fraction form in terms of
  Weierstraß ? functions.

117
Addition formulæ II
118
Jacobian elliptic functions I
119
Jacobian elliptic functions II
120
Jacobian elliptic functions III

• We are hiding the fact that the Weierstraß
  functions depend not only on z but also on
  periods ? and ?

121
Jacobian elliptic functions IV
122
Modular transformations I

• Since the period lattices are the same, elliptic
  functions with these periods are rationally
  expressible in terms of each other

• This is called a first degree transformation

123
Modular transformations II

• Elliptic functions with these periods must also
  be rationally expressible in terms of ones with
  the original periods, although the inverse may
  not be true

• This is called a transformation of degree n

124
Modular transformations III

• This is easily done by taking sn(zk) and scaling
  z by some factor 1/M and choosing a new
  parameter ?

• The periods for z are thus 4LM and 2iLM, so we
  must have LMK and LMK/n

125
Modular transformations IV
126
Modular transformations V

• Likewise, the parameter ? is found by evaluating
  the identity at z K iK/n

127
?????????s formula I

• Fundamental region

• Real period 4K, KK(k)
• Imaginary period 2iK, KK(k) , k2k21

128
?????????s formula II
129
Arithmetico-Geometric mean I

• Long before the work of Jacobi and Weierstraß

130
Gauß transformation
This is just a special case of the degree n
modular transformation we used before
131
Arithmetico-Geometric mean II
132
Arithmetico-Geometric mean III
133
Arithmetico-Geometric mean IV

• We have a rapidly converging recursive function
  to compute sn(uk) and K(k) for real u and 0ltklt1
• Values for arbitrary complex u and k can be
  computed from this using simple modular
  transformations
• These are best done analytically and not
  numerically
• With a little care about signs we can also
  compute cn(uk) and dn(uk) which are then also
  needed
• Testing for ab in floating point arithmetic is
  sinful it is better to continue until a ceases
  to decrease

134
Conclusions II

• One can compute matrix functions exactly using
  fairly low-order rational approximations
• Exactly means to within a relative error of 1
  u.l.p.
• Which is what exactly normally means for scalar
  computations
• One can also use an absolute error, but it makes
  no difference for the sign function!
• For the case of the sign function we can evaluate
  the ????????? coefficients needed to obtain this
  precision on the fly just to cover the spectrum
  of the matrix under consideration
• There seem to be no numerical stability problems

135
Bibliography I

• Elliptic functions
• ???? ????? ??????? (Naum Ilich Akhiezer),
  Elements of the Theory of Elliptic Functions,
  Translations of Mathematical Monographs, 79, AMS
  (1990). ISBN 0065-9282 QA343.A3813.

136
Bibliography II

• Approximation theory
• ???? ????? ??????? (Naum Ilich Akhiezer), Theory
  of Approximation, Constable (1956). ISBN
  0094542007.
• Theodore J. Rivlin, An Introduction to the
  Approximation of Functions, Dover (1969). ISBN
  0-486-64069-8.

137
Hasenbusch Acceleration
Hasenbusch introduced a simple method that
significantly improves the performance of
dynamical fermion computations. We shall consider
a variant of it that accelerates fermionic
molecular dynamics algorithms by introducing n
pseudofermion fields coupled with the nth root of
the fermionic kernel.
138
Non-linearity of CG solver

• Suppose we want to solve A2xb for Hermitian A
  by CG
• It is better to solve Axy, Ayb successively
• Condition number ?(A2) ?(A)2
• Cost is thus 2?(A) lt ?(A2) in general
• Suppose we want to solve Axb
• Why dont we solve A1/2xy, A1/2yb successively?
• The square root of A is uniquely defined if Agt0
• This is the case for fermion kernels
• All this generalises trivially to nth roots
• No tuning needed to split condition number evenly
• How do we apply the square root of a matrix?

139
Rational matrix approximation

• Functions on matrices
• Defined for a Hermitian matrix by diagonalisation
• H UDU-1
• f(H) f(UDU-1) U f(D) U-1
• Rational functions do not require diagonalisation
• ? Hm ? Hn U(? Dm ? Dn) U-1
• H-1 UD-1U-1
• Rational functions have nice properties
• Cheap (relatively)
• Accurate

140
No Free Lunch Theorem

• We must apply the rational approximation with
  each CG iteration
• M1/n ? r(M)
• The condition number for each term in the partial
  fraction expansion is approximately ?(M)
• So the cost of applying M1/n is proportional to
  ?(M)
• Even though the condition number ?(M1/n)?(M)1/n
• And even though ?(r(M))?(M)1/n
• So we dont win this way

141
Pseudofermions

• We want to evaluate a functional integral
  including the fermionic determinant det M

142
Multipseudofermions

• We are introducing extra noise into the system by
  using a single pseudofermion field to sample this
  functional integral
• This noise manifests itself as fluctuations in
  the force exerted by the pseudofermions on the
  gauge fields
• This increases the maximum fermion force
• This triggers the integrator instability
• This requires decreasing the integration step size

143
Hasenbuschs method

• Clever idea due to Hasenbusch
• Start with the Wilson fermion action M1-?H
• Introduce the quantity M1-?H
• Use the identity M M(M-1M)
• Write the fermion determinant as det M det M
  det (M-1M)
• Introduce separate pseudofermions for each
  determinant
• Adjust ? to minimise the cost
• Easily generalises
• More than two pseudofermions
• Wilson-clover action

144
Violation of NFL Theorem

• So lets try using our nth root trick to
  implement multipseudofermions
• Condition number ?(r(M))?(M)1/n
• So maximum force is reduced by a factor of
  n?(M)(1/n)-1
• This is a good approximation if the condition
  number is dominated by a few isolated tiny
  eigenvalues
• This is so in the case of interest
• Cost reduced by a factor of n?(M)(1/n)-1
• Optimal value nopt ? ln ?(M)
• So optimal cost reduction is (e ln?) /?
• This works!

145
RHMC technicalities

• In order to keep the algorithm exact to machine
  precision (as a Markov process)
• Use a good (machine precision) rational
  approximation for the pseudofermion heatbath
• Use a good rational approximation for the HMC
  acceptance test at the end of each trajectory
• Use as poor a rational approximation as we can
  get away with to compute the MD force