SciPost Phys. Core 2, 011 (2020)

Fast and stable determinant quantum Monte Carlo
Carsten Bauer?
Institute for Theoretical Physics, University of Cologne, 50937 Cologne, Germany? bauer@thp.uni-koeln.de

Abstract
We assess numerical stabilization methods employed in fermion many-body quantum
Monte Carlo simulations. In particular, we empirically compare various matrix decomposition and inversion schemes to gain control over numerical instabilities arising in
the computation of equal-time and time-displaced Green’s functions within the determinant quantum Monte Carlo (DQMC) framework. Based on this comparison, we identify
a procedure based on pivoted QR decompositions which is both efficient and accurate
to machine precision. The Julia programming language is used for the assessment andimplementations of all discussed algorithms are provided in the open-source software
library StableDQMC.jl.
Copyright C. Bauer.
This work is licensed under the Creative Commons
Attribution 4.0 International License.
Published by the SciPost Foundation.

Received 14-03-2020
Accepted 26-06-2020
Check forupdatesPublished 29-06-2020
doi:10.21468/SciPostPhysCore.2.2.011

Contents
1 Introduction

1

2 Determinant Quantum Monte Carlo2

3 Numerical instabilities

3

4 Stabilization: time slice matrix multiplications
4.1 Stabilization scheme
4.2 Matrix decompositions
4.2.1 SVD (U DV † )4.2.2 QR (U DT )4.3 Benchmarks
4.3.1 Accuracy
4.3.2 Efficiency

4
4
6
6
6
6
6
7

5 Stabilization: equal-time Green’s function
5.1 Inversion schemes
5.2 Benchmarks
5.2.1 Accuracy
5.2.2 Efficiency

7
8
9
9
10

6 Stabilization: time-displaced Green’s function

101

SciPost Phys. Core 2, 011 (2020)

6.1
6.2

Inversion schemes
Benchmarks
6.2.1 Accuracy
6.2.2 Efficiency

11
12
12
13

7 Discussion

14

A Inversion schemes for time slice matrix stacks
A.1 QR/UDT
A.2 SVD

1415
15

B Spin-fermion model of a metallic quantum critical point

15

References

161

Introduction

Many-fermion systems play an important role in condensed matter physics. Due to their intrinsic correlations they feature rich phase diagrams which can not be captured by purely classicalnor non-interacting theories. Especially at the lowest temperatures, quantum mechanical fluctuations driven by Heisenberg’s uncertainty principle become relevant and lead to novel phasesof matter like superconductivity and states beyond the Fermi liquid paradigm [1,2]. Because of
the presence of interactions, predicting microscopic and thermodynamic properties of fermion
many-body systems is inherently difficult. Analytical approaches are typically doomed to failin cases where one can not rely on the smallness of an expansion parameter [3].
Fortunately, the determinant quantum Monte Carlo (DQMC) method [4–8] overcomes thislimitation. The key feature of DQMC is that it is numerically exact - given sufficient computation time the systematical error is arbitrarily small. Provided the absence of the famous signproblem [9,10], it allows for an efficient exploration of the relevant region of the exponentiallylarge configuration space in polynomial time. It is an important unbiased technique for obtaining reliable insights into the physics of many-fermion systems which, among others, hasbeen applied to the attractive and repulsive Hubbard model [11–13], the Kane-Mele-Hubbard[14], and metallic quantum criticality, including studies of antiferromagnetic [1–3, 15], Isingnematic [16], and deconfined quantum critical points [17] where fermionic matter fields arecoupled to bosonic order parameters.
Although conceptually straightforward, care has to be taken in the implementation of
DQMC because of inherent numerical instabilities arising from ill-conditioned matrix exponentials. Over time, stabilization schemes [5,8,18–20] based on various matrix factorizations,such as singular value decomposition (SVD), modified Gram-Schmidt, and QR decomposition,
have been proposed for lifting these numerical issues. It is the purpose of this manuscript
to review a subset of these techniques and to compare them with respect to accuracy andspeed. Particular emphasis is placed on concreteness and reproducibility: we provide implementations of all discussed algorithms as well as the code to recreate all visualizations in thismanuscript in form of the software library StableDQMC.jl. We choose the open-source, highlevel programming language Julia [21,22] for our assessment which has proven [23,24] to becapable of reaching a performance comparable to established low-level languages in the fieldof numerical computing. Readers are invited to open issues and pull requests at the library

2

SciPost Phys. Core 2, 011 (2020)

repository to discuss, improve, and extend the list of stabilization routines. Beyond reproducibility, the software library will also serve as an important abstraction layer allowing users
to focus on physical simulation instead of numerical implementation details.
Specifically, the structure of the manuscript is as follows. We start by providing a brief
introduction into the DQMC method in Sec. 2. In Sec. 3 we illustrate numerical instabilities
in the DQMC and discuss their origin. Following this, we demonstrate (Sec. 4) how matrix
factorizations can be utilized to remedy these numerical artifacts in chains of matrix products.
In Sec. 5 we present and benchmark different schemes for stabilizing the computation of theequal-times Green’s function, the fundamental building block in DQMC. Lastly, we turn to thecalculation of time-displaced Green’s functions in Sec. 6 before concluding and summarizingin Sec. 7.

2

Determinant Quantum Monte Carlo

We begin by reviewing the essentials of the determinant quantum Monte Carlo method [4–8].
We assume a generic quantum field theory that can be split into a purely bosonic part SB
and a contribution S F from itinerant fermions. The latter comprises both fermion kinetics, T ,and boson-fermion interactions, V . A famous example is given by the Hubbard model after a
decoupling of the on-site interaction U ni,↑ ni,↓ by means of a continuous Hubbard-Stratonovich
or a discrete Hirsch transformation [25]1 . The quantum statistical partition function is givenbyZ
†Z = D ψ, ψ† , φ e−SB [φ]−SF [ψ,ψ ,φ] .
(1)
The first step in DQMC is to apply the quantum-classical mapping [26] and switch from thed dimensional quantum theory above to a D = d + 1 dimensional classical theory. Here, theextra finite dimension of the classical theory is given by imaginary time τ and has an extent
proportional to inverse temperature β = 1/T . Discretizing imaginary time into M slices,
β = M ∆τ, and applying a Trotter-Suzuki decomposition [27, 28] one obtains



