2005 Taipei Summer Institute - PowerPoint PPT Presentation

About This 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:

less

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