SciPost Phys. 8, 012 (2020)

Sparse sampling and tensor network representation of
two-particle Green’s functions
Hiroshi Shinaoka1? , Dominique Geffroy2,3 , Markus Wallerberger4,3 , Junya Otsuki5 ,Kazuyoshi Yoshimi6 , Emanuel Gull4 and Jan Kuneš3,7
1 Department of Physics, Saitama University, Saitama 338-8570, Japan
2 Department of Condensed Matter Physics, Faculty of Science,
Masaryk University, Kotlářská 2, 611 37 Brno, Czech Republic
3 Institute of Solid State Physics, TU Wien, 1040 Vienna, Austria
4 University of Michigan, Ann Arbor, Michigan 48109, USA
5 Research Institute for Interdisciplinary Science,
Okayama University, Okayama 700-8530, Japan
6 Institute for Solid State Physics, University of Tokyo, Chiba 277-8581, Japan
7 Institute of Physics, Czech Academy of Sciences,
Na Slovance 2, 182 21 Praha 8, Czech Republic? shinaoka@mail.saitama-u.ac.jp

Abstract
Many-body calculations at the two-particle level require a compact representation of
two-particle Green’s functions. In this paper, we introduce a sparse sampling scheme
in the Matsubara frequency domain as well as a tensor network representation for twoparticle Green’s functions. The sparse sampling is based on the intermediate representation basis and allows an accurate extraction of the generalized susceptibility from a
reduced set of Matsubara frequencies. The tensor network representation provides asystem independent way to compress the information carried by two-particle Green’s
functions. We demonstrate efficiency of the present scheme for calculations of staticand dynamic susceptibilities in single- and two-band Hubbard models in the frameworkof dynamical mean-field theory.
Copyright H. Shinaoka et al.
This work is licensed under the Creative Commons
Attribution 4.0 International License.
Published by the SciPost Foundation.

Received 18-09-2019
Accepted 15-01-2020
Check forupdatesPublished 28-01-2020
doi:10.21468/SciPostPhys.8.1.012

Contents
1 Introduction

22 Intermediate Representation (IR) for Green’s functions
2.1 One-particle Green’s function
2.2 General form of IR for Two-particle Green’s functions
2.3 Simplified form for fixed bosonic frequency
2.4 Graphical representation

3
3
4
4
6

3 Sparse sampling

6
1

SciPost Phys. 8, 012 (2020)

4 Tensor network representation
4.1 Low-rank tensor decomposition
4.2 Fitting algorithm

88
95 DMFT calculations for single-band Hubbard model on a square lattice10

6 DMFT calculations for two-band Hubbard model

10

7 Summary

15

A Intermediate representation at fixed bosonic frequency

16

B Optimization algorithm for tensor regression
B.1 Alternating least squares
B.1.1 Optimization of x (0)B.1.2 Optimization of x (1)B.1.3 Optimization of x (2)B.1.4 Optimization of x (3)B.2 Acceleration techniques and convergence condition
B.3 Technical details and numerical results

17
17
17
18
18
18
18
19

References

20

1

Introduction

Two-particle (2P) Green’s functions are building blocks of a variety of many-body theories [1].
They are a key element for calculation of susceptibilities in the framework of dynamical meanfield theory (DMFT) [2] as well as for diagrammatic calculations including vertex corrections
and diagrammatic extensions of DMFT [3–16]. In equilibrium finite-temperature formalism,
the 2P quantities depend on three Matsubara frequencies, the low- and high-frequency partsof which must be treated accurately. Storing 2P quantities alone is a challenge, more so when
multiple orbitals, low symmetries or low temperatures are involved.
Several methods have been proposed to address the storage issue. Conventional approaches
are based on a separate treatment of the low- and high-frequency parts [17–21]. The frequency dependence is treated exactly in a small low-frequency box while in the outside region
an asymptotic form is used. This works efficiently at relatively high temperatures. As thetemperature is lowered the size of the low-frequency box grows, until it becomes prohibitively
large.
Recently, the intermediate representation (IR) basis was introduced as a promising solution
to the storage issue [22]. In IR the size of the data grows only logarithmically with the inverse
temperature β and the bandwidth. A fitting procedure allows IR expansion of numerical datain the Matsubara frequency domain. Nevertheless, two obstacles remain: the computational
cost of the IR expansion and the size of the IR tensor.
Regarding the first obstacle, the input data for the fitting scheme is very large, having a
dense support in the Matsubara frequency domain. The fitting procedure thus becomes prohibitively expensive at low T . As for the second one, the IR represents a 2P quantity as ahigh-order tensor involving spin and orbital dimensions, in addition to those for the IR basis itself. Further compactification of the tensor is required for solving realistic multi-orbital
2

SciPost Phys. 8, 012 (2020)

systems at low T .
In this paper, we address these two issues. First, we introduce a sparse grid in the Matsubara frequency domain, which contains the desired information about the 2P Green’s functions
This extends the approach developed in Ref. [23] for one-particle (1P) Green’s functions. We
introduce an efficient tensor network representation of the IR tensor and a fitting (regression)algorithm to determine it. Reduction of the data to be sampled thanks to the sparse grid makes
evaluation of the IR coefficients very efficient and solves the first issue. The tensor regression
provides a model-independent way to compress the IR tensor and tackles the second issue.
We demonstrate the performance of the present method in the context of DMFT. First,
we test the accuracy of sparse sampling and tensor network representation by calculating thestatic susceptibility of the single-band Hubbard model on a square lattice. Next, we show
the efficiency of the present method for dynamical susceptibility calculation for a two-band
Hubbard model with low symmetry.
The paper is organized as follows. In the next section, the IR for 1P and 2P Green’s functions is reviewed. Sparse sampling of 2P Green’s functions is introduced in Sec. 3. In Sec. 4,
the tensor network representation is presented. Its accuracy for computing of static susceptibilities of the single-band Hubbard model is demonstrated in Sec. 5. In Sec. 6, we present
numerical results for dynamical susceptibility calculations, in the more demanding context ofan ordered phase of the two-band Hubbard model. In Sec. 7, we summarize and conclude.

2

Intermediate Representation (IR) for Green’s functions

Here we review the IR for 1P and 2P Green’s functions introduced in [22] and [24]. The reader
may refer to Section 7 of [25] for a review.

2.1

One-particle Green’s function

The IR for 1P Green’s functions was introduced in Ref. [24]. The spectral (Lehmann) representation of the 1P Green’s function G(τ) in the imaginary-time domain reads
Z ωmaxdωK α (τ, ω)ρ(ω),G(τ) = −

(1)

−ωmax

where we assume ħ
h = 1. The superscript α specifies statistics: α = F for fermion and α = B
for boson. The spectrum ρ(ω) is assumed to be bounded within the interval [−ωmax , ωmax ].
The kernel K α (τ, ω) reads
K α (τ, ω) ≡ ωδα,Be−τω
1 ± e−βω(2)

for 0 ≤ τ ≤ β (β is the inverse temperature). Here, the + and − signs are used for fermionsand bosons, respectively. The extra ω factor for bosons in Eq. (2) is introduced in order toavoid a singularity of the kernel at ω = 0.
For fixed values of ωmax and β, the IR basis functions are defined through the singularvalue decomposition (SVD)
K α (τ, ω) =∞X

Slα Ulα (τ)Vlα (ω),

(3)

l=0

Zβ
with
0

dτUlα (τ)Ulα0 (τ) =Z ωmax−ωmax

3dωVlα (ω)Vlα0 (ω) = δll 0 ,SciPost Phys. 8, 012 (2020)