Z
M
X


†
−SB
Z = Dφ e Tr exp −∆τψ T + Vφ ψ .
(2)
l=1

A separation of the matrix exponential then leads to a systematic error of the order O ∆τ2in the partition function,



eA+B ≈ eA e B ,∆τ∆τe−∆τ(T +V ) ≈ e− 2 T e−∆τV e− 2 T + O ∆τ3 ,
 m 
Z
Y

Z=

∆τ

†

†

∆τ

Dφ e−SB Tr

l=1


Bl + O ∆τ2 .

(3)†Here, Bl = e− 2 ψ T ψ e−∆τψ Vφ ψ e− 2 ψ T ψ are imaginary time slice propagators. Note that†the contribution e−∆τψ Vφ ψ depends on the bosonic field φ due to a potential fermion-boson
coupling in V . Rewriting the trace in (3) as a determinant [8] yields the fundamental formZZ = Dφ e−SB det Gφ−1 + O ∆τ2 ,(4)
Depending on the decomposition channel, the bosonic field φ represents either spin or charge fluctuations inthis case.
1

3

SciPost Phys. Core 2, 011 (2020)

where
G = [1 + B M B M −1 · · · B1 ]−1(5)

is the equal-time Green’s function [26]. Accordingly, the Metropolis probability weight is given
by
§
ª
det G
p = min 1, e−∆Sφ.
(6)
det G 0This implies that, considering a generic, global update one needs to compute the Green’s function G and it’s determinant in each DQMC step2 .
Importantly, it is only under specific circumstances, such as the presence of a symmetry,that the integral kernel of the partition function can be safely interpreted as a probability
weight since Gφ and its determinant are generally complex valued. This is the famous sign
problem [29].

3

Numerical instabilities

To showcase the typical numerical instabilities arising in the DQMC framework we consider
the Hubbard model in one dimension at half filling,


X
X
1
1
ni↑ −
ni↓ −
,
(7)
H = −tci† c j + U
22i〈i, j〉

which is free of the sign-problem [29]. We set the hopping amplitude to unity, t = 13 .
As seen from Eq. (5), the building block of the equal-time Green’s function is a matrix
chain multiplication of imaginary time slice matrices. To simplify our purely numerical analysis
below we assume that these slice matrices Bi are independent of imaginary time,B(β, 0) ≡ B M B M −1 · · · B1 = |BB{z· · · B},(8)

M factors

which, physically, amounts to assuming a constant bosonic field φ = const.
First, we consider the non-interacting system, U = 0. As apparent from Fig. 1, a naive
computation of Eq. 8 fails for β ≥ βc ≈ 10. Leaving a discussion of the stabilization of thecomputation for the next section, let us highlight the origin of this instability. The eigenvaluesof the non-interacting system are readily given by
εk = −2t cos(k),(9)

such that energy values are bounded by −2t ≤ εk ≤ 2t. A single positive definite slice matrix
B = e−∆τT therefore has a condition number of the order of κ ≈ e4|t|∆τ and, consequently,
B(τ, 0) has κ ≈ e4|t|M ∆τ = e4|t|β . This implies that the scales present in B(τ, 0) broaden
exponentially at low temperatures T = 1/β leading to inevitable roundoff errors due to finite
machine precision which spoil the result.
We can estimate the expected inverse temperature of this breakdown for the data type
Float64, that is double floating-point precision according to the IEEE 754 standard [30], bysolving κ(β) ∼ 10−17 for βc . One finds βc ≈ 10 in good agreement with what is observed inFig. 1a. Switching to the data type Float1284 (quadruple precision) with βc ≈ 20 in Fig. 1b,
the onset of roundoff errors is shifted to lower temperatures in accordance with expectations.
2

For local updates one can typically avoid those explicit calculations and compute the ratio of determinants in
Eq. (6) directly [2].
3
We will consider the canonical discrete decoupling [25] in the spin channel due to Hirsch in our analysis.
4
The datatype Float128 is provided by the Julia package Quadmath.jl.

4

60

log singular values

log singular values

SciPost Phys. Core 2, 011 (2020)

30
0
-30

stable
naive

-60
0

10

20

30

stable
naive

60
30
0
-30
-60

40

0

inverse temperature β

10

20

30

40

inverse temperature β(a) 64-bit floating point (Float64)(b) 128-bit floating point (Float128)

Figure (1) Numerical instabilities due to finite machine precision arising in thecalculation of the time slice matrix chain product B M B M −1 · · · B1 for model (7). Different lines represent logarithmic singular values as observed in naive (green) andarbitrary precision computations (orange) for a N = 4 system. Due to (quasi) degeneracies only 5 out of 8 singular values are visually distinguishable. The dashed line
(grey) indicates the expected floating point accuracy5 .

4Stabilization: time slice matrix multiplications

4.1

Stabilization scheme

A trivial solution to the issue outlined above is to perform all numerical operations with arbitrary precision. In Julia, this can be realized by means of the BigFloat data type6 . However,this comes at the expense of (unpractical) slow performance due to algorithmic overhead andlack of hardware support. Arbitrary precision numerics is nevertheless a valuable tool and wewill use it to benchmark the accuracy of stabilization methods below7 .
How can we get a handle on the numerical instabilities in a floating point precision computation? As has been realized [18] soon after the introduction of the DQMC method [4],an effective strategy is to keep the broadly different scales in the matrix exponentials separated throughout the computation (as much as possible) and only mix them in a final step, if
necessary. A useful tool for extracting the scale information is a matrix decomposition,
B = U DX .

(10)

