2005 Taipei Summer Institute - PowerPoint PPT Presentation

About This Presentation

2005 Taipei Summer Institute


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


Transcript and Presenter's Notes

Title: 2005 Taipei Summer Institute

Simulation Algorithms for Lattice QCD
  • A D Kennedy
  • School of Physics, The University of Edinburgh

  • 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

Functional Integrals and QFT
  • The Expectation value of an operator ? is
    defined non-perturbatively by the Functional
  • The action is S (?)
  • Defined in Euclidean space-time
  • Lattice regularisation
  • Continuum limit lattice spacing a ? 0
  • Thermodynamic limit physical volume V ? ?

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

Monte Carlo methods II
Central Limit Theorem I
Central Limit Theorem II
  • Note that this is an asymptotic expansion

Central Limit Theorem III
  • Distribution of the average of N samples
  • Connected generating function

Central Limit Theorem IV
Central Limit Theorem V
Importance Sampling I
Importance Sampling II
Importance Sampling III
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)
Importance Sampling V
  • With 100 rectangles we have V 16.02328561
  • But we can do better!
  • With 100 rectangles we have V 0.011642808

Markov Chains I
  • State space ?
  • (Ergodic) stochastic transitions P ? ? ?
  • Deterministic evolution of probability
    distribution P Q ? Q

Convergence of Markov Chains I
  • The sequence Q, PQ, P²Q, P³Q, is Cauchy

Convergence of Markov Chains II
Convergence of Markov Chains III
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
  • Start with an arbitrary state (field
  • Iterate the Markov process until it has converged
  • 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
  • We can compute ratios of integrals

Markov Chains III
  • Integrate w.r.t. y to obtain fixed point
  • Sufficient but not necessary for fixed point
  • Sufficient but not necessary for detailed

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

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
  • The entire sweep therefore has the desired fixed
    point, and is ergodic
  • But the entire sweep does not satisfy detailed
  • 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

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

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

Autocorrelations II
  • Integrated autocorrelations
  • Consider the autocorrelation of some operator O

Autocorrelations III
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

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

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

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

  • 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

Partial Momentum Refreshment
  • The Gaussian distribution of p is invariant under
  • 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

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.

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

Symplectic Integrators I
  • Baker-Campbell-Hausdorff (BCH) formula

Symplectic Integrators II
  • In order to construct reversible integrators we
    use symmetric symplectic integrators

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

Symplectic Integrators IV
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

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

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

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

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

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

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?)

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

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
  • It does not require the matrix a gt 0, unlike its
    bosonic analogue

Dynamical fermions II
  • Pseudofermions
  • Direct simulation of Grassmann fields is not
  • The problem is not that of manipulating
    anticommuting values in a computer
  • The overall sign of the exponent is unimportant

Dynamical fermions III
  • Any operator ? can be expressed solely in terms
    of the bosonic fields
  • E.g., the fermion propagator is

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

Dynamical fermions V
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

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

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

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

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

  • 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

  • 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
  • 1/x is already a rational function, so RHMC
    reduces to HMC in this case
  • Can be made exact using noisy accept/reject step

Chiral Fermions
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

On-shell chiral symmetry II
Neubergers Operator I
  • We can find a solution of the Ginsparg-Wilson
    relation as follows

Into Five Dimensions
  • H Neuberger hep-lat/9806025
  • A Boriçi hep-lat/9909057, hep-lat/9912040,
  • A Boriçi, A D Kennedy, B Pendleton, U Wenger
  • 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

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

Neubergers Operator III
  • Four dimensional space of algorithms
  • Constraint (5D, 4D)
  • Representation (CF, PF, CTDWF)

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

Constraint I
  • So, what can we do with the Neuberger operator
    represented as a Schur complement?
  • Consider the five-dimensional system of linear

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

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

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

Representation Continued Fraction I
Representation Continued Fraction II
Representation Partial Fraction I
  • Consider a five-dimensional matrix of the form
    (Neuberger Narayanan)

Representation Partial Fraction II
  • Compute its LDU decomposition

Representation Partial Fraction III
Representation Cayley Transform I
  • Compute its LDU decomposition
  • Neither L nor U depend on C

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

Representation Cayley Transform III
Representation Cayley Transform IV
Representation Cayley Transform V
  • With some simple rescaling

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

Chiral Symmetry Breaking
  • G is the quark propagator
  • mres is just one moment of ?L

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

Comparison of Representation
Configuration 806, single precision
Matching mp between HS and HW
Computing mres using ?L
mres per Configuration
Cost versus mres
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
  • Tunnelling between different topological sectors
  • Algorithmic or physical problem (at µ0)
  • Reflection/refraction

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.
Polynomial approximation
  • What is the best polynomial approximation p(x) to
    a continuous function f(x) for x in 0,1 ?

Weierstraß theorem
?????????? polynomials
???????s theorem
  • The error p(x) - f(x) reaches its maximum at
    exactly d2 points on the unit interval

???????s theorem Necessity
  • Suppose p-f has less than d2 extrema of equal
  • 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

???????s theorem Sufficiency
  • Therefore p - p must have d1 zeros on the unit
  • Thus p p 0 as it is a polynomial of degree d

??????? polynomials
  • Convergence is often exponential in d
  • The notation is an old transliteration of ???????

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

??????? 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
  • This appears to be numerically stable.

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

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

Liouvilles theorem
Liouvilles theorem Corollary I
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

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

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

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
The only singularity is a double pole at the
origin and its periodic images
Weierstraß elliptic functions III
Weierstraß elliptic functions IV
Simple pole at the origin with unit residue
Weierstraß elliptic functions V
  • That was so much fun well do it again.

No poles, and simple zero at the origin
Expansion of elliptic functions I
  • Let f be an elliptic function with periods ? and
    ? whose poles and zeros are at aj and ßj

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

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

Addition formulæ II
Jacobian elliptic functions I
Jacobian elliptic functions II
Jacobian elliptic functions III
  • We are hiding the fact that the Weierstraß
    functions depend not only on z but also on
    periods ? and ?

Jacobian elliptic functions IV
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

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

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

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

?????????s formula I
  • Fundamental region
  • Real period 4K, KK(k)
  • Imaginary period 2iK, KK(k) , k2k21

?????????s formula II
Arithmetico-Geometric mean I
  • Long before the work of Jacobi and Weierstraß

Gauß transformation
This is just a special case of the degree n
modular transformation we used before
Arithmetico-Geometric mean II
Arithmetico-Geometric mean III
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
  • These are best done analytically and not
  • With a little care about signs we can also
    compute cn(uk) and dn(uk) which are then also
  • Testing for ab in floating point arithmetic is
    sinful it is better to continue until a ceases
    to decrease

Conclusions II
  • One can compute matrix functions exactly using
    fairly low-order rational approximations
  • Exactly means to within a relative error of 1
  • Which is what exactly normally means for scalar
  • 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

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.

Bibliography II
  • Approximation theory
  • ???? ????? ??????? (Naum Ilich Akhiezer), Theory
    of Approximation, Constable (1956). ISBN
  • Theodore J. Rivlin, An Introduction to the
    Approximation of Functions, Dover (1969). ISBN

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.
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?

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

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
  • Even though the condition number ?(M1/n)?(M)1/n
  • And even though ?(r(M))?(M)1/n
  • So we dont win this way

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

  • 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

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
  • Adjust ? to minimise the cost
  • Easily generalises
  • More than two pseudofermions
  • Wilson-clover action

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
  • This is a good approximation if the condition
    number is dominated by a few isolated tiny
  • 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!

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
Write a Comment
User Comments (0)
About PowerShow.com