where the singular values Slα (> 0) decrease with increasing l exponentially.
For fermions, the Green’s function can be expanded as∞
X

G(τ) =

G (l)UlF (τ),

(4)

l=0
G (l) = −SlF ρl ,(5)Rωwhere ρl ≡ −ωmax dωρ(ω)VlF (ω). The exponential decay of SlF ensures a fast decay of Gl ifmax
the spectrum is bounded in [−ωmax , ωmax ]. The accuracy of the expansion can be controlled
by applying a cut-off for the singular values. The Matsubara-frequency representation of theGreen’s function reads
Zβ∞X
dτG(τ)e iωn τ =G(iωn ) =0

l=0with UlF (iωn ) ≡ZβG (l)UlF (iωn )(6)dτUlF (τ)e iωn τ .

0

2.2

General form of IR for Two-particle Green’s functions

The problem of the three- and four-point Green’s functions was considered in Ref. [22]. The
four-point Green’s function can be expressed in the Matsubara domain as
Gi jkl (iωn , iωn0 , iνm )Zβ≡ β −2

0

≡

16 X
X
r=1 l1 l2 l3dτ1 dτ2 dτ3 dτ4 e i(ωn τ12 +ωn0 τ34 +νm τ14 ) 〈Tτ ci† (τ1 )c j (τ2 )ck† (τ3 )cl (τ4 )〉00

0

GO (r, l1 , l2 , l3 )Ulα1 (iω)Ulα2 (iω0 )Ulα3 (iω00 ),

(7)
(8)

where ω, ω0 , ω00 stand for r-dependent combinations of ωn , ωn0 , and νm listed in Table 1.
Similarly, the indices α, α0 and α00 take the value F or B̄ depending on r as indicated in Table 1.
The indices i, j, k, l denote the flavor (combined spin and orbital), while O is a composite index
representing the quadruplet (i, j, k, l). Formally, the tensor G contains full information about
G (meaning in particular, information about its values at all bosonic and fermionic frequencies
for all flavors).

2.3

Simplified form for fixed bosonic frequency

We now derive a variant of Eq. (8), which holds for a fixed bosonic frequency. The principleof the derivation is the same as above. Nevertheless, if one is interested in only a few bosonicfrequencies, the reduced number of degrees of freedom allows us to represent the 2P Green’s
function as a tensor of lower rank. In the conventional particle-hole notation, the Green’s
function depends only on two fermionic frequencies for a fixed bosonic frequency. It is shownin Appendix A that in this case, this frequency dependence can be written as
GO (iωn , iωn0 , iνm ) =12 X
X
r=1 l1 ,l2

(r)

(r)

1

2

GO (r, l1 , l2 ; iνm )Ul (iω)Ul (iω0 ),(9)where ωn and ωn0 are fermionic frequencies, and νm is a bosonic frequency. The index r
relates to the 12 distinct representations generated from the three terms in Eq. (18). The
4

SciPost Phys. 8, 012 (2020)

Table 1: 16 different notations of Matsubara frequencies and statistics for the fourpoint Green’s functions G. Adapted from Ref. [22].
r(iω, iω0 , iω00 )(α, α0 , α00 )

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

(iνm + iωn , −iωn , iωn0 )(iνm + iωn , −iωn , −iνm − iωn0 )(iνm + iωn , iωn0 , −iνm − iωn0 )(−iωn , iωn0 , −iνm − iωn0 )(iνm + iωn , iνm , iνm + iωn0 )(iνm + iωn , iνm , −iωn0 )(iνm + iωn , iνm + iωn + iωn0 , iνm + iωn0 )(iνm + iωn , iνm + iωn + iωn0 , iωn )(iνm + iωn , iωn − iωn0 , −iωn0 )(iνm + iωn , iωn − iωn0 , iωn )(−iωn , iνm , iνm + iωn0 )(−iωn , iνm , −iωn0 )(−iωn , −iωn + iωn0 , iνm + iωn0 )(−iωn , −iνm − iωn − iωn0 , −iωn0 )(iωn0 , iνm + iωn + iωn0 , iνm + iωn0 )(iωn0 , −iωn + iωn0 , iνm + iωn0 )(F,F,F)(F,F,F)(F,F,F)(F,B,F)(F,B,F)(F,B,F)(F,B,F)(F,B,F)(F,B,F)(F,B,F)(F,B,F)(F,B,F)(F,B,F)(F,B,F)(F,B,F)(F,B,F)

W

U

r

r

W

U

r
l1

l2

W
O
r

l1

l2

GO (r, l1 , l2 )Figure 1: Graphical representation for Eq. (9). Each circle or rectangle with N legs
represents a N -way tensor. When two tensors share a leg, the summation over thecorresponding index must be taken. A rectangle with a diagonal line represents an
identity tensor. The N -way identity tensor t is defined as t i1 ···iN = 1 iff i1 = · · · = iN ,and t i1 ···iN = 0 otherwise.

5

SciPost Phys. 8, 012 (2020)

Table 2: 12 different notations of Matsubara frequencies and statistics for four-point
Green’s functions in the particle-hole notation at a fixed bosonic frequency.
r(iω, iω’)(α, α0 )

1
2
3
4
5
6
7
8
9
10
11
12

(iωn , iωn0 )(iωn − iωn0 , iωn0 )(iωn0 − iωn , iωn )(iωn , iωn0 + iνm )(iωn − iωn0 , iωn0 + iνm )(iωn0 − iωn , iωn + iνm )(iωn + iνm , iωn0 )(iωn + iωn0 + iνm , iωn0 )(iωn0 + iωn + iνm , iωn )(iωn + iνm , iωn0 + iνm )(iωn + iωn0 + iνm , iωn0 + iνm )(iωn0 + iωn + iνm , iωn + iνm )(F, F)(B, F)(B, F)(F, F)(B, F)(B, F)(F, F)(B, F)(B, F)(F, F)(B, F)(B, F)

imaginary-time frequencies (iω, iω0 ) and the statistics of the basis functions depend on r as
summarized in Table 2.
This expression is formally similar to Eq. (8), being meant to store the full information fora single bosonic frequency. The same sparse sampling strategy, which we shall introduce inSec. 3, can therefore be employed in both situations.

2.4

Graphical representation

For the sake of clarity, we introduce a graphical representation for the tensor operations involving 2P Green’s functions. As an example, we present in Fig. 1 the diagram corresponding
to the right-hand side of Eq. (9). Each rectangle/circle represents a tensor whose indices aredenoted by legs (the nature of the shape does not matter). The set of expansion coefficientsGO (r, l1 , l2 ; iνm ) appears as a purple four-legged box, representing a rank four tensor. The
green rectangle represents the basis functions terms in Eq. (9). A detailed view of their action
is shown in the upper panel, where W is introduced as a composite index of (iω, iω0 ).

3Sparse sampling

The sparse sampling scheme was originally proposed for a 1P Green’s function in Ref. [23].
For fermions, the expansion of G(iωn ) reads
G(iωn ) =

N
l −1
X
l=0

G (l)UlF (iωn ),(10)

where the number of coefficients Nl = lmax +1 determines the accuracy of the expansion. It was
shown that the full frequency dependence of a 1P Green’s function can be reconstructed fromthe values of the Green’s function on a carefully chosen sparse subset of sampling points in the
Matsubara domain. UlF (iωn ) is real (odd l) or pure imaginary (even l), and oscillates aroundzero. The procedure described in Ref. [23] is based on picking the positions of the extrema of|UlF (iωn )|. The procedure generates Nl (even l) or Nl + 1 (odd l) sampling points. The same
max
procedure generates Nl + 1 (even l) or Nl (odd l) sampling points for bosons.