Here, U and X are matrices of order unity and D is a real diagonal matrix hosting the exponentially spread scales of B. We will refer to the values in D as singular values independentof the particular decomposition. Using Eq. (10), we can stabilize the matrix multiplication ofslice matrices B1 and B2 in Eq. (8) as follows [18], (fact_mult in StableDQMC.jl)
B2 B1 = U2 D2 X 2 U1 D1 X 1 = U2 (D2 ((X 2 U1 )D1 )) X 1 ) = U r Dr X r .
| {z } | {z }
|
{z
}
B2

(11)

U 0 D0 X 0

B1

We estimate the precision as p = log10 (2fraction ), where fraction is the mantissa of a given binary floating point
format. This gives p ∼ 16 for Float64 and p ∼ 34 for Float128.
6
Technically, BigFloat has a finite, arbitrarily high precision.
7
For our non-interacting model system one can alternatively simply diagonalize the Hamiltonian and calculate
the Green’s function exactly.
55

SciPost Phys. Core 2, 011 (2020)

Here, U r = U2 U 0 , Dr = D0 , X r = X 0 X 1 , and U 0 D0 X 0 indicates an intermediate matrix decomposition. If we follow this scheme, in which parentheses indicate the order of operations,largely different scales present in the diagonal matrices won’t be additively mixed throughout
the computation. Specifically, note that the multiplication of the well-conditioned, combined,
unit-scale matrix U = X 2 U1 with D1 and D2 does preserve the scale information: the diagonal
matrices merely rescale the columns and rows of U,


s



s
 s

s
s
s

s


D2 U D1 = 

s

|

s
s
s
s

s
s
s
s

{z

unit scale



s
s 

s
s

s

s



s
 s

s
s
s

s




s

(12)

s

}

s ss s2 s s ss s2 s s=ss ss s2 s
s ss s2 s
 2
s s sss ss2 ss
 sss s2 s ss2 ss =
 ss2 ss2 s3 s2  .
ss ss s2 s 2




s

s

(13)

s

s

s

s

(14)

s

s

s

s

s

Repeating the procedure (11), we obtain a numerically accurate U DX decomposition of
the full time slice matrix chain B(τ, 0), which preserves the scale information as indicated inFig. 1.8 We note in passing that in practice it is often unnecessary to stabilize every individual
matrix product. Instead one typically performs a mixture of naive and stabilized products forthe sake of speed while still retaining numerical accuracy [8].

4.2Matrix decompositions

There are a various matrix decompositions that one could employ to obtain the factorization
B = U DX , Eq (10). In the following we will consider the two most popular choices in DQMC
codes [7, 8, 18].
4.2.1

SVD (U DV † )The singular value decomposition (SVD) is given by
B = USV † ,(15)where U and V † are unitary and S is real and diagonal.
For computing the SVD of a matrix of regular floating point precision (Matrix{Float64}),
Julia utilizes the heavily optimized routines provided by LAPACK9 [31]. Concretely, there exist
three different implementations of SVD algorithms [32]:10• gesdd (default): Divide-and-conquer (D&C)
8
Note that we do not discuss the faster way to calculate B M as U D M X . This is intentional since most real systems
will involve fermion-boson interactions and the slice matrices will depend on φ(τ).
9We will report on results obtained with the LAPACK implementation OpenBLAS that ships with Julia. Qualitatively similar results have been found in an independent test based on Intel’s Math Kernel Library (MKL).
10
Note that the names of LAPACK functions typically encode properties of the input matrix such as realness or
symmetry. In Julia multiple-dispatch takes care of routing different matrix types to different methods. The Julia
function gesdd works for both real and complex matrices, that is there is no (need for) cgesdd.

6

SciPost Phys. Core 2, 011 (2020)

100

log singular values

50
0
-50

exact
QR
SVD (Jacobi)
SVD
SVD (D&C)

-100
-150
0

10

20

30

40

50

inverse temperature β
Figure (2) Comparison of matrix decompositions to heal the numerical instabilities in the calculation of the time slice matrix chain product B M B M −1 · · · B1 formodel (7). Different lines represent logarithmic singular values as observed in stabilized computations. The QR (orange) and Jacobi SVD singular values (yellow)
lie on top of the exact result (red) whereas both the regular SVD (green) and thedivide-and-conquer SVD (purple) show large deviations at low temperatures β ¦ 25(∆τ = 0.1).
• gesvd: Bidiagonal QR iteration (conventional)• gesvj: Jacobi algorithm (through JacobiSVD.jl)
To simplify the manual access to these algorithms we export convenience wrappers of thesame name in StableDQMC.jl. We will compare all three variants below and benchmark
them against an arbitrary precision computation using BigFloat. Since LAPACK doesn’t support special number types, we will utilize the native-Julia SVD implementation provided byGenericSVD.jl in this case.
4.2.2

QR (U DT )A QR decomposition reads
B = QR = U DT,(16)

where Q is unitary, R is upper triangular, and we have split R into a diagonal part D and an
upper triangular part T in the second step. Specifically, U = Q is unitary, D = diag(R) is realand diagonal, and T is upper triangular.
In Julia, one can obtain the QR factored form of a matrix using qr from the standard libraryLinearAlgebra. We will consider the pivoted QR, which is deployed in the public DQMC
implementations ALF [33] and QUEST [34], in form of LAPACK’s geqp3 in our analysis. A
factorization into U DT form is provided by functions udt and udt! in StableDQMC.jl.

7

slowdown factor rel. to QR

SciPost Phys. Core 2, 011 (2020)

25

SVD
SVD (D&C)

20

SVD (Jacobi)
QR

15
10
5
0
100

200

300

400

500

N × N square matrix
Figure (3) Computational efficiency of matrix decompositions. Shown is theruntime cost of the factorization of a complex matrix of size N × N by means ofvarious SVD algorithms relative to the QR decomposition.

4.3

Benchmarks

4.3.1

Accuracy

