# 2005 Taipei Summer Institute - PowerPoint PPT Presentation

Title:

## 2005 Taipei Summer Institute

Description:

### Chiral Fermions. On-shell chiral symmetry. Neuberger's Operator. Into Five ... Chiral Symmetry Breaking. Numerical Studies. Conclusions. Approximation Theory ... – PowerPoint PPT presentation

Number of Views:67
Avg rating:3.0/5.0
Slides: 143
Provided by: tony242
Category:
Tags:
Transcript and Presenter's Notes

Title: 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
• 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
• 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
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
• 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
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
• 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
• Similarly, any elliptic function f can be
expanded in partial fraction form in terms of
Weierstraß ? functions.

117
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
• 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