6

SciPost Phys. 8, 012 (2020)m = −20

400

200

200

0

0

n0

n0

400

−200
−400
−400 −200

−200
0
n

200

−400
−400 −200

400

m = 10

400

200

200

0

0

n0

n0

400

−200
−400
−400 −200

m=0

0
n

200

400

200

400

m = 20

−200
0
n

200

−400
−400 −200

400

0
n

Figure 2: Sampling points generated in the three-frequency space for Λ = 104 and
Nl = 24. We show cuts at m = −20, 0, 10, 20. We show only sampling points at low
frequencies. The actual sampling points are distributed up to |n|, |n0 | ® 1800 and|m| ® 2200.
m = −20

200

100

100

0

0

n0

200

−100

−100

−200

−200

−300
−300 −200 −100 0
n
300

n0

300

m = 10

300

200

200

100

100

0

0

−100

−100

−200

−200

−300
−300 −200 −100 0
n

m=0

−300
−300 −200 −100 0n

100 200 300

n0

n0

300

m = 20−300−300 −200 −100 0n

100 200 300

100 200 300

100 200 300

Figure 3: Sampling points generated in the particle-hole notation independently form = −20, 0, 10, 20 (Λ = 100 and Nl = 19).
7

SciPost Phys. 8, 012 (2020)

We extend this procedure in a straightforward manner for the expansion of the 2P Green’s
function. Each summand indexed by r in Eq. (8) (Eq. (9) ) is handled in turn. For a fixed valueof r, the sets of sampling points relative to each factor in the corresponding product of basis
functions are built. Then, the triplets (pairs) of direct products of such sets are determined,
and make up the r th set of points to be sampled, in Z3 (Z2 ).
For simplicity of implementation, we use the same Nl for fermions and for bosons in the
expansion. In general, for a given value of Λ ≡ βωmax , the singular values decay slower for
fermions than for bosons. In practice, we determine Nl based on a given singular-value cutoff
for fermions, and use the same Nl for bosons.
As an illustration, for Eq. (8), we obtain Nl3 + O(1) sampling points. The final set ofsampling points we need to consider is the union of the sets of sampling points obtained for
all r. The size of the union is less than the sum of the sizes of the individual sets, thanks to
the overlap between them (in particular at low frequencies). In Sec. 6, we will use Λ = 104
and Nl = 24 (cutoff value 10−4 for singular values). For this parameter set, the procedure
generates 165 912 sampling points for Eq. (8), which is slightly smaller than 16Nl3 = 221 184
due to the overlap. Figure 2 shows the distribution of these sampling points. One can see that
their distribution is more dense at low frequencies, getting sparse at high frequencies.
Figure 3 shows the distribution of the sampling points generated for Eq. (9) with Nl = 19
(cutoff value of 10−5 ). We obtain 2 972, 1336, 2 516, 2 972 sampling points for m = −20, 0, 10,
20, respectively. We will use these parameters in Sec. 5.
One technical caveat needs to be pointed out: the expansion of the 2P Green’s functions
involves the so-called “extended” bosonic basis set [22]. A basis function UlB from this set only
exhibits max(0, l − 2) sign changes, due to the extra basis functions at l = 0, 1. The procedure
above would thus yield lmax − 1 = Nl − 2 or Nl − 1 sampling points for this basis set. Therefore
the actual process is slightly altered from the above description. The sampling points relative
to the extended bosonic basis are generated from the extrema of UlB +2 instead of UlB . Thismax
max
ensures that the number of unknown coefficients matches the number of sampling points.
In the following sections, we will demonstrate that the sampling on the sparse grid issufficient to evaluate the 2P Green’s function with the desired precision for any Matsubara
frequency.

4Tensor network representation

In this section, we introduce an efficient fitting algorithm based on a tensor network representation for the IR tensor. We refer the reader who is not familiar with tensor networks to
Refs. [26–28]. In principle, numerical data on the sampling points can be fitted using either
Eq. (8) or Eq. (9) by using the least squares method. The computational load of this naive ap4
4
proach scales as O(Nl6 Norb
) or O(Nl9 Norb
) for two- and three-frequency quantities, respectively.
Here Norb is the number of orbitals, and Nl grows logarithmically with respect to β. The fitting
rapidly becomes too costly at low temperature, especially for three-frequency quantities.

4.1

Low-rank tensor decomposition

We introduce the following low-rank decomposition of G:

GO (r, l1 , l2 , l3 ) ' G̃O (r, l1 , l2 , l3 ) ≡DX(0) (1) (2) (3) (4)x d r x dl x dl x dl x dO ,
1

2

3

(11)

d=1

where r runs over different representations and O ≡ (i, j, k, l). This type of tensor decomposition is known as a Canonical Polyadic (CP) decomposition and is widely used in many fields,
8

SciPost Phys. 8, 012 (2020)

r

l1

l2

O

G(r, l1 , l2 ; O)

'

r

l1

l2

O

x (0)

x (1)

x (2)

x (3)

d

Figure 4: Graphical representation of the tensor decomposition in Eq. (12). Eachcircle or rectangle with N legs represents a N -way tensor. When two tensors share aleg, the summation over the corresponding index must be taken. A rectangle with adiagonal line represents a identity tensor.
e.g., for accelerating quantum chemistry calculations by factorizing Coulomb integrals [29–31]. For the simplified form with fixed bosonic frequency in Eq. (9), the low-rank decomposition reads

GO (r, l1 , l2 ; iνm ) ' G̃O (r, l1 , l2 ; iνm ) ≡

D
X

(0)

(1)

(2)

1

2

(3)

x d r (νm )x dl (νm )x dl (νm )x dO (νm ),

(12)

d=1

which is illustrated in Figure 4.
Both expressions become exact for a large enough D. The decomposition is beneficial forthe fitting procedure if a good approximation of the full tensor is obtained for a reasonably
small value of D. We discuss further in the text how this condition can be checked numerically.
Note that the dependence on the orbital indices is not decomposed and they still appear as
the composite index O in Eqs. (11) and (12). In the CP decomposition, we do not assume any
orthogonality conditions for the decomposed tensors, unlike the SVD of a matrix. Recently,
some of the authors and co-workers have proposed strong-coupling formulas for computing
momentum dependent susceptibilities in DMFT [32]. Their formulas can be regarded as a
special case of Eq. (12) with D = 1 and x d r = δ r1 . Our preliminary results indicate that thefurther decomposition of the dependence on O into individual spin/orbital indices requires
typically an even larger D, which is not beneficial.

4.2

Fitting algorithm

We explain how to fit some existing data of the 2P Green’s function on the sampling points
in the Matsubara domain using Eq. (11) or Eq. (12) without explicitly constructing the big
tensors on the left-hand side. The x tensors in the equation can be regarded as free fitting
parameters. For instance, for Eq. (11), we define a cost function for the fitting asf ({ x(i)}) = ||G(W , O) − G̃(W , O)[x(i)]||22 + α

4X||x (i) ||22 ,

(13)

i=0

where || . . . ||2 denotes the Frobenius norm and W runs over sampling points in the Matsubara
frequency domain. G(W , O) are the data of the Green’s function on the sampling points which
we fit, while G̃(W , O)[x (i) ] denotes the data of the Green’s function evaluated from { x (i) }.
We introduce the small parameter α > 0 in order to regularize the optimization problem.
Without this parameter, the problem would be ill-posed due to the overcompleteness of therepresentation. We have not observed any visible systematic errors in interpolated data forsmall values of α, i.e., α ≤ 10−5 . We used α = 10−8 and 10−5 in Sec. 5 and Sec. 6, respectively.
9