Supplementing our general considerations above, we test the correctness of the matrix productstabilization procedure with respect to varying the concrete SVD and QR factorization algorithms. Fig. 2 shows the logarithmic singular values of the time slice matrix chain B(β, 0) as afunction of inverse temperature β obtained from employing different matrix decompositions.
Clearly, the accuracy of the computed singular values shows a strong dependence on the chosen factorization algorithm. While the results for the QR decomposition and Jacobi SVD seemto fall on top of the exact result, we observe large deviations for the conventional and D&C
SVD algorithms. This effect is particularly pronounced at low temperatues, β ¦ 25. The factthat small scales are lost in these SVD variants, while large ones are still correct, can be understood from LAPACK’s SVD error bounds [35]: The error is bounded relative to the largest
singular value. Thus, large scales are computed to high relative accuracy and small ones may
not be.
4.3.2

Efficiency

Turning to computational efficiency, we illustrate runtime cost measurements for all considered
SVD variants relative to the QR decomposition in Fig. 3. We find that both the conventional
SVD and Jacobi SVD are an order of magnitude slower than the QR decomposition while
only the divide-and-conquer algorithm shows comparable speed. Among the SVD variants,
the Jacobi SVD is the most costly by a large margin, having about twice the runtime of theconventional SVD for small system sizes.

5

Stabilization: equal-time Green’s function

Similar to the considerations above, a naive computation of the Green’s function according
to Eq. (5) is potentially unstable because of numerical roundoff errors due to finite machine
8

SciPost Phys. Core 2, 011 (2020)

precision. In particular, adding the identity to the ill-conditioned slice matrix chain B(τ, 0)will generally wash out small singular values and will lead to a non-invertible result such thatthe subsequent inversion in Eq. (5) is ill-defined. This clearly prohibits a safe calculation of
the equal-time Green’s function and asks for numerical stabilization techniques.

5.1

Inversion schemes

As for the time slice matrix products in Eq. (11), the strategy will be to keep exponentiallyspread scales as separated as possible. A straightforward scheme [7, 8] (inv_one_plus) to
add the unit matrix and perform the inversion of 1 + B(τ, 0) in a stabilized manner is given by
G = [1 + U DX ]−1 = [U (U † X −1 + D) X ]−1 = [(Uu)d(x X )]−1 = U r Dr X r ,
|
{z
}

(17)

ud x

where U r = (x X )−1 , Dr = d −1 , and X r = (Uu)−1 . Here, the intermediate addition (parentheses
in the second line of (17)) of unit scales and singular values is separated from the unitaryrotations such that U † X −1 only acts as a clean cutoff,ss
U † X −1 + D = 
s
s

s
s
s
s

s
s
s
s





s
s  
+
s
s

s



s



s s s s
s

  s
=

s

s
s

s

s
s

s
s
s

s 
.
s
s

(18)

As we will demonstrate for the time-displaced Green’s function in Sec. 6, a procedure like
Eq. (17) based on a single intermediate decomposition will still fail to give accurate results forsome of the matrix decompositions. For this reason, we consider another stabilization procedure put forward by Loh et al. [5, 18] (inv_one_plus_loh), in which one initially separates
the scales of the diagonal matrix D into two factors Dp = max(D, 1) and Dm = min(D, 1),
Dp = 

s

s
s






,


Dm = 



s

s

s


,

s

(19)

s

and performs two intermediate decompositions,
G = [1 + U DX ]−1 = [1 + U Dm Dp X ]−1 = [(X −1 Dp−1 + U Dm )Dp X ]−1= X −1 [Dp−1 (X −1 Dp−1 + U Dm )−1 ] = U r Dr X r ,
{z
}
|
ud x
|
{z
}

(20)
(21)

ud x

where U r = X −1 u, Dr = d, and X r = x.

5.2

Benchmarks

We assess how different matrix decomposition algorithms perform in stabilized computationsof B(β, 0), the Green’s function G both with respect to accuracy and speed. All results arefor the Hubbard model, Eq. 7, with the interaction strength set to U = 0 and U = 1 (alphatransparent in all plots)

9

SciPost Phys. Core 2, 011 (2020)

-5

-5

∆G

0

∆G

0

-10

-10

-15

-15
10

20

30

40

10

inverse temperature β

20

30

40

inverse temperature β(a) inv_one_plus, Eq. (17)(b) inv_one_plus_loh, Eq. 20

Figure (4) Accuracy of the Green’s function obtained from stabilized computations using the listed matrix decompositions and inversion schemes. Shown is∆G = log(max(abs(G − Gexact ))) for U = 0 (solid) and U = 1 (alpha transparent).
5.2.1

Accuracy

Starting from a stabilized computation of B(β, 0), Sec. 4, we calculate the equal-time Green’s
function by performing the inversion according to the schemes outlined above and varying
the applied matrix factorization. In Fig. 4a we show our findings for inv_one_plus, Eq. 17,
where we have taken the maximum absolute difference between the computed and the exact
Green’s function as an accuracy measure. At high temperatures and for U = 0, we observe thatall decompositions lead to a good approximation of Gexact with an accuracy close to floating
point precision. However, when turning to lower temperatures the situations changes dramatically. We find that only the QR decomposition and the Jacobi SVD deliver the Green’s
function reliably. Compared to the other SVD variants, which fall behind by a large marginand fail to reproduce the exact result, they consistently show about optimal accuracy evenin the presence of interactions. As displayed in Fig. 4b, switching to the inversion scheme
inv_one_plus_loh, Eq. 20, does generally improve the accuracy but deviations of the regular SVD and D&C SVD remain of the order of unity at the lowest temperatures.
These findings suggest that only the QR decomposition and the Jacobi SVD, irrespective of
the inversion procedure, are suited for computing the equal time Green’s function in DQMC
reliably.
5.2.2

Efficiency

Independent of the employed inversion scheme, matrix decompositions are expected to be theperformance bottleneck in the Green’s function computation. We hence expect the speed differences apparent in Fig. 3 to dominate benchmarks of the full Green’s function calculation as
well. This anticipation is qualitatively confirmed in Fig. 5, which shows the runtime cost of
the Green’s function computation for both inversion schemes and all matrix decompositions
relative to the QR. While the divide-and-conquer SVD is in the same ballpark as the QR decomposition the other SVD algorithms fall behind by a large margin (an order of magnitude)
for both inversion procedures. Importantly, this apparent runtime difference is increasing
with system size. The observation that the relative slowdown factor is larger for the inversion
scheme inv_one_plus_loh can be understood from the fact that it requires one additional
intermediate matrix decomposition.
10

