On the Fermion Sign Problem in Imaginary-Time Projection Continuum Quantum Monte Carlo with Local Interaction

0.85 MB·

by  Francesco Calcavecchia

On the Fermion Sign Problem in Imaginary-Time Projection Continuum Quantum
Monte Carlo with Local Interaction
Francesco Calcavecchia∗
LPMMC, UMR 5493, Boîte Postale 166, 38042 Grenoble, France
Institute of Physics, Johannes Gutenberg-University, Staudingerweg 7, D-55128 Mainz, Germany and
Graduate School of Excellence Materials Science in Mainz, Staudingerweg 9, D-55128 Mainz, Germany
Markus Holzmann†
LPMMC, UMR 5493 of CNRS, Université Grenoble Alpes, 38042 Grenoble, France and
Institut Laue Langevin, BP 156, F-38042 Grenoble Cedex 9, France
6
1 We use the Shadow Wave Function formalism as a convenient model to study the fermion sign
0 problemaffectingallprojectorQuantumMonteCarlomethodsincontinuumspace. Wedemonstrate
2 thattheefficiencyofimaginarytimeprojectionalgorithmsdecaysexponentiallywithincreasingnum-
ber of particles and/or imaginary-time propagation. Moreover, we derive an analytical expression
r
p that connects the localization of the system with the magnitude of the sign problem, illustrating
A this behavior through numerical results. Finally, we discuss the computational complexity of the
fermion sign problem and methods for alleviating its severity.
6

] INTRODUCTION The fact that it has not been possible to solve not even
h
p oneofthem, justifiesthecommonbeliefthattheseprob-
- Thefermion sign problem isoneofthemostrenowned lem are intractable. As a matter of fact, NPC problems
p
open problems in computational physics. It consists in are so firmly believed to be intractable, that all classi-
m
finding a general algorithm able to determine the exact calencryptionschemesrelyonthisconjecture. However,
o fermionic ground state with a computational cost that despite the importance of this conjecture (known in the
c growsatmostpolynomiallywiththenumberofsimulated literatureasNP(cid:54)=Phypothesis), aproofisstillmissing.
.
s particles. In fact, all known exact algorithms for classi- In2005,TroyerandWieseprovidedademonstrationof
c
i calcomputersscaleexponentiallyexceptforasmallclass the NP-hardness of the Monte Carlo (MC) sign problem
s of model systems [1–9]. In particular, Quantum Monte [15] for a specific Hamiltonian. The NP-hard problems
y
Carlo (QMC) methods [10–13], which are able to pro- are a class of problems which are mappable in polyno-
h
p vide the exact bosonic ground state in polynomial time mialtimeintoNPCsothatsolvinganyNP-hardproblem
[ forawiderangeofHamiltonians,sufferfromasignprob- would provide the solution to all NPC problems. In this
lem when applied to fermions. In the following we will sense NP-hard problems are said to be the hardest ones,
2
often refer to the solution of the fermion sign problem, as they are at least as hard as any other NPC problem
v
8 implicitlymeaningthattheexactfermiongroundstateis (see Fig 1).
5 obtained in polynomial time rather than in exponential
5 time with respect to the number of particles.
1 In complexity theory, the class of decision problems
0
solvableinpolynomialtimeonadeterministic(classical)
.
1 machine are called P (deterministic Polynomial time),
0 whereasproblemswhichcanbeefficientlysolvedbyprob-
6 abilisticalgorithmsarecalledBPP(Bounded-errorProb-
1
abilistic Polynomial time). In respect to this terminol-
:
v ogy,generalfermionicsimulationssufferingfromthesign
i problemseemtoremainoutsidetheP/BPPclasses. This
X
is a reminiscence of Non-deterministic Polynomial Com-
r
a plete (NPC) problems, a set of hundreds of apparently
different problems that, despite many efforts, have not
been solved yet. The most famous of such problems is
the Travelling Salesman problem, formulated in 1930.
Notably, it has been found that all these problems are
mappable one into the other, so that the solution of one
of them would imply the solution of all of them [14].
Figure 1. Schematic representation of computational com-
plexity of QMC simulations.
∗ [email protected] Inthepresentpaper,wefocusonfermionicQMCsimu-
† [email protected] lationsofcontinuumsystemswithlocalinteractions, e.g.
2
non-relativistic electrons interacting via static Coulomb I. THE SHADOW WAVE FUNCTION
potential,anddiscussindetailthecorrespondingfermion FORMULATION
signproblem. Inparticular,westudytheperformanceof
the imaginary-time projection techniques for these sys- The Shadow Wave Function (SWF) is a class of varia-
tems, which provide polynomial scaling solutions to the tional trial wave functions which, by embedding an inte-
correspondingbosonicproblems. Wegiveageneralproof gral,canprofitfromagreatflexibility[16,31]. Ingeneral,
of the exponential scaling of the efficiency of such al- the SWF can be written as
gorithms, both in fermion number and projection time. (cid:90)
Further, we explicitly show that localized orbitals can Ψ (R)=J (R) dSΞ(R,S)ϕ (S)
SWF R T
drastically reduce the sign problem. Our discussion is (2)
(cid:90)
based on the Shadow Wave Function (SWF) formalism
= dSϕ (R,S)
[16]. If the fermion sign problem can be solved for the SWF
SWF, such a solution will be extendable to all the other
where R ≡ (r ,r ,...r ) represents the particle coor-
QMC methods. Whereas, if SWF can be proved to fall 1 2 N
dinates, S labels some 3N-dimensional auxiliary coordi-
in the NP-hard class, the same will apply to the other
nates,J isaJastrow,ϕ isachosentrialwavefunction,
imaginary-time projections QMC methods. R T
and
Weremarkthat,evenifaproblemisNP-hard,onecan Ξ(R,S)=e−C(R−S)2 (3)
obtain approximate solutions, and work on the improve-
ment of the approximations. For example, even though istheso-calledkernel. TheSWFefficientlyremoveslarge
the Travelling Salesman problem cannot be solved ex- part of the bias introduced by the underlying trial wave
actly for many problems of practical interest, there are function, in particular close to phase transitions or in
several methods that provide excellent approximations inhomogeneous systems [32–37]
[17, 18]. Concerning the fermion sign problem, remark- The application of the SWF to fermions requires the
able progress has been obtained over the recent years fulfilment of the Fermi-Dirac statistics by the use of an
by employing more flexible trial wave functions and im- antisymmetric form for spin-like particles. In the follow-
provedoptimizationprocedures,whichsystematicallyre- ing we split the trial wave function into the product of a
duce the deviations from the exact ground state [19–22]. symmetrical part, the Jastrow factor, times an antisym-
metrical part, usually taken as a Slater determinant.
Adifferentstrategywhichwefollowinthispaperisto
A straightforward antisymmetrization of Eq. (2) gives
alleviatethesignproblem,whereoneaimestoreducethe
the Fermionic Shadow Wave Function (FSWF) form
pre-factoroftheultimateexponentialdecayinthesignal
of the desired quantities, similar to released-node algo- (cid:90)
Ψ (R)=J (R) dSΞ(R,S)J (S)Φ(S)
rithm [23] or different exact fermion simulations [24–30]. FSWF R S
Inordertocomparedifferentapproachesinthefollowing, (cid:90) (4)
we refer to the Monte Carlo efficiency, defined as = dSϕ (R,S)
FSWF
where Φ represents a Slater determinant and J a Jas-
S
1 trow. Unfortunately, the FSWF introduces a sign prob-
η ≡ , (1)
cpu time×var lem [38] which emerges by the fact that the product
ϕ∗ (R,S )×ϕ (R,S ) is not necessarily positive.
FSWF 1 FSWF 2
The sign problem of FSWF is avoided using the Anti-
where var is the variance of the computed quantity. The symmetric Shadow Wave Function (ASWF)
SWF is a perfect tool-bench for a systematic study of
(cid:90)
the efficiency of different strategies, as the severity of Ψ (R)=J (R)Φ(R) dSΞ(R,S)J (S)
ASWF R S
thesignproblemcanbecontrolledbykernelparameters,
(5)
(cid:90)
and the method itself is considerably cheaper than most
= dSϕ (R,S),
of imaginary-time projection methods. ASWF
This paper is structured as follows. In Sec. I we in- sinceϕ∗ (R,S )×ϕ (R,S )isgarantiedtobepos-
ASWF 1 ASWF 2
troduce the Shadow Wave Function formalism, and then itive for any R, S and S .
1 2
explain its connection with Diffusion Monte Carlo and Typically, FSWF is considered of higher quality and
Path Integral Ground State in Sec. II. We will charac- expected to yield lower variational energies than ASWF.
terizeitsfermionsignproblemwhenusingaSlaterdeter-
minant of simple plane waves in Sec. III, and generalize
thisresulttoorbitalsofanykindinSec. IV.Wewillthen II. CONNECTION BETWEEN SWF, DMC, AND
review some general method that have been proposed to PIGS
tackle the fermion sign problem in Sec. V, and finally
discuss about its computational complexity in Sec. VI. The SWF can be regarded as a prototype method for
The conclusion are drawn in Sec. VII. any kind of imaginary-time projection method, since the
3
kernel has the form of an approximated Green’s func- In contrast to projector Monte Carlo methods, e.g.
tion with C proportional to the inverse imaginary-time Diffusion Monte Carlo (DMC) or Path Integral Ground
propagation τ. State Monte Carlo (PIGS), SWF can be considered as a
Let us consider a trial wave-function ϕ which is not single step of a chain of small imaginary-time propaga-
T
orthogonal to the exact ground state. It is well-known tions. As a consequence, SWF is in general not exact,
that a propagation in imaginary time will eventually even though it often captures most of the corrections of
project it to the exact ground state, i.e. the imaginary-time propagation, while retaining a low-
computational cost. Further, SWF remains a explicit
e−τH |Ψ (cid:105)−τ−→−−∞→|Ψ (cid:105) (6) trial wave function subject to the Rayleigh-Ritz varia-
T GS
tional principle, and δτ (corresponding to C) can be re-
In order to build a concrete algorithm built upon such
gardedasavariationalparameteratvariancetoprojector
property,wemakeuseoftheSuzuki-Trotterformulaand
MonteCarlomethodswhichmustbeextrapolatedtothe
write
limit δτ →0.
e−τ(T+V) (cid:39)e−τT e−τV , (7)
which is exact in the limit τ →0. III. THE SIGN PROBLEM OF THE
In order to conciliate the necessity of having a large τ FERMIONIC SHADOW WAVE FUNCTION
(Eq. (6)) with the Suzuki-Trotter approximation which
requires a small τ, we can break the propagation into ThesignproblemoftheFermionicShadowWaveFunc-
several small ones. Given a N >1 such that that δτ ≡ tion has been explored numerically in Ref. [38]. In this
τ
τ/N issmallenoughfortheapproximationinEq. (7)to sectionwearegoingtointroducetwodifferentapproaches
τ
hold,theoperatorse−δτTe−δτV canbeappliediteratively which justify qualitatively and quantitatively its occur-
tothetrialwavefunctiontoprojectoutthegroundstate rence.
wave function WewillassumethattheorbitalsoftheSlaterdetermi-
nant are simple plane waves throughout this whole sec-
|Ψ (cid:105)(cid:39)|Ψ [τ](cid:105)≡e−τH |Ψ [0](cid:105)
GS T T tion and generalize the results to any kind of orbitals in
(cid:32)(cid:89)Nτ (cid:33) (8) Sec. IV.
(cid:39) e−δτTe−δτV |Ψ [0](cid:105)
T
i=1
A. Ratio with a positive-definite distribution
The SWF has the functional form of single propagation
δτ, since the kinetic energy propagator is a diffusor term
The expectation value of an operator O computed av-
in the coordinate space, i.e.
eraging over a Shadow Wave Function writes
(cid:104)R|e−δτT |R(cid:48)(cid:105)∝e−(R(cid:48)4−δτR)2 (9) (cid:82) dRdS dS ϕ∗ (R,S )Oϕ (R,S )
(cid:104)O(cid:105)= (cid:82) 1 2 SWF 1 SWF 2
dRdS dS ϕ∗ (R,S )ϕ (R,S )
Using a variation of the Suzuki-Trotter approximation 1 2 SWF 1 SWF 2
which splits the potential propagator on the left and on (cid:82) dRdS dS ϕ∗ (R,S )ϕ (R,S )O (R,S )
the right symmetrically = (cid:82) d1RdS2 dSSWFϕ∗ 1(R,SSWF)ϕ 2(R,SL ) 2
1 2 SWF 1 SWF 2
(13)
e−δτ(T+V) (cid:39)e−δ2τV e−τT e−δ2τV , (10)
where
we can write (assuming that the potential V is real and
Oϕ (R,S)
local) O (R,S)≡ SWF . (14)
L ϕ (R,S)
SWF
ϕ [δτ](R)=(cid:104)Ψ [δτ]|R(cid:105)
T T If we denote by ρ(R,S ,S ) the probability density
=(cid:104)Ψ [0]|e−δ2τV e−δτT e−δ2τV |R(cid:105) function (pdf) that we inte1nd2to sample from, and intro-
T
(cid:90) duce the corresponding weight
= dS(cid:104)Ψ [0]|S(cid:105)(cid:104)S |e−δ2τV e−δτT e−δ2τV |R(cid:105)
T
ϕ∗ (R,S )ϕ (R,S )
(cid:90) w(R,S ,S )= SWF 1 SWF 2 , (15)
= dSϕ (S)e−δ2τV(S)e−(R4−δSτ)2 e−δ2τV(R) 1 2 ρ(R,S1,S2)
T
(11) then we can recast Eq. (13) as
(cid:82)
The generic form of the SWF (Eq. (2)) is the one ob- (cid:104)O(cid:105)= dR(cid:82)dS1dS2ρ(R,S1,S2)w(R,S1,S2)OL(R,S2).
tained by mapping dRdS dS ρ(R,S ,S )w(R,S ,S )
1 2 1 2 1 2
(16)
ϕ4δ1Tτ(S(cid:55)→)eC−δ2τV(S) (cid:55)→ϕT(S) (12) defiIfntithee, pwreodcuacntcϕh∗SoWsFe(Rρ,sSu1c)h×thϕaStWwF(R=,S12)anisdpnoositsiivgen-
e−δ2τV(R) (cid:55)→JR(R) problem will occur.
4
However, if this is not the case, a typical choice is so that we can integrate out dS dS , and obtain
1 2
ρ(R,S1,S2)≡|ϕ∗SWF(R,S1)ϕSWF(R,S2)| (17) e−(cid:80)Ni=4C1k2i (cid:82) dR det(e−ikα·rβ) det(eikα·rβ)
and therefore (cid:104)w(cid:105)FSWF (cid:39) (cid:82) dR det(e−ikα·rβ) det(eikα·rβ)
w(R,S1,S2)=sign(ϕ∗SWF(R,S1)ϕSWF(R,S2))=±1(1.8) =e−(cid:80)Ni=4C1k2i .
introducing a sign problem. The expectation value of O (25)
is then
In the thermodynamic limit we then get
(cid:104)wO (cid:105)
(cid:104)O(cid:105)= L ρ (19)
(cid:104)w(cid:105)ρ (cid:104)w(cid:105)FSWF ∝e−NCρ2/3. (26)
whereby(cid:104)...(cid:105) wemeantheaverageresultingfromsam-
ρ where ρ is the density. Therefore, the efficiency of a
pling the pdf ρ.
QMCsimulationforcomputing(cid:104)O(cid:105)employingtheFSWF
Let us now focus on (cid:104)w(cid:105) .
ρ writes
In the case of PIGS with a projection time τ or Path-
Integral Monte Carlo at finite temperature T = 1/τ, η ∝e−NCρ2/3. (27)
(cid:104)w(cid:105) is equal to the ratio between the fermionic and the
ρ
Assumptionsiandiimayberelaxedforsituationswhere
bosonic partition functions [39] (where the bosonic sys-
reweighting is possible, but in subsection IIIB we will
tem is defined by the positive-definite weight |w|) and
derive the same scaling without relying on them at all.
therefore
(cid:104)w(cid:105) =e−τN∆F (20)
ρ
B. Difference with a positive-definite distribution
where N is the number of particles and ∆F ≥ 0 is the
free energy difference per particle between the fermionic
and the bosonic system. Relation (20) explains the ex- In the spirit of the control variates technique [40, 41],
ponential decay in efficiency, as we recast the FSWF as
(cid:115) (cid:114) (cid:90)
σ((cid:104)w(cid:105)ρ) = (cid:104)w2(cid:105)ρ−(cid:104)w(cid:105)2ρ = 1/(cid:104)w(cid:105)2ρ−1 ∼ eτ√N∆F ΨFSWF(R)= dS (ϕFSWF(R,S)−ϕ˜FSWF(R,S))
(cid:104)w(cid:105)ρ M(cid:104)w(cid:105)2ρ M M (cid:90) (28)
(21) + dSϕ˜ (R,S),
FSWF
whereM isthenumberofsampledpoints. Sincetherela-
tiveerrorof(cid:104)O(cid:105)isgivenbythesumoftherelativeerrors with
of(cid:104)wO (cid:105) and(cid:104)w(cid:105) ,wecanseethatEq. (21)issufficient
to explLainρ the expρonential decay of the efficiency of any ϕ˜ (R,S)=J (R)J˜(R)Ξ(R,S)Φ(S) (29)
FSWF R S
observable with N and τ.
We can now chose the local normalization factor, J˜(R),
In the case of SWF, we can evaluate (cid:104)w(cid:105) explicitly, S
ρ such that
by making two assumptions:
(cid:90) (cid:90)
i The bosonic system is represented by an ASWF; dSϕ˜ (R,S)= dSϕ (R,S) (30)
FSWF FSWF
ii Correlation factors (Jastrow) do not play a crucial
role and therefore can be omitted. Weseethatϕ˜ isafermionicshadowwavefunctionin
FSWF
which the shadow-shadow correlation has been replaced
Under these assumptions, (cid:104)w(cid:105)ρ associated to ΨFSWF by an effective local Jastrow J˜(R). Although needed
S
(cid:82) dRdS dS ϕ∗ (R,S )ϕ∗ (R,S ) for our proof, the reader should bear in mind that the
(cid:104)w(cid:105)FSWF = (cid:82) dRdS1dS2ϕ∗FSWF(R,S1)ϕF∗SWF(R,S2) normalizationfactorinvolvedcannotbeeasilyestimated,
1 2 ASWF 1 ASWF 2 its computation itself will lead to a sign problem.
(22)
We now integrate out the shadows, in order to elimi-
can be simplified using
nate the sign problem in the second integral of Eq. (28):
ϕ (R,S)(cid:39)e−C(R−S)2det(eikα·sβ)
FSWF (23) (cid:90)
ϕ (R,S)(cid:39)e−C(R−S)2det(eikα·rβ). Ψ¯FSWF(R)≡ dSϕ˜FSWF(R,S)
ASWF
wdehteerremeixnaacnttl,yaNndwwavheerveecwtoerassksuamreeoaccsuppinie-dpoinlatrhizeeSdlasytesr- =JR(R)J˜S(R)Φ(R) (cid:16)Cπ(cid:17)32N e−(cid:80)Ni=4C1k2i
tem with N(cid:90) fermions for simplification. We then have =JR(R)J˜S(R)Φ(R)e−(cid:80)Ni=4C1k2i (cid:90) dSΞ(R,S)
dSe−C(R−S)2det(eikα·sβ)= (cid:90)
(24) = dSϕ¯ (R,S)
(cid:16)π(cid:17)32N e−(cid:80)Ni=4C1k2i det(eikα·rβ) FSWF (31)
C
5
where we have formally reintroduced the shadows. We
canseethatΨ¯ canberegardedastheclosestASWF 1000
FSWF
to the given FSWF, i.e. a "bosonized" FSWF. Notice )]
²
that Ψ¯ (R) does not contain a Slater determinant K 100
FSWF c
evaluatedontheshadowcoordinatesanymore,hencethis e
s
wave function is not affected by the sign problem. 1/( 10
The FSWF can now be written as y [
c
(cid:90) n 1
e
Ψ (R)= dS [ϕ (R,S)−ϕ˜ (R,S) ci
FSWF FSWF FSWF fi
Ef 0.1
+ϕ¯ (R,S)]
FSWF
(cid:90)
=J (R) dSΞ(R,S)J˜(R)× 20 30 40 50
R S
N
(cid:20) (cid:18) (cid:19) (cid:21)
Φ(S) JS(S) −1 +Φ(R)e−(cid:80)Ni=4C1k2i
J˜(R) Figure 2. Efficiency of a VMC simulation of liquid 3He that
S
(32) employs the FSWF [38], fitted with an exponential function
∼e−kN, where k is a constant.
(cid:16) (cid:17)
If the term JS(S) −1 were zero, we would have solved
J˜S(R)
the sign problem. Inparticular,wewouldliketoknowhowthesignprob-
Let us now assume that we are able to sample from lem scales with the number of particles. For that we use
the signal, i.e. thefollowingreasoning: Supposethatwehaveperformed
acalculationwithN particleswhichprovidedusanesti-
ρ(R,S ,S )=ϕ¯∗ (R,S )ϕ¯ (R,S ), (33) 0
1 2 FSWF 1 FSWF 2 mate for our observable within a given error bar and the
corresponding efficiency η . We then change the number
then the weight corresponding to our original sampling, 0
of particles to N = κN . To obtain the same accuracy,
Eq. (32), writes 0
we must require that the noise term in the weight is the
(cid:20) (cid:18) (cid:19) (cid:21) same,i.e. (cid:104)noise(cid:105)=(cid:104)noise (cid:105). FromEq. (37),wecanread
w(R,S1,S2)= 1+ ΦΦ((SR1)) JJ˜S((SR1)) −1 e(cid:80)Ni=4C1k2i a dependence from the nu0mber of particles in the terms
×(cid:20)1+ ΦΦ((SR2)) (cid:18)SJJ˜S((SR2)) −1(cid:19) e(cid:80)Ni=4C1k2i(cid:21) (e3−7(cid:80))Ni=4mC1uk2ist(cid:39)beed−ecNCrρe2a/s3e.d Tbyheareffaocrteo,rteh−eκaNCv0eρr2a/g3.esAinssuEmq-.
S ingthevarianceofthenoiseintegrandtobeindependent
(34)
of N, we have σ ∝ (number of sampled points)−1/2 ∝
(cpu time)−1/2 and we can conclude that logη/η ∝−κ.
which can be seen as 0
Similarly, we can derive the dependencies for C and ρ,
weight=signal+noise (35) which gives
N
where signal=1. Hence, we require that logη/η ∝− ρ2/3. (38)
0 C
|(cid:104)noise(cid:105)|<ε|(cid:104)signal(cid:105)| (36) yieldingtheresultoftheprevioussubsectionundermore
general assumptions.
where ε<1 determines the final error of the calculation.
These prediction are confirmed by the numerical re-
We get
sults reported in [38] using FSWF trial wave functions
(cid:12)(cid:28) (cid:18) (cid:19)(cid:29) to compute the VMC energy of unpolarized liquid 3He
(cid:12)(cid:12) Φ(S1) JS(S1) −1 e(cid:80)Ni=4C1k2i (ρ = 0.016588 Å−3) in three dimensions. In Figures 2
(cid:12) Φ(R) J˜(R)
S and 3 these datas are fitted with an exponential func-
(cid:28) (cid:18) (cid:19)(cid:29)
+ Φ(S2) JS(S2) −1 e(cid:80)Ni=4C1k2i tion, demonstrating the exponential dependency on N
Φ(R) J˜(R) and C.
S (37)
(cid:28) (cid:18) (cid:19)
Φ(S )Φ(S ) J (S )
+ 1 2 S 1 −1
Φ2(R) J˜(R)
S IV. GENERALIZATION TO ANY KIND OF
×(cid:18)JS(S2) −1(cid:19)(cid:29) e2(cid:80)Ni=4C1k2i(cid:12)(cid:12)(cid:12)<ε ORBITALS
J˜(R) (cid:12)
S
In this section we are going to generalize the depen-
Equation(37)isanecessaryconditiontoavoidthesign dence of the efficiency on N, C, and ρ to the more gen-
problem, anditcontainsalltheinformationsthatweare eral case of Slater determinants which use any kind of
looking for. orbitals.
6
100000 Following the idea used in subsection IIIB, we write:
²)] ϕFSWF(R,S,K)=JR(R)Ξ(R,S)JS(S)Φ˜(K)Φpw(R,S)
K
ec 10000 ϕ˜FSWF(R,S,K)=JR(R)J˜S(R)Ξ(R,S)Φ˜(K)Φpw(R,S)
1/(s ϕ¯ (R,S,K)=J (R)J˜(R)Ξ(R,S)e−(cid:80)Ni=4C1k2i
cy [ FSWF R ×S Φ˜(K)Φ (R,K).
n 1000 pw
e (43)
ci
fi
Ef
Eq. (32) can be recasted as
100
(cid:90)
1.2 1.4 1.6 1.8 2 Ψ (R)∼J (R) dSdKΞ(R,S)J˜(R)×
1/C = 𝜏 [1/Ų] FSWF R S
(cid:20) (cid:18) (cid:19)
J (S)
Φ˜(K)Φ (S) S −1 (44)
Figure 3. Efficiency of a VMC simulation of liquid 3He that pw J˜(R)
S
employs the FSWF [38], fitted with an exponential function (cid:21)
∼e−k/C, where k is a constant. +Φ˜(K)Φ (R)e−(cid:80)Ni=4C1k2i ,
pw
while Eq. (37) writes
To accomplish this result it is sufficient to work in
Fteorumriienranstpaccaen. bTeheexpmreastsreidx aeslemanenintsteogfratlhoefSalaptreoddudcet- (cid:12)(cid:12)(cid:12)(cid:28)Φpw(S1,K1) (cid:18)JS(S1) −1(cid:19) e(cid:80)Ni=4C1k2i(cid:29)
of matrices over the N-particle momentum space: (cid:12) Φpw(R,K1) J˜S(R)
(cid:28) (cid:18) (cid:19) (cid:29)
+ Φpw(S2,K2) JS(S2) −1 e(cid:80)Ni=4C1k2i
φα(rβ)=(cid:90) dkαe−ikα·rβφ˜α(kα) (cid:28) Φpw(R,K2) J˜S(R) (cid:18) (cid:19) (45)
Φ (S ,K )Φ (S ,K ) J (S )
+ pw 1 1 pw 2 2 S 1 −1
∼(cid:90) dK γ(cid:88)N=1 (cid:0)e−ikγ·rβ(cid:1) (cid:16)φ˜α(kγ)δγα(cid:17) (39) Φpw(R,K1×)Φ(cid:18)pwJ(SR(S,K2)2−) 1(cid:19)J˜eS2(R(cid:80)Ni)=4C1k2i(cid:29)(cid:12)(cid:12)(cid:12)<ε.
(cid:90) J˜(R) (cid:12)
S
= dK(e−ikγ·rβ)·(Iφ˜α(kγ)), From Eq. (45) we see that the role played by e(cid:80)Ni=4C1k2i is
where φ˜(k) is the Fourier transform of the orbital φ(r), now played by
and (cid:28) (cid:29)
Φpw(S,K) e(cid:80)Ni=4C1k2i (46)
Φ (R,K)
I ≡φ˜ (k )δ . pw K
φ˜α(kγ) α γ γα
where (cid:104)...(cid:105) denotes the average over K obtained by
K
Therefore sampling from ϕ¯∗ (R,S ,K )×ϕ¯ (R,S ,K ) for
FSWF 1 1 FSWF 2 2
anygivenR,S andS . Inotherwords,thefactorwhich
1 2
(cid:90)
controlstheefficiencyofthecalculationisnowafunction
det(φα(rβ))∼ dK det(e−ikγ·rβ) det(Iφ˜α(kγ)) of R and S.
(cid:90) N (40) Since we have no simple interpretation of Eq. (46),
= dK det(e−ikγ·rβ) (cid:89)φ˜γ(kγ) wehavenumericallyestimateditsdependenceonthede-
greeoflocalizationoftheorbitalsemployedintheSlater
γ=1
determinant. For doing so we have neglected the Jas-
In the following we will use the notation trow terms, and considered only a kernel Ξ and a Slater
determinant Φ containing gaussian orbitals of the form
Φ(R)≡det(φ (r )) e−G(S−P)2,whereP labelssomelatticepositions. Figure
α β
4 shows our results.
Φ˜(K)≡ (cid:89)N φ˜ (k ) (41) The numerical results demonstrate that the efficiency
γ γ increases exponentially when the orbitals become more
γ=1
localized until it reaches a maximum and finally begins
todecrease. SuchdecreaseatlargeGisduetootherrea-
and
sons than the sign problem, as it affects also the ASWF.
Therefore, in order to isolate the sign problem depen-
Φ (R,K)≡det(e−ikα·rβ) (42) dency, we have introduced the ratio between the FSWF
pw
7
10
1000
Direct (pw)
1
Ry²)] 11000 Ry²)] 0.1 GDiDre (cpt w(1)s)
s s GD (1s)
1/( 1 1/( 0.01
cy [ 0.1 cy [ 0.001
n n
e 0.01 e
fici FSWF fici 0.0001
Ef 0.001 ASWF Ef0.00001
Ratio
0.0001
1x10-6
0.5 1 1.5 2 2.5 0.6 0.8 1 1.2 1.4 1.6 1.8 2
G [1/Bohr²] 1/C = 𝜏 [1/(Bohr²)]
Figure4. EfficiencyasafunctionofthegaussiancoefficientG Figure5. Comparisonbetweenadirectsimulationemploying
which enters in the Slater determinant orbitals. The dataset theFSWFandacalculationwhichemploystheG formula-
det
RatioreportsthedimensionlessratiobetweentheFSWFand tion. TheefficiencyiscalculatedfordifferentC andwithtwo
the ASWF values. Results were obtained by computing the different choices for the orbitals embedded in the Slater de-
expectation value of the energy per particle of the electronic terminant: The 1s orbitals and simple plane-waves (pw). To
structureofanhydrogenbccatomiccrystalwithr =1.7and obtain these results we simulated the electronic structure of
s
periodic boundary conditions. The lattice positions required 3D solid hydrogen with a bcc crystal structure and r =1.8,
s
by the gaussians are the hydrogen’s protons. We have used employingaYukawaJastrowbothforJ andJ . Theresults
R S
C =1. refer to simulations done with 16 atoms.
efficiency and the ASWF one, which are represented by to a direct employment of the FSWF with a value of C
triangular symbols in Fig. 4. Looking at these data, we as small as possible within given computational limits.
can notice that at large G there is a plateau rather than Such methods would greatly benefit from any method
a decay. able to alleviate the severity of the sign problem. In
The threshold at which the ratio reaches a plateau thepastyearssomeremarkableimprovementshavebeen
can be interpreted as the degree of localization at which achieved in this direction.
the fermionic statistics become irrelevant and the quan- We begin by outlining briefly the Gaussian Determi-
tum particles can be conveniently approximated as dis- nant (G ) method [35], introduced for the study of va-
det
tinguishable. canciesinsolid3He,andfurtherinvestigatedin[38]. The
In conclusion we have shown that there is a strong leadingideaistosumoverallthepossibleshadowpermu-
correlation between the localization of the Slater deter- tationsbymeansofananti-symmetricalkernelconsisting
minant orbitals and the sign problem which can be used of a determinant of gaussian orbitals. In practice, it is
to improve the efficiency of fermion simulations. sufficient to replace the kernel with
(cid:16) (cid:17)
Ξ(R,S) → Gdet(R,S)≡det e−C(rα−sβ)2 , (47)
V. ALLEVIATION OF THE FERMION SIGN
PROBLEM AND APPROXIMATED METHODS
where α and β label rows and columns of the matrix.
Figure 5 and 6 demonstrate that the exponential pre-
As we have shown in the previous Sections, Fermi factor is significantly reduced. The extension of the G
det
statistics entail a sign problem which in general makes techniquetoPIGSisstraightforward,asitissufficientto
simulationsoflargenumberoffermionsprohibitive. Nev- replace each gaussian kernel with a determinant of gaus-
ertheless, given the importance of fermionic systems for sians,exactlyasfortheSWF.InanunconstrainedDMC
the comprehension of many natural phenomena (e.g. method, where walkers diffuse according to the standard
quantum chemistry), different methods which cope with gaussianterm,wecanmultiplythebranchingprobability
these difficulties have been devised. Here we will focus by
on methods which involve Monte Carlo algorithms.
A first direct approach is to try to keep τ as small (cid:18) (R−R(cid:48))2(cid:19)
G (R,R(cid:48))/exp − , (48)
as possible, and employ a very good starting trial func- det 4τ
tion. In this way, it is possible to find a lower upper
bound to the exact ground state energy, systematically withtypicalvaluescloseto1forsmallτ. Anegativevalue
improvable by increasing the quality of the trial func- ofG impliesaswitchofthewalkerfromthepopulation
det
tion. This approach goes under the name of release-node carrying a positive sign to the one representing negative
ortransient estimates [23,24,42–45], anditcorresponds contributions.
8
1x106
Direct - C=2.0 (pw) 10
100000
Direct- C=2.0 (1s) AMD
10000 Direct - C=1.5 (pw) Ry²)] 1 Direct
1000 Direct - C=1.5 (1s) s
100 GD - C=2.0 (pw) 1/(
²)] 10 GD - C=2.0 (1s) y [ 0.1
Ry GD - C=1.5 (pw) nc
(s 1 GD - C=1.5 (1s) cie 0.01
1/ 0.1 fi
ncy [ 0.01 Ef0.001
e 0.001
ci
Effi 0.0001 0.7 0.8 0.9 1 1.1 1.2
1/C = 𝜏 [1/Bohr²]
0.00001
1x10-6
Figure7. Performancecomparisonbetweencalculationswith
1x10-7 theFSWFusingthedirectalgorithmandtheAMDmethod,
1x10-8 fordifferentvaluesofC. Theresultsrefertothecomputation
ofthepotentialenergyoftheelectronicstructureof16hydro-
1x10-9
20 40 60 80 100 120 genatomsinabcccrystalstructureatrs =1.8. Weemployed
N aYukawaJastrowbothforJR andJS,and1sorbitalswithin
the Slater determinant. For approximating the marginal dis-
tribution, we used a Slater determinant embedding orbitals
Figure 6. Same as in Fig. 5 but exploring the dependency of
from Quantum Espresso [46]. The AMD results refer to the
theefficiencyonN. Weshowtheresultsfordifferentorbitals
best efficiency attainable by varying the parameters Λ and
and choices of C. The simulated physical system was the
M (see [38]).
S
same as for Fig. 5.
10
A second method is the Approximated Marginal Dis- 1 AMD
)]
tribution Method (AMD) [38]. Here, the leading idea is Ry² 0.1 Direct
to sample the R coordinates from a modified sampling
s
distributionwhichisintendedtoaccountfortheintegra- 1/( 0.01
trieopnreosevnertatthioensohfatdhoewms,aargnidnatlhdeirsetfroirbeutpioronvfiodresRa. Abtetttheer ncy [ 0.001
e
sametime,theweightswhichcarrythesignarecomputed ci1x10-4
bsaymspulminmgionfgthuepsthhaedocwonst,raibccuotriodninsgctoomtihneggfrrooumpimngultteicphle- Effi1x10-5
nique. The improvements attainable with this method 1x10-6
can be seen in Fig. 7 and 8. We point out that this
20 40 60 80 100 120
method is valuable when the sign problem is "strong", N
whereasthedirectalgorithmoutperformsitintheoppo-
site limit. Figure 8. Same as in Fig. 7 but changing the number of
Several approximated techniques which avoid the sign simulated particles N. We used C =1.3 Bohr−2.
problem are used routinely. Probably the most common
oneisthefixed-nodeapproximation[42,43,47–49]which
imposes a nodal surface, forcing the solution to be anti- terminant of the FSWF) and to require that S and its
symmetrical. Concretely, this is accomplished by taking imaginary time projection R are in the same nodal re-
an antisymmetrical wave function ϕ and use the DMC gion.
T
scheme restricted to the positive (or negative) domain In the following we outline three general approaches
of ϕ : Whenever a walker crosses the nodal surface, it is which have been devised in the past and aimed to an ex-
T
suppressed. Throughthisprocedureitispossibletofilter act solution of the fermion sign problem. Even if they
the best antisymmetrical ground state within the given have not been able to fully overcome the exponential be-
nodal surface. Therefore, the nodal surface is the input haviour, they might have the potential to alleviate this
that will determine the quality of the final result. The trend.
fixed-nodeapproximationhasbeensuccessfullyextended In 1982, a Green’s function Monte Carlo algorithm for
to Path Integral Monte Carlo [50]. For transferring the fermions based on a cancellation process has been intro-
fixed-node method to FSWF calculations, it is sufficient duced[51]. Theproposedalgorithmisbasedontheintu-
to choose a nodal surface (typically using the Slater de- itive idea of having two different populations of walker,
9
one carrying a positive sign, and the other one a nega- accurate methods employed in quantum chemistry [58–
tive sign, which will cancel each other if they get "close 62].
enough". The method has been employed for few par-
ticle systems, but it suffers from a substantial drawback
that prevents its application to many-body system. The VI. COMPUTATIONAL COMPLEXITY OF THE
reason behind this is that the algorithm requires a high FERMION SIGN PROBLEM
density of walkers, and unfortunately such requirement
implies an exponential growth of the computational cost In this section we discuss the computational complex-
proportional to the number of simulated particles, be- ity of the fermion sign problem in QMC. Before going
cause of the increased dimensionality of the problem. into the details, we would like to briefly introduce the
Further works in the same direction [25, 26, 52, 53] have reader to the complexity classes P, NP, NPC, and NP-
shown that calculations for small system sizes are feasi- hard. The interested reader can refer to [14, 63, 64] for
ble,butthefermionsignproblemisnotsolvedingeneral. an exhaustive introduction to the topic.
In 1985 a method for treating the fermion sign prob- Pistheclassofproblemswhicharesolvableinpolyno-
lem involving mirror potentials has been devised [45], mial time: Provided an input problem of size n, it exists
where a fictitious repulsive interaction is used to keep an algorithm which can solve it in O(nk) time where k
the distribution of positive and negative walkers apart is a constant. In our specific case, the input is provided
from each other, and hence avoids the collapse into the bythephysicalparametersofthesystem,thevariational
same bosonic ground state. By increasing the repulsion parameters to be employed, and the Hamiltonian of the
between the two populations this method reduces to the system. The problem’s size is given by the number of
fixed-node approximation, whereas when it is set to zero simulated particles.
one obtains a release-node simulation. The mirror po- NP (Nondeterministic Polynomial) problems are the
tential method allows exact fermion calculations, but is ones which can be verified in polynomial time: Given a
limited to a small number of particles due to the expo- solution (certificate) of the problem, it exists an algo-
nential growth in number of walkers required to describe rithm that can verify its correctness in O(nk) time. We
themirrorpotential,similarlytowhathappenedwiththe remark that any problem in P belongs also to NP, be-
cancellation idea. To overcome this difficulty, it is pos- cause if it is possible to solve a problem in polynomial
sible to make use of a trial wave function. However this time, such a solution can be used to verify a certificate.
will lead to an approximate result, although potentially InordertoillustratetheNPCandtheNP-hardclasses
more accurate than the fixed-node one. we need to further introduce the concept of reducibility.
TheproblemAissaidtobepolynomial-time reducible to
InthecontextofPathIntegralMonteCarlo, therehas
problem B if it is possible to find a map A(cid:55)→B which is
been an interesting attempt towards the solution of the
computable in polynomial time.
fermionsignproblemin1998-2000[54,55]. Theproposed
The NPC (NP-Complete) problems are the NP prob-
approachgoesunderthenameofmultilevel blocking,and
lems which have the additional property of being
consists of distributing the integrals for the propagation
polynomial-timereducibletoanyotherNPproblem. Ifa
ofasingleimaginarytimestep∆τ inanelaboratedpyra-
problempossessesthelatterpropertybutnotnecessarily
midal structure, solving it in a bottom-up fashion. First
the first one, then it is said to be NP-hard.
one computes the N /2 integrals at the bottom, and
τ
IntheworkofTroyerandWiese[15],ithasbeenshown
then use these informations to compute the integrals at
that a general QMC algorithm which is able to compute
a coarser level, i.e. for 2∆τ. This procedures is repeated
the most general partition function can also be used to
untiltheintegralforthefullimaginary-timepropagation
solveaNPCproblem. ThereferenceNPCproblem[65]is
N ∆τ isfound. However,wehaveseenthatevenasingle
τ
the following: Given a classical 3D Ising spin glass with
integrationstepalreadyintroducesasignproblemwhich
Hamiltonian
scales exponentially in the number of particles, so that
the proposed scheme will eventually scale exponentially, (cid:88)
H =− J S S (49)
too. i1i2 i1 i2
<i1,i2>
In 2009 a new method, FCIQMC [56], for treating
fermion very accurately has been introduced, based on and a bound energy E , does a spin configuration with
0
a Monte Carlo imaginary-time projection technique per- energy E ≤ E exist? The interaction matrix J has
0
formed in the space of Slater determinants, in the spirit values j, 0, or −j chosen randomly, and the spins S can
of Full Configuration Interaction (FCI). In this method, have values ±1.
if the number of walkers is sufficient to populate such The latter decision problem is connected to Monte
space, the FCI wave function will emerge from the cal- Carlo calculations by the fact that provided a large
culation within a less severe computational cost com- enough inverse temperature β, the average energy of the
paredtothetraditionalone[57]. However,theuseofsuch spin glass system will be less than E + j/2 if a con-
0
variational space for treating particle correlation implies figuration with E ≤ E exists, and larger than E +j
0 0
a size-extensivity problem. Nevertheless, this method otherwise (basically, the simulated annealing minimiza-
has demonstrated to be competitive with other highly- tion method). As a consequence, the computation of E
10
with a Monte Carlo simulation would provide an answer method [70], whose efficiency will decrease exponentially
to the given 3D Ising spin glass NPC problem. In other with the given statistical error.
words, such Monte Carlo average is necessarily at least Noticethatthisdecisionalproblemdoesnotadmitthe
as hard as a NPC problem, i.e. it is NP-hard. existence of a certificate, as it is not in the form "does
Since the classical NPC problem in Eq. (49) can be ... such that ... exists?". Therefore one should find a
mapped in a quantum one simply by replacing classical corresponding problem which will allow to identify it as
spins with quantum ones, a general algorithm to solve aNP-hardproblem(ifitissuch). Unfortunately,wehave
quantumproblems(includingthosewithasignproblem) to leave this question open.
willprovidealsoasolutiontoourNPCproblemandthus
willbeNP-hard[15]. However,notallmany-fermionsys-
temsposeNPhardproblemsanditremainsopenifthere VII. CONCLUSION
is a criterion that allows us to immediately distinguish a
NP-hard situation from P or BPP. We have presented the SWF formalism and showed
It would be of great importance to be able to better that the fermion sign problem appearing in typical
characterise the fermion sign problem complexity. Is it imaginary-timeprojectioncontinuumQMCmethodscan
possibletodeviseacriterioninordertopredictwhenthe be reduced to the sign problem of the Fermionic Shadow
QMCsimulationwillbeNP-hard? Isitpossibletofinda Wave Function. This formalism was used to characterise
caseinwhichtheNP-hardnessoriginatefromthefermion bothanalyticallyandnumericallythefermionsignprob-
statistic? Some progress in this direction has been made lem,demonstratingitsdependenceonthenumberofpar-
recently [9, 27, 66–69] identifying sets of Hamiltonians ticles, length of the imaginary-time projection, and lo-
wihout sign problem. calisation of the system. Even though it seems that the
In the following we discuss a definition of the fermion exponential decay of the simulation efficiency cannot be
signproblemconcerningpurelycontinuumQMCsimula- overcome, we have shown that some methods can lead
tion with local interactions, using the SWF formalism. toasignificativereductionof itsexponentialfactor, thus
Decisional fermion sign problem: Is the value of the extending the applicability of exact QMC methods for
integral fermions. Concerning the complexity class of imaginary-
timeprojectionQMCalgorithms,aseparateproofofthe
(cid:90) NP-hardness of 2D and 3D fermionic systems with local
dSe−C(R−S)2det(φ (s )) J (S) (50)
α β S interactions in continuum space is still lacking.
strictly greater than zero?
If one could answer to this decisional problem, exact ACKNOWLEDGMENTS
fermion calculations in polynomial time would be possi-
ble by employing the fixed-node algorithm. If the pro- F.C.wouldliketoacknowledgetheNanosciencesFoun-
videdanswerisaffectedbyastatisticalerror,asitwould dationofGrenobleforfinancialsupportandT.D.Kühne
be the case by using a Monte Carlo technique, then the forallowingustoaccesstheMogonHPCwhichhasbeen
fixed-node method will have to make use of the penalty used for most of the numerical calculations.
[1] R.Blankenbecler,D.J.Scalapino,andR.L.Sugar,Phys. [8] P. Broecker and S. Trebst, Cond-Mat, arXiv:1511.02878
Rev. D 24, 2278 (1981), URL http://link.aps.org/ (2015).
doi/10.1103/PhysRevD.24.2278. [9] L.Wang,Y.-H.Liu,M.Iazzi,M.Troyer,andG.Harcos,
[2] J. E. Hirsch and R. M. R Martin Fye, Phys. Rev. Lett. Phys.Rev.Lett.115,250601(2015),URLhttp://link.
56, 2521 (1986). aps.org/doi/10.1103/PhysRevLett.115.250601.
[3] D. M. Ceperley, Journal of Statistical Physics 63, 1237 [10] M.A.Morales,R.Clay,C.Pierleoni,andD.M.Ceperley,
(1991), ISSN 1572-9613, URL http://dx.doi.org/10. Entropy 16, 287 (2014).
1007/BF01030009. [11] B. M. Austin, D. Y. Zubarev, and W. A. Lester, Chem.
[4] S. Chandrasekharan and U.-J. Wiese, Phys. Rev. Lett. Rev. 112, 263 (2012).
83, 3116 (1999), URL http://link.aps.org/doi/10. [12] J. Kolorenc and L. Mitas, Rep. Prog. Phys. 74, 026502
1103/PhysRevLett.83.3116. (2011).
[5] R.K.Kaul,Phys.Rev.B91,054413(2015),URLhttp: [13] R. J. Needs, M. D. Towler, N. D. Drummond, and P. L.
//link.aps.org/doi/10.1103/PhysRevB.91.054413. Rios, J. Phys.: Condens. Matter 22, 023201 (2010).
[6] M. G. Alford and A. Kryjevski, Journal of Physics G: [14] T.H.Cormen,C.E.Leiserson,R.L.Rivest,andC.Stein,
Nuclear and Particle Physics 37, 025002 (2010), URL Introduction to Algorithms (The MIT Press, 2009), 3rd
http://stacks.iop.org/0954-3899/37/i=2/a=025002. ed.
[7] Z.-X. Li, Y.-F. Jiang, and H. Yao, Phys. Rev. B [15] M. Troyer and U.-J. Wiese, Phys. Rev. Lett. 94,
91,241117(2015),URLhttp://link.aps.org/doi/10. 170201 (2005), URL http://link.aps.org/doi/10.
1103/PhysRevB.91.241117. 1103/PhysRevLett.94.170201.