SciPost Phys. 8, 012 (2020)

The minimization of this cost function is a non-convex optimization problem. We foundthat despite its non-convex nature, the cost function can be minimized efficiently using standard methods starting from randomly initialized parameters. In some cases, we observed theexistence of multiple solutions, being different only slightly in terms of the cost function. Thisissue thus does not matter in practice. We refer the interested reader to Appendix B for moredetails on the optimization method.

5

DMFT calculations for single-band Hubbard model on a squarelattice

As a test bed for our method, we first consider a single-band Hubbard model on a square latticeat half filling for U = 12t (1.5× the bandwidth W = 8t), where the hopping t = 1 sets the unitof energy. The inverse temperature β = 2.5 is slightly above the antiferromagnetic transition.
We use (approximate) Hubbard-I solver, which provides semi-analytic representation of local
susceptibilities and thus allows precise analysis of our data compression approach. We firstcompute and fit the local (impurity) generalized susceptibility X loc by subtracting the relevant
disconnected parts from the local 2P Green’s function. Interpolating the local generalized
susceptibility in the Matsubara frequency domain, we solve the Bethe-Salpeter equation (BSE)and compute the static DMFT susceptibilities of the model.
We compute X loc for all spin-orbital components on 1 336 sampling points generated forΛ = 100 and Nl = 19 at zero bosonic frequency m = 0. The distribution of the sampling
points is shown in the right top panel of Fig. 3. Then, we fit the data using Eq. (12) (G issimply replaced by X loc ). Figure 5 shows how the fitting errors decay as D is increased. We
found that the residual of the fit vanishes quickly with increasing D. Figure 6 compares theexact and interpolated values of the local susceptibility in Matsubara frequency space. One
can see that for D = 15, the fit matches the exact values on the sampling points and preciselyinterpolates the data. Increasing D further does not improve the fit substantially, which may
be due to the truncation of the basis.
The inversion of the Bethe-Salpeter equations (for the determination of the lattice susceptibilities) directly in the IR and tensor network format is still an open question (see the discussion
in Sec. 7). In this study, we execute the inversion itself using the Matsubara representation,
based on the interpolation of the generalized susceptibility in a box of width [−Nω , Nω −1] forfermionic frequencies. Corrections from higher frequencies outside the box are treated using
the procedure described in Appendix B of Ref. [32]. The corresponding physical quantities are
obtained by summation over the fermionic frequencies.
Figure 7 shows the physical lattice susceptibility using the fitted results for D = 5 andiν = 0. For the calculations in this section, we took Nω = 100. We found that the results for
D = 5 are already indistinguishable from the exact ones at the scale of the figure. There is apronounced peak at M = (π, π), corresponding to an antiferromagnetic spin order. The otherthree eigenvalues, which correspond to charge susceptibility, are not enhanced.
The storage of X loc computed on the sampling points takes up 340 kB (including all spin
sectors), while the compressed data for D = 5 only 5.2 kB.

6

DMFT calculations for two-band Hubbard model

Next, we apply the present IR approach to a state-of-art linear response DMFT calculation. To
this end, we choose the two-band Hubbard model in an intermediate coupling regime and low
spontaneously broken symmetry, a problem that some of us have studied recently [33]. The
10

SciPost Phys. 8, 012 (2020)

10−1
10−2
|∆X loc |max
|X loc |max

10−3
10−4
10−5
10−6

0

10

20

30

D

Figure 5: Fitting errors of X loc for the single-band Hubbard model on the squarelattice at U = 1.5W and β = 2.5.
model Hamiltonian reads
X
X
X
 ∆X
†
†
b
b
a jσ + t b biσ
b jσ +
t a aiσ(naiσ − niσ)+Unαi↑ nαi↓ +(U 0 − Jδσσ0 )naiσ niσ
H=
0,
2
i j,σ
i,σ
i,α
i,σσ0
†
†
where aiσand biσare fermionic operators that create electrons with the respective orbital
flavors and spin σ at site i of a square lattice. The first term describes the nearest neighbor†
ciσ , captures the crystalhopping. The rest, expressed in terms of local densities nci,σ ≡ ciσfield ∆, the Hubbard interaction U and Hund’s exchange J in the Ising approximation. We
use the same hopping parameters t a = 0.4118, t b = −0.1882 as in [33], but choose weaker
interaction U = 2 (J = U/4, U 0 = U − 2J) and ∆ = 1.55. At the sstudied inverse temperature
β = 60, an ordered phase called polar excitonic condensate [34] is realized. We follow thesame the algorithmic programme as in Ref. [33], except for the representation of the 2P Green’s
function.
The local 2P Green’s function G loc is sampled on a non-uniform grid in the Matsubara
frequency domain, shown in Fig. 2 (Λ = 104 , Nl = 24), using a modified version of the
ALPS/CT-HYB impurity solver [35, 36] based on the continuous-time hybridization expansion
algorithm [37, 38]. The regression (11,13) for G loc provides us with the x (i) , i ∈ {0, . . . , 4}coefficients, Eq. (11), which in turn are used to interpolate the local 2P Green’s function at all
Matsubara frequencies. The generalized local susceptibility X loc is then evaluated by subtracting the disconnected part from the local 2P Green’s function.
Figure 8 illustrates the convergence of the process with increasing value of D in this situation. The solid lines represent the real part of G loc in the Matsubara representation, obtained
from Eq. 8. Both panels show the profile for a fixed bosonic frequency in the particle-holenotation (m = 0 and m = 10), along a cut in the two-dimensional fermionic frequency space,slightly away from the main diagonal (shifted by ten Matsubara frequencies). The crosses represent the actual sparse frequencies, where the data was sampled. The structure of the data
at zero bosonic frequency is relatively simple, so that the fit is excellent for D = 30, but thetenth bosonic frequency requires D = 60.
Using the above data we can solve the BSE to get the linear response of our system. As an

11

SciPost Phys. 8, 012 (2020)(a) O = (↑, ↓, ↑, ↓)(b) O = (↑, ↑, ↑, ↑)

n0

0
−20
−20

0
n

20

0.02

20

0.01
0

n0

0.06
0.04
0.02
0.00
−0.02
−0.04
−0.06

20

0.00
−0.01

−20
−20

n0 = 0

10−4

10−8
−100

|X loc|10−4|X loc|

10−2

10

D=5
D = 10
D = 15
Exact

10

0
n

−6

10−8
−100

100

n = n0

10−4

10−8
−100

D=5
D = 10
D = 15
Exact

0
n

0
n

|X loc|10−4|X loc|

10−2

10

D=5
D = 10
D = 15Exact

100

n = n0

10−2

−6

20

n0 = 0

10−2

−6

0
n

−0.02

10

10−8
−100

100

D=5
D = 10
D = 15
Exact

−6

0
n

100

Figure 6: Comparison of interpolated results and exact values of the generalized
susceptibility X loc for the single-band Hubbard model. The panels (a) and (b) show
the data for O = (↑, ↓, ↑, ↓) and (↑, ↑, ↑, ↑), respectively. The 2D colormaps show theinterpolated data for D = 15.
example we calculate the diagonal susceptibilities for local operators [33]
X 

Sz =
n↑c − n↓c ,
c=a,b

Rφγ =
I φγ =

X
αβ

X
αβ