slowdown factor rel. to QR

slowdown factor rel. to QR

SciPost Phys. Core 2, 011 (2020)

10.0
7.5
5.0
2.5
100

200

300

400

500

40
30
20
10
0
100

number of sites N

200

300

400

500

number of sites N(a) inv_one_plus, Eq. (17).

(b) inv_one_plus_loh, Eq. (20).

Figure (5) Efficiency of the stabilized Green’s function calculation using thelisted matrix decompositions and inversion schemes. Shown are results for U = 0.
Combined with the accuracy results these findings suggest that among the QR decomposition and the Jacobi SVD, which are found to be reliable in both inversion schemes, the QR
decomposition has a significantly lower runtime cost and is therefore to be preferred in DQMC.

6

Stabilization: time-displaced Green’s function

In this section, we turn to the stabilization of time-displaced Green’s functions. While these
are not required in the basic DQMC, that is for generating a representative Markov chain ofconfigurations, they are central to the measurement of time-displaced correlation functions
such as pairing correlations and the superfluid density [6, 7].
First, we generalize our definition of the Green’s function, Eq. 5, to include imaginary timeτ = l∆τ,G(τ) = 〈ci c †j 〉φl = [1 + Bl−1 . . . B1 B M . . . Bl ]−1 .

(22)Note that G ≡ G1 = G M +1 = [1 + B M . . . Bl ]−1 due to fermionic boundary conditions. The
time displaced Green’s function can now be defined in terms of the time ordering operator T
as [7, 8]
Gl1 ,l2 ≡ G(τ1 , τ2 ) ≡ 〈Tci (τ1 )c †j (τ2 )〉ϕ .
More explicitly, this reads¨
G(τ1 , τ2 ) =Bl1 · · · Bl2 +1 Gl2 +1 ,
τ1 > τ2 ,
−1− 1 − Gl1 +1 Bl2 · · · Bl1 +1, τ2 > τ1 .

(23)

In principle, this gives us a prescription for how to calculate G(τ1 , τ2 ) from the equal time
Green’s function discussed in Sec. 5. However, when |τ1 − τ2 | is large a naive calculationof the matrix products in Eq. 23 would be numerically unstable, as seen in Sec. 4. More
importantly, by first calculating the equal-time Green’s function one already mixes (and looses)
scale information in the last recombination step, G = U DX . It is therefore advantageous to11

SciPost Phys. Core 2, 011 (2020)

compute the time-displaced Green’s function directly as (assuming τ1 > τ2 for simplicity)G(τ1 , τ2 ) = Bl1 · · · Bl2 +1 Gl2 +1
−1= Bl1 · · · Bl2 +1 1 + Bl2 . . . B1 B M . . . Bl2+1

−1


= Bl−1+1 · · · Bl−1 + Bl2 . . . B1 B M . . . Bl1+1 
2
1
|
{z
}
|
{z
}

(25)

(26)

UR DR X R

U L DL X L= [U L DL X L + UR DR X R ]

6.1

(24)

−1

.

(27)

Inversion schemes

As for the equal time Green’s function (Sec. 5), one must be careful to keep the scales in DL
and DR separated when performing the summation and inversion to avoid unnecessary floating
point roundoff errors. As a first explicit procedure, we consider a simple generalization of
Eq. 17 (inv_sum),G(τ1 , τ2 ) = [U L DL X L + UR DR X R ]−1= [U L (DL X L X R−1 + U L† UR DR ) X R ]−1
|
{z
}
ud x

= [(U L u)d−1(x X R )]−1(28)

= U r Dr X r ,where U r = (x X R )−1 , Dr = d −1 , and X r = (U L u)−1 .
Analogously, we can generalize the scheme by Loh et al. [5], Eq. 20, in which we split thescales into matrix factors Dm = min(D, 1), Dp = max(D, 1), (inv_sum_loh)
G(τ1 , τ2 ) = [U L DL X L + UR DR X R ]−1= [U L DLm DLp X L + UR DRm DRp X R ]−1

−1




DLm
DRm
†
−1

=  U L DL p
X L X R + U L UR
X R DRp 

DRp
DLp
|
{z
}
ud x

= X R−1

11[ud x]−1
U†
DRp
DLp L
|
{z
}

(29)

ud x

= U r Dr X r ,where U r = X R−1 u, Dr = d, and X r = x U L† .
We note that Hirsch [7, 36] has proposed an alternative method for computing the timedisplaced Green’s function based on a space-time matrix formulation of the problem. Although
this technique has been successfully deployed in many-fermion simulations we won’t discuss
it here because of its subpar computational scaling: for a system composed of N lattice sites,
fermion flavors f , and imaginary time extent M one has to invert (naively a O(x 3 ) operation)a matrix which takes up O((N M f )2 ) memory. Similarly, Assaad et al. [8] have described
an approach to compute both equal time and time-displaced Green’s functions in one step.
However, this requires to work with extended matrices of doubled linear dimension comparedto the regular Green’s functions.
12

SciPost Phys. Core 2, 011 (2020)

5

0

0

∆G

∆G

-5
-5

-10

-10

-15

-15
0

10

20

30

0

40

10

20

30

40

imaginary time τ

imaginary time τ(a) inv_sum, Eq. (28).

(b) inv_sum_loh, Eq. (29).

Figure (6) Accuracy of the time-displaced Green’s function obtained from stabilized computations using the listed matrix decompositions and inversion schemes.
Shown is ∆G = log(max(abs(G(τ, 0) − Gexact (τ, 0)))) for β = 40.

6.2

Benchmarks

6.2.1

Accuracy

In Fig. 6, we show the logarithmic, maximal, absolute deviation of the time-displaced Green’s
function from the arbitrary precision result as a function of time-displacement τ at inverse
temperature β = 40. Focusing on the inversion scheme inv_sum first, Fig. 6a, both regularand D&C SVD clearly fail to capture the intrinsic scales sufficiently and errors much beyond
floating point precision are visible. Compared to these SVD variants, the QR decomposition
systematically leads to equally or more accurate results. However, it clearly fails to be reliable
at long times τ ∼ β/2 (the Green’s function is anti-periodic in τ). Only the Jacobi-method
based SVD produces accurate Green’s function values for all considered imaginary times.
Switching to the inversion scheme inv_sum_loh, the situation changes, as illustrated inFig. 6b. While the non-Jacobi SVDs still have insufficient accuracy, the result for the QR decomposition improves dramatically compared to inv_sum and leads to stable Green’s function
estimates up to floating point precision along the entire imaginary time axis.
Similar to our findings for the equal-time Green’s function, this suggests that only the
Jacobi SVD and the QR decomposition are reliable for DQMC. The latter, however, must be
paired with the inv_sum_loh inversion scheme to be reliable. To the best of our knowledge,
this finding has not yet been mentioned in the literature.
6.2.2

Efficiency

Finally, we compare the computational runtime cost associated with both stable approaches:
the Jacobi SVD combined with the regular inversion and the QR decomposition paired with
inv_sum_loh. As shown in Fig. 7, we find that the latter is consistently faster for all considered system sizes. In relative terms, the SVD based approach falls behind by at least afactor of two and seems to display inferior scaling with the chain length N . This indicates that
QR decompositions should be preferred over singular value decompositions when computing
time-displaced Green’s functions, in spite of the need to use an inversion scheme of higher
complexity.

13

SciPost Phys. Core 2, 011 (2020)

60
rel. to QR

computation time

QR
SVD (Jacobi)

40

10
8
6
4
2
0

100

200

300

20

0
0

100

200

300

number of sites N
Figure (7) Efficiency of the time-displaced Green’s function obtained from stabilized computations using the QR decomposition in combination with the inversion
scheme inv_sum_loh, Eq. (29) and the Jacobi SVD paired up with the regular inversion scheme inv_sum, Eq. (28). Measurements are taken over multiple runs at
τ = β/2 = 20. The inset show the slowdown of the Jacobi SVD relative to the QR
based approach.

7

Discussion

Numerical instabilities arise naturally in finite machine-precision quantum Monte Carlo simulations of many-fermion systems. Different schemes based on matrix factorizations have been
proposed to handle the intrinsic exponential scales underlying these instabilities in a stablemanner. As we have shown in this manuscript, these techniques can have vastly different
accuracy and efficiency rendering them more or less suited for determinant quantum Monte
Carlo simulations.
For our test system, the one-dimensional Hubbard model, we find that conventional anddivide-and-conquer based singular value decompositions consistently fail to produce accurate
equal time and time-displaced Green’s functions, in particular at the lowest considered temperatures, β ∼ 40. Only the QR decomposition and the Jacobi-method based SVD are able tostabilize the computation and produce reliable results. Importantly, we observe that in caseof the time-displaced Green’s function, the QR must be paired with an inversion scheme put
forward by Loh et al. [18], an observation that, to the best of our knowledge, has not been
mentioned in the literature before. No such qualitative dependence on the inversion procedure
is observed for the Jacobi SVD.
In terms of efficiency, we find that the QR decomposition outperforms the Jacobi SVD by alarge margin when utilized in stable Green’s function computations. While expected from thefact that QR decompositions are computationally cheaper than SVDs, this difference is evenapparent when the QR factorization is employed in a inversion scheme of higher complexity
involving two (rather than one) matrix decompositions.
In summary, our empirical assessment demonstrates that, among the considered matrix
factorizations and algorithms, the QR decomposition paired with the appropriate inversion
schemes is the most efficient stabilization method for Green’s function calculations in DQMC.
Finally, let us remark that the performance of any stabilization scheme is, in principle,model (parameter) dependent. While a systematic theoretical investigation is beyond the14

SciPost Phys. Core 2, 011 (2020)

scope of this manuscript, we include a brief analysis of a spin-fermion model for antiferromagnetic metallic quantum criticality in App. B. We therefore believe that our major conclusions
bear some universality and can serve as a useful guide.

Acknowledgements
We thank Peter Bröcker, Yoni Schattner, Snir Gazit, and Simon Trebst for useful discussions
and Frederick Freyer for identifying a few typos in an early version of this manuscript. We
acknowledge partial support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Project No. 277101999–TRR 183 (project B01).

AInversion schemes for time slice matrix stacks

In practical DQMC implementations one typically stores intermediate decomposed time slice
matrix products B(τ1 , τ2 ) in a stack for reuse in future equal time Green’s function calculations[7, 8, 37]. In this case, the inversion schemes in Eq. (17) needs to be prefixed by a stableprocedure to combine two elements U L , DL , X L and UR , DR , X R from the stack, corresponding toBl−1 . . . B1 and B M . . . Bl in Eq. 22. Below we describe the latter for both matrix decompositions
considered in the main text11 .

A.1

QR/UDT(StableDQMC.jl: inv_one_plus(::UDT, ::UDT))−1
G = 1 + U L DL TL (UR DR TR )†

−1


= 1 + U L DL TL TR† DR UR† 
|
{z
}

(30)

ud t

= [1 + U DT ]−1,

with U = U L u, D = d, and T = t UR† . This U DT factorization may then be substituted into
Eq. (17).

A.2

SVD

(StableDQMC.jl: inv_one_plus(::SVD, ::SVD))
−1
G = 1 + U L DL VL† UR DR VR†

−1
 
= 1 + U L DL VL† UR DR VR† 
|
{z
}

(31)

ud v †

= 1 + U DV


† −1

,

with U = U L u, D = d, and V = vVR . This SVD factorization may then be substituted into
Eq. (17).
11

In StableDQMC.jl, Julia’s multiple dispatch will automatically select the correct method based on the numberof provided UDT factorizations.

15

SciPost Phys. Core 2, 011 (2020)

150

QR
SVD (Jacobi)
SVD
SVD (D&C)

log singular values

100
50
0
-50
-100
0