γ
σαβ aα† bβ + bα† aβ ,γiσαβ aα† bβ − bα† aβ ,which capture the low-energy dynamics of the polar condensate. We have chosen the ordered
12

SciPost Phys. 8, 012 (2020)

6

D=5
Exact

5

χ

4
3
2
1
0
X

Γ

M

X

M0

Figure 7: Physical lattice susceptibility χ(q) computed for the single-band Hubbard
model on the square lattice at U = 1.5W and β = 2.5. We used D = 5. Different
branches correspond to different eigen modes of χ(q , 0), i.e., the spin susceptibility
and the charge susceptibility.

(a) O = (a ↑, a ↑, a ↑, a ↑), m = 0
100
60

0

0.0

−20
−40

−0.2

−80

−0.4

−60
−100
−100−80−60−40−20 0 20 40 60 80 100n

Re X loc

20
n0

0.4

0.2

40

D=8
D = 30D = 60

0.6

0.4

80

0.2
0.0
−0.2
−0.4

−100

−50

0
n

50

100

(b) O = (b ↓, b ↓, a ↑, a ↑) , m = 10
100
80

0.015

60
0.02

n0

20
0
−20

0.00

−40

−0.02

−80

−0.04

−60
−100
−100−80−60−40−20 0 20 40 60 80 100
n

0.010
Re X loc

40

D=8
D = 30D = 60

0.04

0.005
0.000
−0.005
−100

−50

0
n

50

100

Figure 8: Local 2P Green’s function G loc interpolated at zero bosonic frequency m = 0[(a)] and a finite bosonic frequency m = 10 [(b)] for the two-band Hubbard model.
The 2D color maps show results of the local 2P Green’s function for D = 60. The
right-hand halves of the panels (a) and (b) show the convergence of the interpolated
values with respect to D on the line of n0 = n + 10.

13

SciPost Phys. 8, 012 (2020)

Sz

y

40
30
20
10

8
6
4
2
(, )

(0, 0)( , 0)y(, )(0, 0)( , 0)

4
3
2
1
(, )

(0, 0)( , 0)

Figure 9: Selected components of the static susceptibility χ(q ) throughout the Brillouin zone, for the two-band Hubbard model described in the text, for β = 60.
phase such that only 〈Rφ x 〉 is non-zero. In this set-up Rφ y generates a spin rotation of theorder parameter (Goldstone mode) and I φ y couples to Sz (only in the ordered phase).
In Fig. (9) we show the static susceptibilities on a fine grid in the 2D Brillouin zone. Unlike in the strong coupling case of Ref. [33], the spin susceptibility is dominated by a Fermi
surface nesting present both in the ordered and disordered (not shown here) phases, which
gives rise to a peak at an incommensurate vector on the X -M zone boundary. The response
corresponding to the Goldstone mode in the middle panel exhibits the expected divergence atthe ordering wave vector (q = 0). The excitonic susceptibility in the right panel exhibits, inaddition to the main (finite) peak at q = 0, additional peaks that coincide with the maxima of
the spin susceptibility. This reflects the coupling between Sz and I φ y induced by the symmetry breaking. While in the strong coupling regime of Ref. [33] Sz was a slave to the dynamics
of I φ y , here we can see that Sz affects I φ y . In Fig. (10) we show the absorptive (imaginary)parts of the dynamical susceptibilities obtained by the analytic continuation described in Supplemental Material of Ref. [33]. It reveals the Goldstone nature of the Rφ y response andcomplex nature of the spin response. Interestingly, the sharp response in the vicinity of Γ does
not reflect formation of a bound state, but is a consequence of parallel bands upon opening of
the excitonic gap. The low-energy hot spot on the X -M linearand its counterpart in the static
susceptibility reflect the vicinity of an antiferromagnetic phase [39].
The presented susceptibilities obtained from the IR inputs are in excellent agreement withbenchmark data obtained using the Legendre representation used in Ref. [33]. The current
setup is more flexible, insofar as it is not limited by the sampling window, neither in the bosonicnor the fermionic Matsubara frequency domain. It is also very compact. The sparse grid solely
saves computational time and memory footprints for QMC significantly. The tensor regression
further compresses the sparsely sampled data by several orders of magnitude: The measured
data is 700 MB large on the sparse grid, while the tensor network representation takes up only
330 kB for D = 60.
In general, if a too large D is employed, one could overfit QMC noise, giving a rise tooscillatory behavior between the sampling points in the interpolated data. A practical recipe
for avoiding overfitting is to use the value of D that minimizes test errors rather than training
errors. In the present study, however, we did not observe overfitting behavior. This may be
because the fitting parameters is still much smaller than the fitted QMC data in size. Moredetailed analysis of the stability of the fitting procedure is a topic of future studies.

14

SciPost Phys. 8, 012 (2020)

Sz

1.0

y

y

0.8
0.6
0.4
0.2
0.0

X

M

X

M

X

M

Figure 10: Selected modes of the dynamical susceptibilityy for the two-band Hubbard
model described in the text, in the presence of the excitonic condensate, β = 60.

7

Summary

Based on the IR basis, we have introduced a procedure for generating sparse grids in the
Matsubara frequency domain and a fitting algorithm based on a tensor network representation.
These two enable an efficient transformation of numerical data from Matsubara to IR domain.
The tensor network representation provides a model-independent way to compress the IR
expansion coefficients (IR tensor) by decoupling the frequency and spins/orbital dependence.
Low-temperature calculations for multi-orbital systems benefit from this compression.
We have demonstrated the efficiency and accuracy of the present method in DMFT calculations: static susceptibility calculations for single-band Hubbard model and dynamic susceptibility calculations for two-band Hubbard model with low symmetry. We have shown that
accurate susceptibilities can be obtained already with low-rank approximation of the IR tensor.
The sparse sampling and the tensor network decomposition are independent proceduresthat are controlled separately. The size of sparse sampling grid, and thus its computationalcost, depends only on temperature, the energy window and the desired accuracy. The “compression rate” and the accuracy of tensor network representation are controlled by the rankof decomposition D. In the present work, we have demonstrated that the local 2P Green’s
function can be compressed from 700 MB to 330 kB for the two-band Hubbard model. The
concept of tensor network representation is flexible and further compression may be possiblefor different tensor network topologies. The choice of ideal tensor network topology requires
an extensive experience with the performance of the method for various models and is beyond
the scope of the present work.
Potential applications of the present scheme include DFT+DMFT calculations for realistic
multi-orbital models and diagrammatic extensions of DMFT. It is highly desirable to develop
efficient methods for solving equations at the 2P level such as Bethe-Salpeter and parquet
equations directly in the tensor network format with the sparse sampling. This requires efficient evaluation of contractions of 2P quantities, e.g., a vertex function and a generalized
susceptibility. Potentially useful techniques for manipulating matrix product states and tensor
networks have already been developed in other fields of condensed matter theory [26,27,40].

Acknowledgments
HS thanks Lei Wang for stimulating discussions on tensor networks and machine learning.
Part of the calculations were run on the facilities of the Supercomputer Center at the Institutefor Solid State Physics, University of Tokyo, using codes based on ALPSCore [41, 42]. We
used the irbasis library [43] for computing IR basis functions. We used DCore [44] based
15

SciPost Phys. 8, 012 (2020)

on TRIQS [45] and TRIQS/DFTTools [46] for DMFT calculations of the single-band Hubbard
model.
Funding information H.S, J.O and K.Y were supported by JSPS KAKENHI Grant
No. 18H01158. H.S. was supported by JSPS KAKENHI Grant No. 16K17735. J.O. was supported by JSPS KAKENHI Grant No. 18H04301 (J-Physics). K.Y. was supported by Building ofConsortia for the Development of Human Resources in Science and Technology, MEXT, Japan.
D.G and J.K were supported by the ERC Grant Agreements No. 646807 under EU Horizon
2020. D.G was supported by the Czech Science Foundation (GAČR) under Project No. GA1916937S. This work was supported by The Ministry of Education, Youth and Sports from theLarge Infrastructures for Research, Experimental Development and Innovations project"IT4Innovations National Supercomputing Center – LM2015070". MW and EG were supportedby the Simons Foundation via the Simons Collaboration on the many-electron problem. Access to computing and storage facilities owned by parties and projects contributing to the
National Grid Infrastructure MetaCentrum provided under the programme “Projects of Large
Research, Development, and Innovations Infrastructures” (CESNET LM2015042), and by theAustrian Federal Ministry of Science, Research and Economy through the Vienna Scientific
Cluster (VSC) Research Center, is greatly appreciated.

A

Intermediate representation at fixed bosonic frequency

The frequency dependence of G(iωn , iωn0 , iνm ) can be decomposed into 16 distinct components shown in Table. 1. For instance, the first component (r = 1) depends on the three frequencies through iνm + iωn , −iωn , iωn0 , which defines the structure of discontinuity planes in
the imaginary-time domain. To be more specific, as discussed in Ref. [22], the first component
can be discontinuous at three equal-time planes: τ1 = τ4 , τ2 = τ4 and τ3 = τ4 (mod β).
Such a function with this discontinuity structure may be well approximated byρ (1) (ε1 , ε2 , ε3 ),(iνm + iωn − ε1 )(−iωn − ε2 )(iωn0 − ε3 )(14)

where ρ (1) (ε1 , ε2 , ε3 ) is an auxiliary spectrum bounded in [−ωmax , ωmax ]. Applying the sameprocedure to all the 16 components, we obtain the assumption that
16 Z ωmaxXρ (r) (ε1 , ε2 , ε3 )G(iωn , iωn0 , iνm ) 'dε1 dε2 dε3.
(15)0 − ε )(iω00 − ε )
(iω
−
ε
)(iω
1
2
3
−ω
r=1
max

For instance, the term for r = 1 can be decomposed as1(i(νm + ωn ) − ε1 )(−iωn − ε2 )(iωn0 − ε3 )


1
1
1
1
=
+
iνm − ε1 − ε2 i(νm + ωn ) − ε1 −iωn − ε2 (iωn0 − ε3 ) F
1=K (iωn + iνm , ε1 ) − K F (iωn + ε2 ) K F (iωn0 , ε3 ).
iνm − ε1 − ε2(16)

One can see the first and second terms in the last line are the products of two fermionic kernels.
The first term can be represented compactly as
X SlF SlF VlF (ε1 )VlF (ε3 )1
2
1
2
UlF (iωn + iνm )UlF (iωn0 ),
(17)
1
2
iν
−
ε
−
ε
m
1
2
l ,l
1 2

16

SciPost Phys. 8, 012 (2020)

where the coefficient in the parenthesis decays as fast as the singular values with respect to l1
and l2 .
Applying the same procedure to all the terms in Eq. (15), one obtains a compact overcomplete representation
∞ ¦
X
(1)
Gss0 l l UslF (iωn )UsF0 l (iωn0 )X

G(iωn , iωn0 , iωm ) =1 2

s,s0 =0,1 l1 ,l2 =0

1

2

(2)

+ Gss0 l l UlB s (iωn + (−1)s+1 iωn0 )UlF s0 (iωn0 )
1
2
1 2
©
(3)
B
s+1
F
+ Gss0 l l Ul s (iωn0 + (−1) iωn )Ul s0 (iωn ) .
1 2

1

2

(18)

Here we definedUslα (iωn ) ≡ Ulα (iωn + siωm )(19)

for s = 0, 1.
We can recast Eq. (18) into a form more analogous to Eq. (8) asG(iωn , iωn0 , iωm ) =12 X
X
r=1 l1 ,l2

B

(r)

(r)

1

2

G (r, l1 , l2 ; iν)Ul (iω)Ul (iω0 ).

(20)

Optimization algorithm for tensor regression

We minimize the cost function in Eq. (13) by means of a accelerated alternating least squares(ALS) method. The essential idea of ALS is to optimize each tensor in {x (0) , x (1) , · · · } at one
time. The optimization of a single tensor reduces to a convex optimization problem. In ALS,
we sweep through all the tensors until the value of the cost function is converged. In addition,
we introduce recently proposed acceleration techniques to improve the convergence of ALS.
In the following, we detail the procedure of ALS and the acceleration techniques.

B.1

Alternating least squares

We explain how to optimize the tensors in Eq. (13) by alternating least squares.
B.1.1

Optimization of x (0)The minimization of Eq. (13) with respect to x (0) can be recast intoX(0)min || y(W, O) −A(W, O; d, r)x d r ||2 + α||x (0) ||2 ,
x (0)

(21)

r,d

which is a regularized convex optimization in the well-known form of Ridge regression. The
tensor A reads(3)

A(W, O; d, r) ≡ B(W, r, d)x dO ,

(22)

where
B(W, r, d) ≡

X

(r)

(1)

1

1

Ul (iω)x dl

l1

X
l2

17

(r)

(2)

2

2

Ul (iω0 )x dl .

(23)

SciPost Phys. 8, 012 (2020)

In matrix form, Eq. (21) reads
min ||y − Ax (0) ||2 + α||x (0) ||2 ,

(24)

x (0)

where y and x (3) are flatted 1D arrays. The matrix is size of (NW NO , Nr D), where Nr is thenumber of different representations (=12), NO is the size of the combined index for spin andorbitals.
In practice, we solve Eq. (24) by an iterative method, LSQR [47], without constructing thematrix A explicitly. In LSQR, the matrix A is used only to compute Av and A† u for various v
and u. Hence we store the (precomputed) tensor B in memory, and compute these products bymeans of tensor contractions. We illustrate the tensor contractions for computing Av and A† u
in Figs. 11(a) and 11(b), respectively. This approach not only reduces memory footprints but
also reduces the computational complexity from O(NW DNr NO ) to O(NW DNr ) + O(NW DNO ).
B.1.2

Optimization of x (1)The optimization of x (1) can be done in a way very similar to that of x (0) . Thus, we focus only
on the differences. The reduced least squares problem reads
X(1)
min || y(W, O) −A0 (W, O; d, l1 )x dl ||2 + α||x (1) ||2 ,
(25)
x (1)

1

d,l1

where
(1)

A0 (W, O; d, l1 ) ≡ B 0 (W, d, l1 )x dO ,

(26)

where
(0)

B 0 (W, d, l1 ) ≡ x d r

X

(r)

(1)

1

1

Ul (iω)x dl .

(27)

l1
†

As illustrated in Fig. 11, we can readily compute A0 v and A0 u by tensor contractions.
B.1.3

Optimization of x (2)The tensor x (2) can be optimized in exactly the same way as x (1) . Thus, we do not describe
the optimization of x (2) for simplicity.
B.1.4

Optimization of x (3)The optimization of x (3) is rather simple. The reduced least squares problem reads
X(3)
min || y(W, O) −A00 (W ; d, O)x dO ||2 + α||x (3) ||2 ,
x (3)

(28)

d,O

where the tensor A00 can be stored in memory. Furthermore, this optimization problem isseparable with respect to D and thus can be solved independently.

B.2

Acceleration techniques and convergence condition

In the previous subsection, we have explained how to perform one sweep through the tensors.
This single ALS sweep corresponds to the function ALS in Algorithm 1. The function ALS takes
an array obtained by flattening tensors of fitting parameters as input. After a single sweep, theupdated tensors are returned as a flattened array. Here, flattening means recasting multiple
18

SciPost Phys. 8, 012 (2020)

tensors of complex numbers into a single one-dimensional array of real numbers (the order isarbitrary).
The whole procedure of the accelerated ALS is illustrated in Algorithm 1. The main difference from the plain ALS is the existence of β. In the second last line, β(x k − x k−1 ) acts asa momentum term for β > 0. Although this momentum term accelerates the convergence by
updating the parameters aggressively, this sometimes leads to oscillatory behavior or divergence. We use a restarting mechanism to stabilize the accelerated ALS. In practice, when therestarting condition f (x k ) > f (x k−1 ) is met ( f is the cost function), a ALS sweep is forced
by setting β = 0 (see the comment in Algorithm 1). For more details on the acceleration
techniques, please refer to Ref. [48].
The loop is exited when a convergence condition is met. δ in Algorithm 1) is a relative
tolerance.