10

20

30

40

50

inverse temperature β
Figure (8) Comparison of matrix decompositions in the calculation of the timeslice matrix chain product B M B M −1 · · · B1 near a metallic antiferromagnetic quantum critical point (the spin-fermion model considered in Ref. [1]). Different lines
represent a selection12 of logarithmic singular values as observed in the stabilized
computations. The divide-and-conquer SVD (purple) shows large deviations at low
temperatures β ¦ 8 (∆τ = 0.1) while the other matrix decompositions are in approximate agreement (on top of each other).

B

Spin-fermion model of a metallic quantum critical point

To assess the generality of the findings of the main text, we examine another model in whichitinerant electrons are coupled to an antiferromagnetic order parameter. Specifically, we consider the two-dimensional square lattice spin-fermion model of Ref. [1] hosting a metallic
quantum critical point (QCP), marking the onset of antiferromagnetic spin-density wave order at T = 0. As shown in extensive DQMC studies in Ref. [1], this system exhibits hightemperature superconductivity and features a non-Fermi liquid state in which fermionic quasiparticles lose their coherence.
We will focus our analysis on the vicinity of the QCP (r = −1.74) where quantum criticalfluctuations are most pronounced and interactions are strong and relevant (in the renormalization group sense). The parameters are chosen as for the phase diagram in Fig. 2b of Ref. [1]with the exception of L = 10.
Analogous to Fig. 2 of the main text, Fig. 8 shows logarithmic singular values of the imaginary time slice matrix product chain B M · B M −1 · · · B1 stabilized by various matrix decompositions. Here, in contrast to Fig. 2, we drop the approximation Bi = B and retain the fullimaginary time dependence of the slice matrix factors. While the QR, Jacobi SVD, and regular
SVD variants are capable of accurately capturing all intrinsic scales, the divide-and-conquer
based SVD fails to reliably stabilize the matrix products and displays finite precision artifacts
at low temperatures T = 1/β ® 0.125.
Comparing these results to the findings of the main text — for the one-dimensional Hubbard model — we observe a qualitative deviation: The regular SVD appears to be more stable
for the spin-fermion model. This is in spite of the fact that slice matrices near the metallic
12
The four flavor (two spins, two bands) spin-fermion model with L = 10 leads to 400 singular values si . For
simplicity, Fig. 2 shows only 10 of those, namely s1 , s41 , s81 , ..., s361 .

16

SciPost Phys. Core 2, 011 (2020)

QCP are less well conditioned. A systematic investigation of the implementation details of
LAPACK’s regular SVD with respect to this disparity in accuracy for the two systems would
be desirable, and is left for future work. Importantly, given that the SVD comes with a muchhigher computational cost, the major conclusion of the main text still holds for the strongly
coupled spin-fermion model: Of all considered matrix factorizations, the QR decomposition isclosest to the sweet spot of combined performance and accuracy.

References[1] C. Bauer, Y. Schattner, S. Trebst and E. Berg, Hierarchy of energy scales in an O(3) symmetric antiferromagnetic quantum critical metal: A Monte Carlo study, Phys. Rev. Research 2,023008 (2020), doi:10.1103/PhysRevResearch.2.023008.
[2] Y. Schattner, M. H. Gerlach, S. Trebst and E. Berg, Competing orders ina nearly antiferromagnetic metal, Phys. Rev. Lett. 117, 097002 (2016),
doi:10.1103/PhysRevLett.117.097002.
[3] E. Berg, S. Lederer, Y. Schattner and S. Trebst, Monte Carlo studies of quantum critical metals, Annu. Rev. Condens. Matter Phys. 10, 63 (2018), doi:10.1146/annurev-conmatphys031218-013339.
[4] R. Blankenbecler, D. J. Scalapino and R. L. Sugar, Monte Carlo calculations of coupledboson-fermion systems. I, Phys. Rev. D 24, 2278 (1981), doi:10.1103/PhysRevD.24.2278.
[5] E. Y. Loh, J. E. Gubernatis, R. T. Scalettar, S. R. White, D. J. Scalapino and R. L. Sugar,
Numerical stability and the sign problem in the determinant quantum Monte Carlo method,Int. J. Mod. Phys. C 16, 1319 (2005), doi:10.1142/S0129183105007911.
[6] D. J. Scalapino, S. R. White and S. Zhang, Insulator, metal, or superconductor: The criteria,
Phys. Rev. B 47, 7995 (1993), doi:10.1103/PhysRevB.47.7995.
[7] R. R. dos Santos, Introduction to quantum Monte Carlo simulations for fermionic systems,
Braz. J. Phys. 33, 36 (2003), doi:10.1590/S0103-97332003000100003.
[8] F. Assaad, Quantum Monte Carlo methods on lattices: The determinantal approach, Lecture
Notes, J. Grotendorst, D. Marx, A. Muramatsu (Eds.), John von Neumann Institute for
Computing, Jülich, NIC Series 10, ISBN 3000090576 (2002).
[9] E. Y. Loh, J. E. Gubernatis, R. T. Scalettar, S. R. White, D. J. Scalapino and R. L. Sugar,
Sign problem in the numerical simulation of many-electron systems, Phys. Rev. B 41, 9301
(1990), doi:10.1103/PhysRevB.41.9301.
[10] M. Troyer and U.-J. Wiese, Computational complexity and fundamental limitations
to fermionic quantum Monte Carlo simulations, Phys. Rev. Lett. 94, 170201 (2005),
doi:10.1103/PhysRevLett.94.170201.
[11] J. E. Hirsch, Two-dimensional Hubbard model: Numerical simulation study, Phys. Rev. B
31, 4403 (1985), doi:10.1103/PhysRevB.31.4403.
[12] S. R. White, D. J. Scalapino, R. L. Sugar, E. Y. Loh, J. E. Gubernatis and R. T. Scalettar,
Numerical study of the two-dimensional Hubbard model, Phys. Rev. B 40, 506 (1989),
doi:10.1103/PhysRevB.40.506.

17