B.3

Technical details and numerical results

We parallelize the whole fitting procedure by MPI with respect to frequencies. This parallelization is efficient particularly for G. We parallelize the LSQR implementation in the SciPy Python
packageusing MPI with respect to sampling frequencies.
Figure 12 shows the convergence of the root squared errors of the local susceptibility forthe single-band Hubbard model analyzed in Sec. 5. One can see that the fitting errors quickly
converge. The small oscillatory behavior is due to the acceleration.
(a)

W

W

O

B or B 0

O
d

x (3)

O(NW DNr/l )r or l1x (3)uO(NW NO )

W

d

O

(b)
v

W

O

B or Br or l1

O(NW NO )0Wd

B or B 0
r or l1x (3)dO(NW Nr/l )r or l1dFigure 11: Graphical representations of tensor contractions for solving Eqs. (21)and (25), iteratively. The panel (a) illustrates the order of tensor contractions forcomputing Au and A0 u. The panel (b) illustrates the procedure for computing A† v
and (A0 )† v . The computational complexity for each operation is shown.

19

SciPost Phys. 8, 012 (2020)

Algorithm 1 Accelerated alternating least squares
function ALS(x ){x (i) } ← UNFLATTEN(x )for i = 0, 1, 2, 3 do
Construct A(i)x (i) ← LSQR( y, A(i) , α)
end for
Return FLATTEN({x (i) })end function. Unflatten input array. Flatten updated tensors

Root squared error

x 0 ← Random initial values
x 1 ← ALS(x 0 )for k = 1, 2, ... do
if f (x k ) > f (x k−1 ) thenx k ← x k−1
β ←0
else
β ←1
end ifx k+1 ← ALS(x k + β(x k − x k−1 ))if | f (x k+1 ) − f (x k )| < δ × f (x k ) thenExit loop
end ifend for. Force ALS at this iteration

100

D=5
D=10
D=15
D=30

10-1
10-2
10-3
10-4
10-5

0

50

100
150
Iterations

200

250

Figure 12: Convergence of root squared errors for the local susceptibility X loc of thesingle-band Hubbard model analyzed in Section 5.

References[1] G. D. Mahan, Many-Particle Physics, Kluwer Academic/Plenum Publishers, New York(2000).
[2] A. Georges, G. Kotliar, W. Krauth and M. J. Rozenberg, Dynamical mean-field theory ofstrongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys.
68, 13 (1996), doi:10.1103/RevModPhys.68.13.
20

SciPost Phys. 8, 012 (2020)[3] A. Toschi, A. A. Katanin and K. Held, Dynamical vertex approximation: A
step beyond dynamical mean-field theory, Phys. Rev. B 75, 045118 (2007),
doi:10.1103/PhysRevB.75.045118.
[4] A. Valli, T. Schäfer, P. Thunström, G. Rohringer, S. Andergassen, G. Sangiovanni, K. Held and A. Toschi, Dynamical vertex approximation in its parquet implementation: Application to Hubbard nanorings, Phys. Rev. B 91, 115115 (2015),
doi:10.1103/PhysRevB.91.115115.
[5] A. Galler, P. Thunström, P. Gunacker, J. M. Tomczak and K. Held, Ab initio dynamical vertex
approximation, Phys. Rev. B 95, 115107 (2017), doi:10.1103/PhysRevB.95.115107.
[6] A. N. Rubtsov, M. I. Katsnelson, A. I. Lichtenstein and A. Georges, Dual fermion approachto the two-dimensional Hubbard model: Antiferromagnetic fluctuations and Fermi arcs,
Phys. Rev. B 79, 045133 (2009), doi:10.1103/PhysRevB.79.045133.
[7] H. Hafermann, G. Li, A. N. Rubtsov, M. I. Katsnelson, A. I. Lichtenstein and H. Monien,
Efficient perturbation theory for quantum lattice models, Phys. Rev. Lett. 102, 206401(2009), doi:10.1103/PhysRevLett.102.206401.
[8] A. E. Antipov, E. Gull and S. Kirchner, Critical exponents of strongly correlated fermion
systems from diagrammatic multiscale methods, Phys. Rev. Lett. 112, 226401 (2014),
doi:10.1103/PhysRevLett.112.226401.
[9] J. Otsuki, H. Hafermann and A. I. Lichtenstein, Superconductivity, antiferromagnetism,
and phase separation in the two-dimensional Hubbard model: A dual-fermion approach,Phys. Rev. B 90, 235132 (2014), doi:10.1103/PhysRevB.90.235132.
[10] H. Kusunose, Influence of spatial correlations in strongly correlated electron systems: Extension to dynamical mean field approximation, J. Phys. Soc. Jpn. 75, 054713 (2006),
doi:10.1143/JPSJ.75.054713.
[11] G. Rohringer, A. Toschi, H. Hafermann, K. Held, V. I. Anisimov and A. A.
Katanin, One-particle irreducible functional approach: A route to diagrammatic extensions of the dynamical mean-field theory, Phys. Rev. B 88, 115112 (2013),
doi:10.1103/PhysRevB.88.115112.
[12] C. Taranto, S. Andergassen, J. Bauer, K. Held, A. Katanin, W. Metzner, G. Rohringer andA. Toschi, From infinite to two dimensions through the functional renormalization group,
Phys. Rev. Lett. 112, 196402 (2014), doi:10.1103/PhysRevLett.112.196402.
[13] M. Kitatani, N. Tsuji and H. Aoki, FLEX+DMFT approach to thed-wave superconducting
phase diagram of the two-dimensional Hubbard model, Phys. Rev. B 92, 085104 (2015),
doi:10.1103/PhysRevB.92.085104.
[14] T. Ayral and O. Parcollet, Mott physics and spin fluctuations: A unified framework, Phys.
Rev. B 92, 115109 (2015), doi:10.1103/PhysRevB.92.115109.
[15] T. Ayral and O. Parcollet, Mott physics and collective modes: An atomic approximation of the four-particle irreducible functional, Phys. Rev. B 94, 075159 (2016),
doi:10.1103/PhysRevB.94.075159.
[16] G. Li, Hidden physics in the dual-fermion approach: A special case of a nonlocal expansion
scheme, Phys. Rev. B 91, 165134 (2015), doi:10.1103/PhysRevB.91.165134.

21

SciPost Phys. 8, 012 (2020)[17] J. Kuneš, Efficient treatment of two-particle vertices in dynamical mean-field theory, Phys.
Rev. B 83, 085102 (2011), doi:10.1103/PhysRevB.83.085102.
[18] G. Rohringer, A. Valli and A. Toschi, Local electronic correlation at the two-particle level,
Phys. Rev. B 86, 125114 (2012), doi:10.1103/PhysRevB.86.125114.
[19] G. Li, N. Wentzell, P. Pudleiner, P. Thunström and K. Held, Efficient implementation of theparquet equations: Role of the reducible vertex function and its kernel approximation, Phys.
Rev. B 93, 165103 (2016), doi:10.1103/PhysRevB.93.165103.
[20] J. Kaufmann, P. Gunacker and K. Held, Continuous-time quantum Monte Carlocalculation of multiorbital vertex asymptotics, Phys. Rev. B 96, 035114 (2017),
doi:10.1103/PhysRevB.96.035114.
[21] N. Wentzell, G Li, A. Tagliavini, C. Taranto, G. Rohringer, K. Held, A. Toschi and S. Andergassen, High-frequency asymptotics of the vertex function: diagrammatic parametrization
and algorithmic implementation (2016), arXiv:1610.06520.
[22] H. Shinaoka, J. Otsuki, K. Haule, M. Wallerberger, E. Gull, K. Yoshimi and M. Ohzeki,
Overcomplete compact representation of two-particle Green’s functions, Phys. Rev. B 97,
205111 (2018), doi:10.1103/PhysRevB.97.205111.
[23] J. Li, M. Wallerberger, C.-N. Yeh, N. Chikano, E. Gull and H. Shinaoka, Sparse sampling
approach to efficient ab initio calculations at finite temperature, Phys. Rev. B 101, 035144(2020), doi:10.1103/PhysRevB.101.035144.
[24] H. Shinaoka, J. Otsuki, M. Ohzeki and K. Yoshimi, Compressing Green’s function using intermediate representation between imaginary-time and real-frequency domains, Phys. Rev.
B 96, 035147 (2017), doi:10.1103/PhysRevB.96.035147.
[25] J. Otsuki, M. Ohzeki, H. Shinaoka and K. Yoshimi, Sparse modeling in quantum manybody problems, J. Phys. Soc. Jpn. 89, 012001 (2020), doi:10.7566/JPSJ.89.012001.
[26] U. Schollwöck, The density-matrix renormalization group in the age of matrix productstates, Ann. Phys. 326, 96 (2011), doi:10.1016/j.aop.2010.09.012.
[27] J. C. Bridgeman and C. T. Chubb, Hand-waving and interpretive dance: an introductory course on tensor networks, J. Phys. A: Math. Theor. 50, 223001 (2017),
doi:10.1088/1751-8121/aa6dc3.
[28] https://www.tensors.net.
[29] F. Hummel, T. Tsatsoulis and A. Grüneis, Low rank factorization of the Coulomb
integrals for periodic coupled cluster theory, J. Chem. Phys. 146, 124105 (2017),
doi:10.1063/1.4977994.
[30] M. Motta, J. Shee, S. Zhang and G. Kin-Lic Chan, Efficient ab initio auxiliary-field quantum
Monte Carlo calculations in Gaussian bases via low-rank tensor decomposition, J. Chem.
Theory Comput. 15, 3510 (2019), doi:10.1021/acs.jctc.8b00996.
[31] V. Khoromskaia, B. N. Khoromskij and R. Schneider, Tensor-structured factorized calculation of two-electron integrals in a general basis, SIAM J. Sci. Comput. 35, A987 (2013),
doi:10.1137/120884067.
[32] J. Otsuki, K. Yoshimi, H. Shinaoka and Y. Nomura, Strong-coupling formula formomentum-dependent susceptibilities in dynamical mean-field theory, Phys. Rev. B 99,
165134 (2019), doi:10.1103/PhysRevB.99.165134.
22

SciPost Phys. 8, 012 (2020)[33] D. Geffroy, J. Kaufmann, A. Hariki, P. Gunacker, A. Hausoel and J. Kuneš, Collective modes
in excitonic magnets: Dynamical mean-field study, Phys. Rev. Lett. 122, 127601 (2019),
doi:10.1103/PhysRevLett.122.127601.
[34] J. Kuneš, Phase diagram of exciton condensate in doped two-band Hubbard model, Phys.
Rev. B 90, 235140 (2014), doi:10.1103/PhysRevB.90.235140.
[35] H. Shinaoka, E. Gull and P. Werner, Continuous-time hybridization expansion quantum
impurity solver for multi-orbital systems with complex hybridizations, Comput. Phys. Commun. 215, 128 (2017), doi:10.1016/j.cpc.2017.01.003.
[36] H. Shinaoka, F. Assaad, N. Blümer and P. Werner, Quantum Monte Carlo impurity solvers
for multi-orbital problems and frequency-dependent interactions, Eur. Phys. J. Spec. Top.
226, 2499 (2017), doi:10.1140/epjst/e2017-70050-x.
[37] P. Werner and A. J. Millis, Hybridization expansion impurity solver: General formulation
and application to Kondo lattice and two-orbital models, Phys. Rev. B 74, 155107 (2006),
doi:10.1103/PhysRevB.74.155107.
[38] P. Werner, A. Comanac, L. de’ Medici, M. Troyer and A. J. Millis, Continuoustime solver for quantum impurity models, Phys. Rev. Lett. 97, 076405 (2006),
doi:10.1103/PhysRevLett.97.076405.
[39] S. Hoshino and P. Werner, Superconductivity from emerging magnetic moments, Phys. Rev.
Lett. 115, 247001 (2015), doi:10.1103/PhysRevLett.115.247001.
[40] H.-J. Liao, J.-G. Liu, L. Wang and T. Xiang, Differentiable programming tensor networks,
Phys. Rev. X 9, 031041 (2019), doi:10.1103/PhysRevX.9.031041.
[41] A. Gaenko et al., Updated core libraries of the ALPS project, Comput. Phys. Commun. 213,
235 (2017), doi:10.1016/j.cpc.2016.12.009.
[42] M. Wallerberger et al.,arXiv:1811.08331.

Updated core libraries of the ALPS project (2018),[43] N. Chikano, K. Yoshimi, J. Otsuki and H. Shinaoka, irbasis: Open-source database andsoftware for intermediate-representation basis functions of imaginary-time Green’s function,
Comput. Phys. Commun. 240, 181 (2019), doi:10.1016/j.cpc.2019.02.006.
[44] https://github.com/issp-center-dev/DCore.
[45] O. Parcollet, M. Ferrero, T. Ayral, H. Hafermann, I. Krivenko, L. Messio and P. Seth, TRIQS:
A toolbox for research on interacting quantum systems, Comput. Phys. Commun. 196, 398(2015), doi:10.1016/j.cpc.2015.04.023.
[46] M. Aichhorn et al., TRIQS/DFTTools: A TRIQS application for ab initio calculations of correlated materials, Comput. Phys. Commun. 204, 200 (2016),
doi:10.1016/j.cpc.2016.03.014.
[47] C. C. Paige and M. A. Saunders, LSQR: An algorithm for sparse linear equations and sparse
least squares, ACM Trans. Math. Softw. 8, 43 (1982), doi:10.1145/355984.355989.
[48] D. Mitchell, N. Ye and H. De Sterck, Nesterov acceleration of alternating least squares
for canonical tensor decomposition: Momentum step size selection and restart mechanisms(2018), arXiv: 1810.05846.

23