SciPost Phys. Core 2, 011 (2020)[13] A. Moreo and D. J. Scalapino, Two-dimensional negative-U Hubbard model, Phys. Rev.
Lett. 66, 946 (1991), doi:10.1103/PhysRevLett.66.946.
[14] M. Hohenadler, Z. Y. Meng, T. C. Lang, S. Wessel, A. Muramatsu and F. F. Assaad, Quantum phase transitions in the Kane-Mele-Hubbard model, Phys. Rev. B 85, 115132 (2012),
doi:10.1103/PhysRevB.85.115132.
[15] M. H. Gerlach, Y. Schattner, E. Berg and S. Trebst, Quantum critical properties of a metallic spin-density-wave transition, Phys. Rev. B 95, 035124 (2017),
doi:10.1103/PhysRevB.95.035124.
[16] Y. Schattner, S. Lederer, S. A. Kivelson and E. Berg, Ising nematic quantum critical point in a metal: A Monte Carlo study, Phys. Rev. X 6, 031028 (2016),
doi:10.1103/PhysRevX.6.031028.
[17] S. Gazit, M. Randeria and A. Vishwanath, Emergent Dirac fermions and broken symmetries in confined and deconfined phases of Z2 gauge theories, Nat. Phys. 13, 484 (2017),
doi:10.1038/nphys4028.
[18] E. Y. Loh, J. E. Gubernatis, R. T. Scalettar, R. L. Sugar and S. R. White, Stable
matrix-multiplication algorithms for low-temperature mumerical simulations of fermions,in Interacting Electrons in Reduced Dimensions, Springer US, Boston, MA, ISBN
9781461278696 (1989), doi:10.1007/978-1-4613-0565-1_8.
[19] Z. Bai, C. Lee, R.-C. Li and S. Xu, Stable solutions of linear systems involving
long chain of matrix multiplications, Linear Algebra Its Appl. 435, 659 (2011),
doi:10.1016/j.laa.2010.06.023.
[20] S. Sorella, S. Baroni, R. Car and M. Parrinello, A novel technique for the simulation of interacting fermion systems, Europhys. Lett. 8, 663 (1989), doi:10.1209/02955075/8/7/014.
[21] J. Bezanson, A. Edelman, S. Karpinski and V. B. Shah, Julia: A fresh approach to numericalcomputing, SIAM Rev. 59, 65 (2017), doi:10.1137/141000671.
[22] J. M. Perkel, Julia: come for the syntax, stay for the speed, Nature 572, 141 (2019),
doi:10.1038/d41586-019-02310-3.
[23] C. Rackauckas, ODE solver multi-language wrapper package work-precision benchmarks (MATLAB, SciPy, Julia, deSolve (R)), https://benchmarks.juliadiffeq.org/html/
MultiLanguage/wrapper_packages.html, [Online; accessed 4-March-2020].
[24] X.-Z. Luo, J.-G. Liu, P. Zhang and L. Wang, Yao.jl: Extensible, efficient framework for quantum algorithm design (2019), arXiv:1912.10877.
[25] J. E. Hirsch, Discrete Hubbard-Stratonovich transformation for fermion lattice models,
Phys. Rev. B 28, 4059 (1983), doi:10.1103/PhysRevB.28.4059.
[26] S. Sachdev, Quantum phase transitions, In Quantum Phase Transitions, Cambridge University Press, ISBN 9780521514682 (2011).
[27] H. F. Trotter, On the product of semi-groups of operators, Proc. Am. Math. Soc. 10, 545(1959), doi:10.2307/2033649.
[28] M. Suzuki, Quantum statistical monte carlo methods and applications to spin systems, J.
Stat. Phys. 43, 883 (1986), doi:10.1007/BF02628318.
18

SciPost Phys. Core 2, 011 (2020)[29] Z.-X. Li and H. Yao, Sign-problem-free fermionic quantum Monte Carlo: Developments andapplications, Annu. Rev. Condens. Matter Phys. 10, 337 (2019), doi:10.1146/annurevconmatphys-033117-054307.
[30] D. Goldberg, What every computer scientist should know about floating-point arithmetic,
ACM Comput. Surv. 23, 5 (1991), doi:10.1145/103162.103163.
[31] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz,
A. Greenbaum, S. Hammarling, A. McKenney and D. Sorensen, LAPACK Users’ Guide, Society for Industrial and Applied Mathematics, Philadelphia, ISBN 0-89871-447-8 (1999).
[32] J. Dongarra, M. Gates, A. Haidar, J. Kurzak, P. Luszczek, S. Tomov and I. Yamazaki, The
singular value decomposition: Anatomy of optimizing an algorithm for extreme scale, SIAM
Rev. 60, 808 (2018), doi:10.1137/17M1117732.
[33] M. Bercx, F. Goth, J. S. Hofmann and F. Assaad, The ALF (Algorithms for Lattice Fermions)project release 1.0. Documentation for the auxiliary field quantum Monte Carlo code, SciPost
Phys. 3, 013 (2017), doi:10.21468/SciPostPhys.3.2.013.
[34] R. C. Lee, S. Chiesa, C. N. Varney, E. Khatami, Z. Bai, E. F. D’Azevedo, M. Jarrell, T. Maier,
S. Y. Savrasov, R. T. Scalettar and K. Tomko, Quest: Quantum electron simulation toolbox(2010), https://web.cs.ucdavis.edu/~bai/PETAMAT/scidac10paper.pdf .
[35] S. Blackford, Error bounds for the singular value decomposition, http://www.netlib.org/
lapack/lug/node96.html, [Online; accessed 4-March-2020] (1999).
[36] J. E. Hirsch, Stable Monte Carlo algorithm for fermion lattice systems at low temperatures,Phys. Rev. B 38, 12023 (1988), doi:10.1103/PhysRevB.38.12023.
[37] P. Broecker and S. Trebst, Numerical stabilization of entanglement computation in
auxiliary-field quantum Monte Carlo simulations of interacting many-fermion systems,
Phys. Rev. E 94, 063306 (2016), doi:10.1103/PhysRevE.94.063306.

19