SciPost Phys. Lect.Notes 63 (2022)

Efficient ab initio many-body calculations
based on sparse modeling of Matsubara Green’s function
Hiroshi Shinaoka1? , Naoya Chikano1 , Emanuel Gull2 , Jia Li2 , Takuya Nomoto7 ,Junya Otsuki4 , Markus Wallerberger5 , Tianchun Wang3 and Kazuyoshi Yoshimi6
1 Department of Physics, Saitama University, Saitama 338-8570, Japan
2 University of Michigan, Ann Arbor, Michigan 48109, USA
3 Department of Applied Physics, University of Tokyo, Japan
4 Research Institute for Interdisciplinary Science, Okayama University,
Okayama 700-8530, Japan
5 Institute of Solid State Physics, TU Wien, 1040 Vienna, Austria
6 Institute for Solid State Physics, University of Tokyo, Chiba 277-8581, Japan
7 Research Center for Advanced Science and Technology, University of Tokyo,
Tokyo 153-8904, Japan? h.shinaoka@gmail.com

Abstract
This lecture note reviews recently proposed sparse-modeling approaches for efficient
ab initio many-body calculations based on the data compression of Green’s functions.
The sparse-modeling techniques are based on a compact orthogonal basis, an intermediate representation (IR) basis, for imaginary-time and Matsubara Green’s functions. A
sparse sampling method based on the IR basis enables solving diagrammatic equations
efficiently. We describe the basic properties of the IR basis, the sparse sampling method
and its applications to ab initio calculations based on the GW approximation and theMigdal–Eliashberg theory. We also describe a numerical library for the IR basis and thesparse sampling method, sparse-ir, and provide its sample codes. This lecture notefollows the Japanese review article [H. Shinaoka et al., Solid State Physics 56(6), 301(2021)].
Copyright H. Shinaoka et al.
This work is licensed under the Creative Commons
Attribution 4.0 International License.
Published by the SciPost Foundation.

Received 25-06-2021
Accepted 20-07-2022
Check forupdatesPublished 22-09-2022
doi:10.21468/SciPostPhysLectNotes.63

Contents
1 Introduction

22 Intermediate-representation (IR) of Matsubara Green’s functions

3

3 Sparse sampling
3.1 Dyson equation
3.2 Sparse sampling and Fourier transform
3.3 Notes on calculations using the sparse sampling method

6
6
7
10

1

SciPost Phys. Lect.Notes 63 (2022)

3.4
3.5

Related approaches
Discrete Lehmann representation (DLR)

10
11

4 Review of applications
4.1 Noble-gas atoms, silicon crystals: Self-consistent GW and second-order perturbation (GF2) calculations
4.2 First-principles calculations of superconducting transition temperature
based on Migdal–Eliashberg equation
4.3 Calculation of magnetic interaction based on Liechtenstein formula
4.4 Other applications

12

5 sparse-ir library
5.1 Installation
5.2 Usage
5.3 Example: Second-order perturbation theory

16
17
17
19

6 Summary and future directions

21

A Some properties of the IR basis functions

22

B Numerical stability of sparse sampling

23

References

25

1

12
14
15
16

Introduction

Perturbation and quantum field theories based on Green’s functions are widely used for abinitio and quantum many-body calculations. The imaginary-time formalism based on Matsubara Green’s functions is well known for its simplicity in numerical treatment. It has beenused for ab initio calculations and model calculations based on various theories such as theGW approximation [1], the random phase approximation (RPA) [2], the fluctuation exchange
approximation (FLEX) [3, 4], the dynamical mean-field theory (DMFT) [5], and quantum
Monte Carlo methods [6]. Most readers of this Lecture Note may have had experience using such computational methods. Applications of imaginary-time formalism are not limited tocondensed matter physics but is also used in quantum chemistry and high energy physics [7–9].
Although the Matsubara Green’s function is used in a wide range of fields, its efficient numerical handling has yet to be fully established. For example, in ab initio calculations, it isnecessary to simultaneously deal with multiple energy scales that differ by orders of magnitude. Typical energy scales are the width of low-energy bands (several eV to several tens of eV)and the low temperature at which physical phenomena occur (1 K ' 0.1 meV). In such cases,
the number of imaginary (Matsubara) frequencies required in the calculations increases, andthe computation time and the amount of memory required increase to an unmanageable level.
Furthermore, calculations at the two-particle level (e.g., susceptibility calculations using
the Bethe-Salpeter equation) require the handling of vertex functions that depend on multiple
imaginary frequencies. Therefore, calculations with full frequency dependence are extremely
costly even for simple model calculations, and applications to realistic materials are impractical. In this Lecture Note, we review efficient many-body and ab initio calculation methods
based on the data compression of Matsubara Green’s functions. A compact basis for Matsubara
2

SciPost Phys. Lect.Notes 63 (2022)

Green’s function named intermediate representation (IR) was proposed in 2017 [10, 11]. In
2020, an efficient calculation method based on sparse sampling was developed [12]. Since
then, its applications to ab initio calculations have been rapidly spreading [13–18].
This Lecture Note summarizes the theoretical progress made during the last few years.
A comprehensive overview is given of 1) the basic properties of the IR basis, 2) the sparse
sampling method for fast computation of diagrammatic equations, and 3) numerical librarysparse-ir. We hope that this Lecture Note help readers start using the IR basis and thesparse sampling method in various fields.

2Intermediate-representation (IR) of Matsubara Green’s functions

The intermediate representation (IR) is a basis set in which the imaginary frequency and timedependence of the two-point imaginary-time/-frequency correlation function can be expanded
compactly [10, 19]. Let us summarize its definition and properties.
The general definition of the imaginary-time two-point correlation function is:G(τ) = −〈Tτ Aα (τ)B α (0)〉 ,(1)

where β = 1/T is the inverse temperature, τ is the imaginary time, Tτ is the time ordering
operator, and α is either F (fermion) or B (boson). Here, we take kB = 1. In addition, A and
B are the operators. In the important special case of a one-particle Green’s function, A is anannihilation operator, and B is a creation operator.
The spectral representation of the two-point correlation function is given as follows:Z ωmaxdω0 K α (τ, ω0 )A(ω0 ) ,G(τ) = −

(2)

−ωmax

where A(ω) is a corresponding spectral function and the kernel K α (τ, ω) is defined as
 −τω

 e
(α = F) ,αe−βωK (τ, ω) = 1 +−τωe(α = B) ,1 − e−βω(3)

for 0 < τ < β. We take ωmax to be sufficiently large such that the spectral function A(ω) isnonzero only for ω ∈ [−ωmax , ωmax ]. We can regard Eq. (2) as an integral equation which connects the spectral function A(ω) to the correlation function G(τ) through the kernel K α (τ, ω).
To avoid the divergence of the bosonic kernel at ω = 0, we reformulate Eq. (2) as [20]Z ωmaxdω0 K(τ, ω0 )ρ(ω0 ) ,G(τ) = −

(4)

−ωmax

where we call K(τ, ω)[= K F (τ, ω)] a logistic kernel, and ρ(ω) is the modified spectral function
defined byA(ω)(α = F) ,(5)
ρ(ω) ≡A(ω)(α = B) .
tanh(βω/2)This allows to use the same kernel for fermions and bosons. Note that a different regularization
was used in the original proposal of IR [10].
3SciPost Phys. Lect.Notes 63 (2022)(b) 0.4

β = 100
β = 1000
β = 10000

Ul (τ )

0.2

Sl

10−4
10−8
10−12

(c)4

l=0
l=1
l = 10Vl (ω)

(a) 1

0.0

−0.2

0

25

50
l

75

0.0

100

0.2

0.4

0.6

0.8

1.0

τ /β

l=0
l=1
l = 10

2
0

−2
−1.0

−0.5

0.0
ω

0.5

1.0

Figure 1: (a) Singular value Sl computed for various values of β. (b), (c) IR basis
functions Ul (τ) and Vl (ω) computed for l = 0, 1, 10 and β = 100, respectively. ωmax
is fixed at ωmax = 1.
The task at hand is to find a compact representation of the imaginary-frequency (Matsubara) Green’s function G(iω) that retains full information of the real-frequency spectral function
ρ(ω). For this, we perform a singular value expansion of the logistic kernel K for a given βand ωmax [10]:∞X
K(τ, ω0 ) =Ul (τ)Sl Vl (ω0 ) .
(6)
l=0
0

The Fourier transform of K(τ, ω ) with respect to τ yields the singular value expansion in theimaginary-frequency domainαK̂ (iω

(α)

0,ω ) ≡ −Zβ0dτK(τ, ω )e0iω(α) τ=−∞X

Ûlα (iω(α) )Sl Vl (ω0 ) ,

(7)

l=0

where
Ûlα (iω(α) ) ≡ZβdτUl (τ)eiω

(α)

τ

.

(8)

0

Here, iω(α) is the imaginary frequency corresponding to the statistics α and we attached a
hatˆto the quantities defined in the imaginary-frequency domain. The minus sign in Eq. (7)originates from the convention K(τ, ω) > 0. In the following, we will denote iω(α) as iω forsimplicity.
The singular value expansion has similar properties as the singular value decompositionof matrices: The singular values Sl satisfy S0 > S1 > . . . > 0. The left singular functions{U0 (τ), U1 (τ), . . .} form an orthonormal set on the imaginary time axis τ ∈ [0, β], while theright singular functions {V0 (ω), V1 (ω), . . .} form an orthonormal set on the real frequency axisω ∈ [−ωmax , ωmax ].
The so-called IR basis functions are nothing but the left singular functions Ul (τ) [Ûlα (iω)]
and the right singular functions Vl (ω). It should be emphasized that the IR basis functions only
depend on the inverse temperature β, statistics α, and the energy (frequency) cutoff ωmax anddo not depend on the details of the system.
Figure 1 shows the singular values Sl and the IR basis functions Ul (τ) and Vl (ω). The IR
basis functions have the following interesting properties:
Property 1 The singular values Sl are non-degenerate, non-negative, and monotonically decrease faster than exponentially with respect to l 1 .
1

For numerical evidence, refer to Appendix D of [21].

4SciPost Phys. Lect.Notes 63 (2022)

Property 2 The number of numerically significant singular values (e.g., Sl /S0 ≥ 10−15 ) is
determined by the dimensionless quantity Λ ≡ βωmax and increases only logarithmicallywith respect to Λ [21].
Property 3 Ul (τ) and Vl (ω) can be chosen to be real functions, and they then become even(odd) functions for even (odd) l. In this convention, Ûlα (iω) is purely imaginary or real.
Property 4 Ul (τ) and Vl (ω) have l roots. In addition, in the limit of βωmax → 0 (the hightemperature limit), Ul (τ(x)) and Vl (ω(x)) coincide with the Legendre polynomial Pl (x)up to a constant [τ(x) ≡ β(x + 1)/2, ω(x) ≡ xωmax for −1 < x < 1] 2 .
We sketch mathematical proofs of Properties 1, 3, and 4 in Appendix A. Property 2 has beenverified only numerically3 .
Now, let us expand an “arbitrary” two-point correlation using a finite number (L) of the IR
basis functions:
L−1
X
Ĝ(iω) =Ûlα (iω)Gl + ε L ,(9)l=0where ε L is a truncation error. If the corresponding spectral function ρ(ω) is finite only in theinterval [−ωmax , ωmax ], the expansion coefficient Gl can be evaluated using the orthonormalityof the IR functions as
Z ωmaxGl = −Sldωρ(ω)Vl (ω) ≡ −Sl ρl .

(10)

−ωmax

From Eq. (10), it can be seen that Gl vanishes at least as fast as Sl , i.e., faster than the exponential function from Property 1.
Figure 2 shows the convergence of the expansion coefficients ρl and Gl for some simple
spectral functions ρ(ω) 4 The convergence of ρl strongly depends on the shape of ρ(ω).
In general, ρl decays exponentially in the case of a smooth spectral function, but does notdecay at all in the case of a δ function spectrum. By contrast, Gl always decays as fast as oreven faster than the singular values, regardless of the spectral functions. We would like toemphasize again that this convergence is determined by the singular values of the kernel anddoes not depend on the details of the system.
In Eq. (10), we can see that when we transform ρl into Gl , the contribution of ρl is suppressed by the singular values exponentially at a large l. This means that some of the information in the spectral function is lost in the transformation to the imaginary-time domain.
Thus, it is usually difficult to recover the exact spectral function from the numerical data of
the Matsubara Green’s function. This is because, in the inverse transformation ρl = −(Sl )−1 Gl ,the error of Gl is amplified by a small singular value at a large l. Many approximate methods
based on different principles have been proposed to mitigate this problem, such as the maximum entropy method [23] and the sparse modeling approach (SpM) [11, 19]. However, we
will not discuss this topic here.
The main focus of this article is to take the advantage of the fact that the Matsubara Green’s
function has less information than its real-frequency counterpart, and to speed up the computation within the imaginary-time/-frequency domain. In this respect, the IR basis Ûlα (iω) andUl (τ) efficiently represent the remaining information in imaginary-time quantities.
2
As a compact basis for imaginary-time Green’s functions, the Legendre basis [22] has been widely used in thecontext of DMFT and QMC. This fact is interesting because it shows that the Legendre basis corresponds to thehigh-temperature limit of the IR basis.
3
Note that the analytical form of the IR basis has yet to be determined as well. If you are good at mathematics,
please give it a try.
4
The convergence of ρl requires to establish ∃Vmax ∈ R ∀l ∈ N : ||Vl ||∞ ≤ Vmax . This is a conjecture supported
by numerical evidence.

5

SciPost Phys. Lect.Notes 63 (2022)

(a)

(b)

ρ(ω)

-1

+1

0

-1

10−1

10−1

10−4

10−4
|Gl |

10−7
−10

10

10

20

−10

10

30

+1

0

|Gl |

10−7

|ρl |
0

ρ(ω)

|ρl |
0

10

l

20

30

l

Figure 2: IR expansion coefficients Gl (Green’s function) and ρl (spectral function)for (a) a semicircular spectral function and (b) a delta function. We take β = 100,
ωmax = 1, and α = F. In (a), Gl and ρl for an odd l are not plotted because they arevalues of zero.
Finally, we present a numerical example that demonstrates the advantages of the IR basis
at low temperatures. Figure 3 shows the β-dependence of the number of expansion coefficients, N , required to reconstruct the Green’s function within a given accuracy for the modelin Fig. 2(b). Three kinds of expansions are compared: the IR-basis expansion, the Legendre
expansion, and the Matsubara (imaginary-frequency) expansion. As expected from Fig. 1(a),N increases logarithmically in the IR basis. By contrast, in the cases of thep imaginary-frequency
representation and the Legendre basis [22], N scales as O(β) and O( β), respectively. The
advantage of the IR basis becomes more pronounced at lower temperatures.

3

Sparse sampling

3.1

Dyson equation

In the previous section, we showed that the Matsubara Green’s function can be represented
compactly using the IR basis. A natural question is whether we can efficiently solve diagrammatic equations in the imaginary-time formalism. As an example, let us consider the Dyson
equation defined asĜ(iω) =1iω + µ − H − Σ̂(iω)

.

(11)

The self-energy Σ̂(iω) can also be compactly expanded with the IR basis:Σ̂(iω) 'L−1
X

Σl ÛlF (iω) .

(12)

l=0

Here, we excluded the Hartree term for simplicity. The question is whether the expansion
coefficients of Eq. (11), Gl , can be computed efficiently from the given Σl . In the following
demonstration, for simplicity, we set H = 0 and µ = 0.
6

SciPost Phys. Lect.Notes 63 (2022)

105
IR
Legendre basis
Matsubara

N

104
103
102
101 1
10

102

103

104

105

106

β
Figure 3: β–dependence of the number of basis functions, N , needed to reconstruct
the Green’s function with a high accuracy. We used the same spectral function as inFig. 2(b). We counted the minimum number of basis functions required to represent
G(τ = 0) with an accuracy of 10−8 . For the the imaginary-frequency representation,
we count the number of nonnegative imaginary frequencies. A second-order highfrequency expansion is applied in the Fourier transform into G(τ).
Solving the Dyson equation directly in the IR basis is not the most computationally efficientmethod. To see this, we transform the Dyson equation (11) into a linear equation for Ĝ(iω)
asXAiω,iω0 Ĝ(iω0 ) = 1 .
(13)
iω0

Here, the linear operator A is represented as a diagonal matrix defined byAiω,iω0 ≡ δiω,iω0 (iω − Σ̂(iω)). In principle, we can transform Aiω,iω0 and Ĝ(iω0 ) into the IR
basis. However, in the transformed linear equation, the linear operator is no longer diagonal,
and we need to solve the linear equation of size L × L. Therefore, the computational cost3
3of the Dyson equation scales as O(Nband
L 3 ) = O(Nband(log(βωmax ))3 ) for multiband systems.
This computational complexity grows more slowly than any power of βωmax . Although thisis asymptotically faster than the conventional approach, this method is not advantageous fortypical parameters in ab initio calculations, Nband = 100 and L ' 100.
How can we simultaneously benefit from the compactness of the representation and the
diagonality of the Dyson equation? We will answer this question in the next section.

3.2

Sparse sampling and Fourier transform

The sparsity of information can be exploited to reduce the computational complexity of the
diagrammatic equations either by the sparse sampling method [12] or the minimax isometry
method [24]. Before introducing the sparse sampling method, we discuss how to extract
information efficiently from physical quantities on the imaginary frequency axis.
The inverse transformation of Eq. (9) is given as
Gl =

i∞
X

ÛlF (iω)∗ Ĝ(iω) .

(14)

iω=−i∞

Here, we ignore ε L in Eq. (9). In principle, to determine Gl from this equation, we need toknow Ĝ(iω) in the whole imaginary-frequency domain. However, we do not need to evaluate
7

SciPost Phys. Lect.Notes 63 (2022)

(a)

(b)
2
F
(iωnF)Û39

U39(τ )

50
−5−G(τ )

10−20−1−102−101−10−10
0 10−1 100 101 102

τ
−G(τ
)

100

|∆G(τ )|

−4

|ImĜF(iωnF)|

10

0

10−5

10−8

10

|ImĜF(iωnF)||∆ImĜF(iωnF)|

10−10

−12

10−16

0

0

25

50
τ

75

10−15 2 1 0−10 −10 −10

100

0
ωnF

100 101 102

Figure 4: (Upper panel) The basis function at l = 39 and the sampling points (redcrosses) determined from the basis function. (Lower panel) The Green’s function
reconstructed by sparse sampling (solid line) and the deviation from the exact value(dashed line). The exact Green’s function is not shown because it is almost identicalto the reconstructed value. (a) Imaginary time domain and (b) imaginary frequency
domain. We used the semicircular model shown in Fig. 2(a) (β = 100, ωmax = 1).
Revised version of Fig. 2 in Ref. [12].

Sparse sampling points in imaginary timeG(τ̄k ) , Σ(τ̄k ) , P (τ̄k ) , W (τ̄k )
IR

F

F+
Arbitrary imaginary
time/frequency

Gl , Σl , Pl , Wl
F̂α

F̂+
α

Sparse sampling points in imaginary frequency

G(iω̄kF ), Σ(iω̄kF ), P (iω̄kB ), W (iω̄kB )Figure 5: Overview of the sparse sampling method. One can efficiently transform
data between sparse meshes in the imaginary frequency and time domains throughthe IR basis. The data can be evaluated at any frequency and time once the IR expansion coefficients are known.

8

SciPost Phys. Lect.Notes 63 (2022)

this equation in practice because the frequency dependence of Ĝ(iω) has only L degrees of
freedom at most, as shown in Fig. 2. Therefore, if we choose the appropriate L imaginary
frequencies iω and know Ĝ(iω) at those points, we can determine the values of Gl for all
l = 0, · · · , L − 1. We call these frequencies “sampling frequencies.”
The selection of sampling frequencies is not unique. In particular, we use the fact thatÛlα (iω) has a sign structure similar to a polynomial, and choose the maximum and minimumαneighborhoods of Û L−1(iω) (see Fig. 4):¯ α} .
¯ α , . . . , iω¯ α , iωW α = {iω
L
2
1

(15)

Given the values of the Green’s function on the sampling frequencies, we can easily calculate Glα using the least-squares fitting 5 as
Gl = argmin

XĜ(iω) −Gl=

iω∈W α+F̂α ĝ l .

L−1
X

2

Ûlα (iω)Gl

l=0

(16)

Here, F̂+α is the Moore-Penrose pseudo-inverse of F̂α , where F̂α is a matrix consisting of the¯ α ). By contrast, ĝ isvalues of the basis functions at the sampling frequencies: (F̂α )kl = Ûlα (iωk¯ α ). Because
a vector consisting of the values of the Green’s function to be fitted: (ĝ )k = Ĝ α (iωkF̂α is a small matrix of at most 100 × 100, as long as the sampling frequencies are chosen as
described above, the fitting procedure is numerically stable (refer to Sec. 3.3 and Appendix B).
Once F̂+α is computed, Eq. (16) is a matrix-vector product: The computational complexity is
2O(L ), and fast libraries such as BLAS can be used.
Similar sampling points in imaginary time can also be constructed by considering the sign
structure of Ul (τ). For the convention in Eq. (4), one can use the same sampling points forfermions and bosons:Gl = argmin

XGl

Ĝ(τ̄k ) −

L−1
X

2

Ûl (τ̄k )Gl

l=0

k


= F+ g l ,(17)where τ̄k are the sampling points and (F)kl ≡ Ul (τ̄k ).
Now, let us numerically demonstrate the accuracy of the sparse sampling method. The
upper panels in Figs. 4(a) and (b) show the IR basis functions and the sampling points at
ωmax = 1 and β = 100. It can be clearly seen that the basis functions have sign changes. The
imaginary-time sampling points are densely distributed near τ = 0 and β, where G(τ) varies
significantly [see the upper panel in Fig. 4(a)]. In the imaginary-frequency domain, the distribution of sampling points is sparse at high frequencies where the Green’s function essentiallyhas no structure [see the upper panel in Fig. 4(b)]. The lower panels in Figs. 4(a) and (b)show the Green’s functions reconstructed from Gl . The reconstructed Green’s functions agreewith the exact values with an accuracy of approximately 15 digits, indicating the accuracy andnumerical stability of the sparse sampling method.
Some may wonder why we can determine the expansion coefficient Gl of l < L − 1 acα
curately using the sampling points determined from the structure of U L−1 (τ) and U L−1(iω).
In fact, the distribution of the roots (zeros) of the IR basis functions for l < L − 1 is always
sparser than that of the basis functions for l = L − 1 (Figure 1)6 . Therefore, the sampling
5

Note that the L2 norm is not only the choice in the fitting. It is claimed that using the infinite norm leads tosmaller errors [25].
6More precisely, for l 0 < l, between two adjacent roots of Ul 0 (τ), Ul (τ) contains one or more roots. In the caseof l = l 0 + 1, this is mathematically proven. In a general case, no counterexample was found in the numericalexperiments.

9

SciPost Phys. Lect.Notes 63 (2022)

Table 1: Comparison between the sparse sampling method and the conventional
method (FFT). The power γ (≥ 1) depends on the details of the high-frequency expansion.

Sparse sampling
Conventional method (FFT)Memory cost
L = O(log(ωmax β))Nω = O(ωmax β)

Truncation errorO(e−αL )γ
O(1/Nω )Runtime cost
O(L 2 )
O(Nω log Nω )αpoints determined from U L−1 (τ) and U L−1(iω) capture the sign changes of the basis functions
for l < L − 1. For this reason, they serve as good sampling points for the truncated basis ofsize L.
We now have all ingredients of the sparse sampling method. Henceforth, the sampling
points at imaginary frequency and imaginary time will be referred to as the “sparse mesh”7 .
The sparse sampling method is summarized in Fig. 5. We can efficiently transform the dataon the sparse mesh (either in imaginary time or imaginary frequency) to the IR basis, and viceversa. Once the expansion coefficients Gl in the IR basis are obtained, the Green’s function
can be reconstructed at any frequency or time. In this sense, sparse sampling is a physicallymotivated interpolation method for imaginary-time and imaginary-frequency quantities.
The sparse sampling method allows us to solve the Dyson Eq. (21) by considering only the
sampling frequencies. The exact imaginary frequency/time dependence can be reproduced
from these values on the sparse mesh.
Traditionally, when one wants to transform data between imaginary frequency and time,
the fast Fourier transform is the standard method. Table 1 shows a comparison between theconventional method and the sparse sampling method. As the main feature of the sparse
sampling method and based on the data size, the computational complexity increases more
slowly than any power of β and bandwidth ωmax . Note that the sparse sampling method canalso be used in combination with conventional polynomial bases [12].

3.3Notes on calculations using the sparse sampling method

In practical applications using the sparse sampling method, one should make sure if ωmax and
L are large enough. If ωmax or the basis size L is not large enough, this may introduce large
systematic errors in the basis representation as defined in Eq. (9), such as large truncation
errors ε L . Such a systematic error may be amplified in the “fitting” procedure of Eqs. (16) and(17), leading to a large numerical error. Therefore, to ensure stable numerical calculations
with the sparse sampling method, one should choose the appropriate basis parameters L andΛ such that the basis representation stays accurate. The choice of the basis parameters (L,Λ) can be verified by checking if expansion coefficients Gl decay as fast as the singular values
down to the noise level. Otherwise, one should increase the value of Λ, e.g., by a factor of
10 or increase L. The basis size L must be chosen so that the S L−1 /S0 is at least the desired
accuracy for results multiplied by the condition number of the fitting matrix.

3.4

Related approaches

After establishing the idea of compressing the finite-T Green’s function using the spectral representation in [10], related approaches have been proposed. In this susbsection, we give a
brief overview on these approaches.
M. Kaltak and G. Kresse independently proposed a similar sampling method [24] by extending previously proposed zero-T optimal time and frequency grids [26, 27] to finite tem7

The minimax isometry method [24] uses a similar sparse mesh.

10

SciPost Phys. Lect.Notes 63 (2022)

perature. Their choice of sampling times and frequencies is based on a different principle(minimax isometry method). One interesting difference from the sparse sampling is that theirsampling points are not associated with an orthogonal basis.
An interesting follow-up study of IR has been done recently by J. Kaye et al. [20]. They
proposed to use simple exponentials in the time domain corresponding to δ functions in thereal-frequency domain as well as a systematic approach to choose sampling frequencies andtimes. Their “discrete Lehmann representation (DLR)” is generated by essentially the sameidea as of the IR, i.e., on a low rank decomposition of the kernel. DLR is slightly less compact
than IR and non-orthogonal, but has the same scaling with respect to βωmax . IR is unique fora given kernel and β, while DLR depends on the pivoting scheme and is not unique.
A detailed comparison of these related approaches would be an interesting topic for future
study.

3.5

Discrete Lehmann representation (DLR)

We explain the implementation of the DLR in our library. The poles on the real-frequency
axis selected for the DLR are based on a rank-revealing decomposition, which offers accuracy
guarantees. In our library, we instead select the pole locations based on the zeros of the IR
basis functions on the real axis. This is a heuristic motivated by an observation in a compact
discretization of an electron bath of a quantum impurity model [28]. We do not expect thatdifference to matter, which is however to be confirmed in a future study.
We first model the spectral function as
ρ(ω) =LX

c p δ(ω − ω̄ p ) ,

(18)

p=1
0
where sampling frequencies {ω̄1 , · · · , ω̄ L } are chosen to be the extrema of VL−1(ω). This
choice is heuristic but allows us a numerically stable transform between ρl and c p through therelation

ρl =

L
X

Vl p c p ,

(19)

p=1

where the matrix Vl p [≡ Vl (ω̄ p )] is well-conditioned. Figure 6(a) shows the sampling frequencies generated for β = 100 and ωmax = 1 (L = 40).
It is clear that for any given {g l } (l = 0, · · · , L − 1), one can find {ρl } (l =P
0, · · · , L − 1) thatL−1approximates the Green’s function within a desired accuracy by minimizing l=0 |g l + Sl ρl |2 .
Then, the pole coefficients {c p } are obtained from {ρl } using Eq. (19).
Figures 6(b) and 6(c) demonstrate the accuracy of DLR for the semicircular-DOS model[see Fig. 2(a)]. Figure 6(b) shows the expansion coefficients in DLR, c p , computed by fitting
the IR expansion coefficients Gl . The DLR expansion coefficients c p exhibits a dip aroundω = 0 where the distribution of the sampling frequencies is dense. As shown in Fig. 6(c),the IR coefficients reconstructed from the DLR coefficients match the numerically exact values
precisely.
Fitting a Green’s function with such a pole basis is equivalent to the infamously difficult
analytic continuation problem, i.e., reconstructing a spectral function. If the input is noisy,
one needs an appropriate regularization.

11

SciPost Phys. Lect.Notes 63 (2022)(a)

VL−1(ω)2
0−2

DLR coefficients cp−1.0(b) 0.08

0.0
ω

0.5

1.0

−0.5

0.0
ω

0.5

1.0

0.06
0.04
0.02
0.00
−1.0

(c)

−0.5

10−2
10−8

|Gl |

exact
DLR
error

10−14
10−20
0

5

10

15

20

25

l

Figure 6: Discrete Lehmann representation (DLR). (a) VL−1 (ω) computed for
β = 100 and ωmax = 1 (fermion, L = 40). The vertical lines denote sampling
frequencies. (b) Expansion coefficients in DLR for the Green’s function of the semicircular DOS model [see Fig. 2(a)]. (c) Comparison of IR expansion coefficientsreconstructed from the DLR data in (b) and the numerically exact result.

4

Review of applications

In the previous section, we showed that the expansion coefficients of the IR basis functions canbe calculated with high accuracy from the values of Green’s functions for a small number ofsampling points. This makes it possible to efficiently transform data between the imaginarytime and imaginary-frequency domains. In this section, we show how the sparse sampling
method can be used to efficiently solve the diagrammatic equations.

4.1

Noble-gas atoms, silicon crystals: Self-consistent GW and second-order perturbation (GF2) calculations

From Ref. [12], which proposed a sparse sampling method, we present the results of a selfconsistent GW calculation and a benchmark calculation using second-order perturbations(GF2). The purpose of this study was to verify the stability and accuracy of the numericalcalculations. Thus, they chose the physically “trivial” but technically difficult noble gas atoms
and silicon crystals. In this article, we present the results for noble gases. Please refer to theoriginal paper [12] for the results for silicon crystals.
12

SciPost Phys. Lect.Notes 63 (2022)

We consider the Hamiltonian
X
1 XX
†
† †
H=
hi j ciσ
c jσ +
Vi jkl ciσ
ckσ0 clσ0 c jσ ,2i jσi jkl σσ0(20)where hi j is the non-interacting Hamiltonian; Vi jkl is the Coulomb integral; i, j, k, and l areorbital indices; and σ and σ0 are spin indices. In this case, the Dyson equation reads:
F −1ˆ
,Ĝ(iωFn ) = [(iωFn + µ)I − F − Σ̃(iω
n )]

(21)

where the Fock matrix is given by F = h + ΣHF . The Hartree-Fock term can be computed asΣHF
i j = (2Vi jkl − Vilk j )ρkl ,

(22)

ˆ
F
where we used Einstein’s notation. We use different Σ̃(iωn ) for different approximations.
Based on the GF2 theory, we use(2)Σ̃i j (τ) ' Σi j (τ) = −Gkl (τ)Gqm (τ)Gnp (−τ)× Vikpq (2Vl jmn − Vm jln ) .

(23)

In the GW approximation, we takeΣ̃i j (τ) ' ΣGW
i j (τ) = −Glk (τ)W̃ilk j (τ) .

(24)Here, W̃ = W − V is the dynamic part of the screened interaction W :Ŵi jkl (iωBn ) = Vi jkl + Vi jpq P̂qpsr (iωBn )Ŵrskl (iωBn ) .

(25)

The bare polarization function is given by
Pi jkl (τ) = −Gil (τ)G jk (−τ) .

(26)

Equations (21)–(26) are diagonal in the imaginary frequency or imaginary time. Therefore, if the IR expansion coefficients of the response functions appearing on the right-handside of each of are known, the values of the left-hand side can be calculated quickly on thesampling points and then immediately converted into the IR basis. Thus, by going through the
IR basis, we can conduct self-consistent calculations while keeping the data compact.
Figure 7 shows the results of the calculation using the all-electron basis (cc-pVDZ) for
noble-gas atoms (β = 1000 Eh −1 , T ' 316 K). Owing to the deep core state of noble gas atoms,
it is difficult to apply all-electron calculations on a conventional polynomial basis, which requires thousands of basis functions even with the introduction of an effective potential. However, by combining the IR basis with sparse sampling, we can conduct all-electron calculations
with as few as 100 basis functions.
As the number of basis functions L increases, the total energy converges exponentially,
with relative errors easily converging below the level of 10−10 . For heavier noble-gas atoms,
the core states deepen. The IR basis can easily handle these deep core states by increasingΛ = βωmax . The results of the GW method, which involves the treatment of the response
functions of both fermions and bosons, also show convergence comparable to that of GF2.
This indicates the high numerical stability of the method.

13

SciPost Phys. Lect.Notes 63 (2022)

GF2
50

GW
100

50

100
He Λ = 104
Ne Λ = 105
Ar Λ = 106
Kr Λ = 106

10−1|∆Etot | [Eh ]

10−3

10−1
10−3
10−5

10−7

10−7

10−9

10−9

10−1

10−1

10−3

10−3

10−5

10−5

10−7

10−7

10−9

10−9

10−11

10−11

|∆Etot |/|Etot |

10−5

10−13

50

100

50

N

100

10−13

N

Figure 7: Convergence of total energy in GF2 (left panel) and GW (right panel) calculations for noble gas atoms. The horizontal axis is the number of IR bases. The
upper panel shows the absolute error of the energy (in Hartree Eh ), and the lower
panel shows the relative error. From Ref. [12].

4.2First-principles calculations of superconducting transition temperature
based on Migdal–Eliashberg equation

The transition temperature (Tc ) of a conventional (phonon-mediated) superconductor was recently obtained from first-principles calculations. The main theoretical approaches are the superconducting density functional theory and Migdal–Eliashberg theory. The latter has the advantage that the effective mass renormalization—owing to the electron-phonon interaction—
can be determined self-consistently. By contrast, a realistic prediction of the transition temperature is extremely difficult at low temperatures because it requires simultaneous treatmentof a phonon-mediated attraction and Coulomb repulsion, the energy scales of which are twoor three orders of magnitude apart. In this section, we introduce a recent study [29], whichtackles this problem using the IR basis.
In the framework of the Migdal-Eliashberg theory, Tc is obtained by solving the following
linearized gap equation:Xλ̃∆i (iωFn ) = −T
Ki j (iωFn − iωFn0 )F j (iωFn0 ) .
(27)
j,ωF 0nHere, subscript i denotes a combined index of the momentum and the energy band, ∆i isel−ph
the superconducting gap function, and Ki j = Ki j
+ Kiel−el
denotes the electron-phonon andj(screened) Coulomb interaction. Because the anomalous Green’s function Fi is given by
Fi (iωFn ) = |Gi (iωFn )|2 ∆i (iωFn ) .

14

(28)

SciPost Phys. Lect.Notes 63 (2022)

0.95(b) 0.95

0.9

0.9

0.85

0.85

~
λmax

~
λmax

(a)

0.8

0.75

0.75
0.7

0.8

LaH10
Nb
32

64

128

256

512

1024

2048

Number of Matsubara frequencies

0.7

LaH10
Nb
101

102

103

104

105

106

107

Λ

Figure 8: The convergence of the maximum eigenvalue λ̃max of the linearized gap
equation is just above the transition temperature. (a) Imaginary frequency + FFT and(b) IR basis + sparse modeling. The chain line in the figure shows the value obtainedusing the IR basis. This plot is essentially the same as Fig. 2 of Ref. [29]. The minor
differences in the values of λ̃max originate from differences in implementation details.
Equation (27) is an eigenvalue problem for a matrix, and λ̃ is its eigenvalue. The generalframework of the theory aims to solve Eq. (27) using Gi , including the self-energy correctionel−phby Ki j , to find Tc as the temperature at which the maximum value of λ̃ reaches 1.
Because the convolution integrals of the boson and fermion Green’s functions are diagonal
in imaginary time, the IR basis and the sparse sampling method can be applied to this problem.
In general, because realistic materials contain a large number of bands, and the number of kpoints required for convergence is large at low temperatures, it is difficult to solve Eq. (27).
Using the IR basis, we can expect a dramatic improvement, as shown in Table 1.
The actual eigenvalues of the linearized gap equation λ̃ for BCC-Nb and cubic LaH10 under high pressure (250 GPa) are shown in Fig. 8. The calculated temperatures are 19.7 and271.6 K, respectively, which are slightly higher than the respective transition temperatures.
Figure 8(a) shows the conventional method (imaginary frequency + FFT), which converges
roughly at Nω = 512 for LaH10 , but still does not converge at Nω = 2048 for Nb. As a roughestimate, if we assume that the upper limit of the band is ωmax ∼ 10 eV, which leads toβωmax ' 104 , we will need a higher Nω to converge for Nb. However, with the IR basis +sparse sampling method shown in Figure 8(b), we already have a fully converged value atΛ = 104 . Comparing the two calculations with Nω = 4096 (shown in Ref. [29]) and Λ = 105
as an example, we found that the number of bases is L = 138, which saves approximately
30-times the amount of memory and improves the actual calculation speed by approximately
20-fold. We note that we did not successfully obtain converged self-consistent solutions forΛ < 102 (LaH10 ) and Λ < 104 (Nb), which may be because Λ was not large enough (refer
to Appendix B). In the present calculation, the static shielding Coulomb interaction is used as
the electron-electron interaction, but it can be extended to include the plasmon effect in the
random phase approximation.

4.3

Calculation of magnetic interaction based on Liechtenstein formula

As is well known, density functional theory (DFT) is applicable only to the ground state, and
it is impossible to calculate the properties of finite temperatures without modification. In
the case of magnetic materials, the magnitude of the magnetic moment at absolute zero canbe predicted; however, the transition temperature cannot be calculated. Although there areseveral approaches to this problem, the most convenient and commonly used approach is to15

SciPost Phys. Lect.Notes 63 (2022)

derive an effective spin model from the results of the electronic structure calculation, andto obtain the transition temperature by a calculation using the spin model (e.g., mean-fieldapproximation and Monte Carlo calculation). Among the mappings to the spin model [30], a
method called the Liechtenstein method, has been found to significantly improve the efficiencyof the calculation by using the IR basis + sparse sampling method, which is briefly introducedin this subsection [13–15].
In thisPtheory, we compute the interaction parameters Ji j in the Heisenberg-type spin model
H = −2 〈i, j〉 Ji j Si · S j as follows:
Ji j = −
TX Tr Gi j (iωFn )M j G ji (iωFn )Mi .
2 ω

(29)

n

Here, Gi j is Green’s function for an effective tight-binding model derived from DFT calculaitions, and i, j are site indices. When the number of orbitals belonging to site i is Norb
, Gi j is an
j

i
Norb
× Norb matrix. In addition, although Mi is a perturbation term to the Hamiltonian owing
to the infinitesimal rotation of the spins, here it should simply be regarded as a constant matrixiiof Norb
× Norb, which does not depend on ωn . Moreover, Tr in Eq. (29) denotes the trace regarding the orbital component. Considering the computational cost required to evaluate thisterm, the dominant factor is the Fourier transform of the Green’s function from the k space
to the real space, which is approximately O(cNω ) using the usual Matsubara frequency mesh.
2Here, the coefficient factor c is approximately c = O((Nk log Nk )3 Norb
) when the number of kmeshes in one dimension is Nk . By contrast, if Green’s function calculation is limited to thesparse meshes, the computational cost becomes approximately O(c L), which is much lower at
low temperatures. In the actual calculation, the trace of Eq. (29) is calculated on the sparsemesh and finally converted into the IR basis. Subsequently, Eq. (29) is obtained by recovering
the values at the imaginary time τ = 0 from the expansion coefficients in the IR basis. In thismethod, the loss in the basis conversion is almost negligible, and the computation time can be
reduced by a factor of Nω /L.
In this study, although we used the Liechtenstein method as an example, we can generally
apply this method
for the free energy F [A] of a system described by the second-order HamilP
tonian H = 12 A12 c1† c2 , with the derivative δA12 F, δA12 δA34 F, · · · with respect to A. Although
this computation is not extremely heavy, it is a good example showing that the IR basis canbe used as a common tool. As we have seen, the convergence of the calculations using the IR
basis is extremely high, and the number of parameters required to check the convergence canbe reduced, which is always a significant advantage in research using numerical calculations.

4.4

Other applications

In addition to those introduced in this section, there are many other applications in the fieldof ab initio calculations, such as the calculation of the superconducting transition temperatureusing the FLEX approximation [17, 18], ab initio calculations of NiO compounds based on the
self-energy embedding method [16], and an ab initio estimation of the magnetic interaction
constant of NdNiO2 [15]. The exponential convergence of the calculation accuracy with respect to the basis size is also useful in model calculations, such as self-energy denoising [31]and efficient discretization of quantum impurity problems [28].

5 sparse-ir library
In the previous sections, we described the IR basis and sparse sampling method. We have developed the Python library sparse-ir [32], which is an updated version of irbasis [33]. This
16

SciPost Phys. Lect.Notes 63 (2022)

library can be used from Julia and Fortran as well. This library allows these new techniques
to be easily used without knowing the details of the computation of the IR basis functions andsampling points. Once the values of the basis functions and the sampling points are written
to a file, diagram calculations using sparse sampling can be conducted in any programming
language (the superconductivity calculations described above were actually implemented in
Fortran).
sparse-ir [32], has a more user-friendly interface than its predecessor irbasis [33].
While irbasis stores (precomputed) tabulated data of the IR basis functions for several valuesof Λ, sparse-ir allows on-the-fly computation of the IR basis functions at an arbitrary valueof Λ. The technical details will be explained in more depth in a follow-up paper.

5.1

Installation

Because sparse-ir is registered to PyPI, its installation is quite simple and can be achieved
by running the following command from a shell:$ pip install -U sparse-ir xprec
The -U option aims to update to the latest version of sparse-ir if it is already installed.
Though this is not strictly required, we strongly recommend installing the xprec package
alongside sparse-ir as it allows to compute the IR basis functions with greater accuracy.

5.2

Usage

sparse-ir can construct a basis object for given β and ωmax . Below is an example code that
construct a basis for Λ = 1000 (fermion) and β = 100 then evaluates the basis functions atcertain values of τ, ω, and imaginary frequency n.
# Compute IR basis for fermions and \ Lambda = \ beta * \ omega_max =1000
import sparse_ir
import numpy as np
lambda_ = 1000
eps = 1e - 8 # cut - off value for singular values
b_xy = sparse_ir . Dimensionles sBa si s ( ’F ’ , lambda_ , eps )# All singular values
print ( " singular values : " , b_xy . s )print ( " u_0 ( 0 . 1 ) " , b_xy . u [ 0 ] ( 0 . 1 ) )print ( " v_0 ( 0 . 1 ) " , b_xy . v [ 0 ] ( 0 . 1 ) )print ( "n - th derivative of u_l ( x ) and v_l ( y ) " )for n in range (1 , 3 ) :u_n = b_xy . u . deriv ( n )v_n = b_xy . v . deriv ( n )print ( " n = " , n , u_n [ 0 ] ( 0 . 1 ) )print ( " n = " , n , v_n [ 0 ] ( 0 . 1 ) )# Compute u_ { ln } as a matrix for the first# 10 non - nagative fermionic Matsubara frequencies# Fermionic / bosonic frequencies are denoted by odd / even integers .
hatF_t = b_xy . uhat ( 2 * np . arange ( 10 ) + 1 )print ( hatF_t . shape )When running this, you will instantly obtain the following results.
singular values : [6.93676632 e -02 6.39029304 e -02 4.73525981 e -02
3.78765556 e -02

17

SciPost Phys. Lect.Notes 63 (2022)

2.69709397 e -02 1.97919952 e -02 1.38987851 e -02
6.77586537 e -03 4.66563008 e -03 3.18095991 e -03
1.44908097 e -03 9.68433898 e -04 6.43222579 e -04
2.78944576 e -04 1.82210441 e -04 1.18411677 e -04
4.92756267 e -05 3.15628325 e -05 2.01254341 e -05
8.07529293 e -06 5.08266091 e -06 3.18587477 e -06
1.23669448 e -06 7.65992957 e -07 4.72633737 e -07
1.77929561 e -07 1.08573268 e -07 6.60143404 e -08
2.41478607 e -08 1.45293028 e -08 8.71232105 e -09
3.10136291 e -09 1.84126342 e -09 1.08960675 e -09
u_0 (0.1) 0.2740189634895232
v_0 (0.1) 0.7859154340971233
n - th derivative of u_l ( x ) and v_l ( y )n = 1 0.04705147674680676
n = 1 -5.044715458614174
n = 2 0.48534261358041686
n = 2 76.97205691155098(44 , 10)9.79170114 e -03
2.15479043 e -03
4.24765651 e -04
7.65700139 e -05
1.27758728 e -05
1.98887606 e -06
2.90529576 e -07
3.99959952 e -08
5.20671801 e -09
6.42730925 e -10]

One can evaluate the basis functions for multiple values of x, y, or n with a single function
call (corresponding to numpy broadcasting, see the examples below).
sparse-ir can construct a basis object for given β and ωmax directly, allows us to evaluate
Ul (τ), Vl (ω), Û α (iωn ) directly. This is more convenient in practical calculations. The interface
is quite similar to that for the dimensionless presentation. An example code is shown below.
# Compute IR basis for fermions and \ beta = 100 and \ omega_max = 10
lambda_ = 1000
beta = 100
wmax = lambda_ / beta
eps = 1e - 8 # cut - off value for singular values
b = sparse_ir . FiniteTempBasis ( ’F ’ , beta , wmax , eps = eps )x = y = 0.1
tau = 0 . 5 * beta * ( x + 1 )omega = wmax * y# All singular values
print ( " singular values : " , b . s )print ( " U_0 ( 0 . 1 ) " , b . u [ 0 ] ( tau ) )print ( " V_0 ( 0 . 1 ) " , b . v [ 0 ] ( omega ) )print ( "n - th derivative of U_l ( tau ) and V_l ( omega ) " )for n in range (1 , 3 ) :u_n = b . u . deriv ( n )v_n = b . v . deriv ( n )print ( " n = " , n , u_n [ 0 ] ( tau ) )print ( " n = " , n , v_n [ 0 ] ( omega ) )# Compute u_ { ln } as a matrix for the first# 10 non - nagative fermionic Matsubara frequencies# Fermionic / bosonic frequencies are denoted by odd / even integers .
hatF_t = b . uhat ( 2 * np . arange ( 10 ) + 1 )print ( hatF_t . shape )When running this, you will instantly obtain the following results.
singular values : [1.55110810 e +00 1.42891296 e +00 1.05883628 e +00
8.46945531 e -01
6.03088545 e -01 4.42562468 e -01 3.10786283 e -01 2.18949094 e -01
1.51512956 e -01 1.04326660 e -01 7.11284259 e -02 4.81825788 e -02
3.24024355 e -02 2.16548403 e -02 1.43828941 e -02 9.49804870 e -03
6.23739033 e -03 4.07434932 e -03 2.64776559 e -03 1.71215756 e -03

18

SciPost Phys. Lect.Notes 63 (2022)

1.10183651 e -03 7.05766389 e -04 4.50018387 e -04
1.80569039 e -04 1.13651753 e -04 7.12383254 e -05
2.76533293 e -05 1.71281232 e -05 1.05684116 e -05
3.97862594 e -06 2.42777207 e -06 1.47612553 e -06
5.39962581 e -07 3.24885087 e -07 1.94813421 e -07
6.93485829 e -08 4.11719016 e -08 2.43643475 e -08
U_0 (0.1) 0.038752133451430165
V_0 (0.1) 0.24852828200268673
n - th derivative of U_l ( x ) and V_l ( y )n = 1 0 .0 001 330 81 673 09 003 30 5
n = 1 -0.15952790996681684
n = 2 2.745512426092119 e -05
n = 2 0.24340701602860684(44 , 10)5.3

2.85677201 e -04
4.44726207 e -05
6.49643881 e -06
8.94337640 e -07
1.16425754 e -07
1.43719004 e -08]

Example: Second-order perturbation theory

As a simple example, we calculate the self-energy of the two-dimensional Hubbard model
within the range of second-order perturbations. The following code is available as a Jupyter
notebook [35].
First, for β = 103 and ωmax = 102 (Λ = 105 ), we construct a basis object (with 138 basis
functions). The generation of sampling points and the computation of F̂α and F as well as theirpseudo-inverse matrices are done by MatsubaraSampling and TauSampling.
import sparse_ir
import numpy as np
from numpy . fft import fftn , ifftn
beta = 1e + 3
lambda_ = 1e + 5
wmax = lambda_ / beta
eps = 1e - 15
print ( " wmax " , wmax )b = sparse_ir . FiniteTempBasis ( ’F ’ , beta , wmax , eps = eps )print ( " Number of basis functions " , b . size )# Sparse sampling in tau
smpl_tau = sparse_ir . TauSampling ( b )# Sparse sampling in Matsubara frequencies
smpl_matsu = sparse_ir . MatsubaraS ampling ( b )We now compute the Green’s function for a band dispersion:ε(k) = −2(cos k1 + cos k2 ) .

(30)

Note that the Green’ function
G(k, iω) =1iω − ε(k)is computed only on the sampling frequencies.
# Parameters
nk_lin = 64
U , kps = 2 .0 , np . array ( [ nk_lin , nk_lin ] )nw = smpl_matsu . sampling_points . size
ntau = smpl_tau . sampling_points . size

19

(31)

SciPost Phys. Lect.Notes 63 (2022)# Generate k mesh and non - interacting band energies
nk = np . prod ( kps )kgrid = [ 2 * np . pi * np . arange ( kp ) / kp for kp in kps ]k1 , k2 = np . meshgrid ( * kgrid , indexing = ’ ij ’)ek = - 2 * ( np . cos ( k1 ) + np . cos ( k2 ) )iw = 1j * np . pi * smpl_matsu . sampling_points / beta
# G ( iw , k ) : ( nw , nk )gkf = 1 . 0 / ( iw [ : , None ] - ek . ravel () [ None , : ] )At this point, gkf is a two-dimensional array of the number of sampling frequencies × thenumber of k-points. In imaginary time and real spaces, the second-order self-energy is givenbyΣ(τ, r) = U 2 G 2 (τ, r)G(β − τ, r) ,(32)

where r represents a position in real space. Thus, one can evaluate the self-energy by transforming Green’s function from the sampling frequencies to the IR basis and then to the sampling τ points as follows:# G (l , k ) : (L , nk )gkl = smpl_matsu . fit ( gkf )# G ( tau , k ) : ( ntau , nk )gkt = smpl_tau . evaluate ( gkl )# G ( tau , r ) : ( ntau , nk )grt = np . fft . fftn ( gkt . reshape ( ntau , * kps ) , axes = (1 , 2 ) ) .\reshape ( ntau , nk )# Sigma ( tau , r ) : ( ntau , nk )srt = U * U * grt * grt * grt [ : : -1 , : ]# Sigma (l , r ) : (L , nk )srl = smpl_tau . fit ( srt )(i , k = 0)

0.10
0.05
0.00
0.05
0.10
2000 1000

01000 2000

Figure 9: Self-energy computed for the Hubbard model on the sampling frequencies(k = 0). In sparse sampling methods, the self-energy is evaluated at a few sampling
frequencies.
20

SciPost Phys. Lect.Notes 63 (2022)# Sigma ( iw , r ) : ( nw , nk )srf = smpl_matsu . evaluate ( srl )At this point, the computed self-energy Σ(τ, r) is stored in a two-dimensional array srt.
Finally, we transform the self-energy into sampling frequencies (through the IR basis) andmomentum space as follows:# Sigma (l , r ) : (L , nk )srl = smpl_tau . fit ( srt )# Sigma ( iw , r ) : ( nw , nk )srf = smpl_matsu . evaluate ( srl )# Sigma ( iw , k ) : ( nw , kps [ 0 ] , kps [ 1 ])srf = srf . reshape ( nw , * kps )skf = ifftn ( srf , axes = (1 , 2 ) ) / nk ** 2Figure 9 plots the computed self-energy. The self-consistent calculations can be performedon sampling frequencies/imaginary times without using a dense mesh.

6

Summary and future directions

In this Lecture Note, we introduced a compact representation of Green’s function and a sparse
sampling scheme for solving diagrammatic equations. Using the IR basis, the data size increases only logarithmically with respect to the inverse temperature. Furthermore, the sparse
sampling method enables us to efficiently solve the diagrammatic equations, keeping the numerical data compact and numerically accurate. These methods are extremely useful in ab
initio calculations, where we need to deal with wide bandwidths and low temperatures. As
examples, we present the results of ab initio calculations using the GW approximation and theMigdal-Eliashberg theory. Furthermore, we reviewed the numerical library sparse-ir. We
hope that this Lecture Note will serve as a starting point for applications in various fields.
Finally, we discuss future directions. The frontier of method developments has alreadyshifted towards two-particle quantities. There are various theories at the two-particle level:
diagrammatic calculations with vertex corrections, dynamical susceptibility calculations basedon dynamical mean field theory, and nonlocal extensions of dynamical mean field theory that
can handle unconventional superconductivity (such as the dual fermion method [36] andthe dynamical vertex approximation approach [37]). These two-particle quantities have not
only three augments for imaginary frequencies but also arguments for orbitals, spins, andwavenumbers. To carry out calculations using these theories, we need an efficient numericaltreatment of two-particle quantities, which remains a challenging issue. The IR basis and thesparse sampling method have been extended to the two-particle quantities to address suchproblems. The interested readers are referred to recent articles [38–40].
We hope that the development of such fundamental numerical techniques will lead to technological breakthroughs and a better understanding of the properties of strongly correlated
compounds based on accurate ab initio calculations.

Acknowledgments
We would like to thank our collaborators in the research described in this article. Our collaboration with Prof. Masayuki Ohzeki, which aims to apply sparse modeling to quantum
many-body calculations, led to the discovery of the IR basis.

21

SciPost Phys. Lect.Notes 63 (2022)

Funding information This work was supported by the following Grant-in-Aid for Scientific Research: No. 15H05885 (J-Physics), No. 16H01064 (J-Physics), No. 18H04301 (JPhysics), No. 18H01158, No. 16K17735, No. 19K14654, No. 19K14654, No. 21H01003,
No. 21H01041. HS was supported by JST, PRESTO Grant No. JPMJPR2012. MW acknowledges support by the FWF (Austrian Science Funds) through Project No. P30997. EG and JL
were supported by the Simons foundation via the Simons Collaboration on the Many-Electron
Problem. TN was supported by JST, PRESTO Grant No. JPMJPR20L7.

A

Some properties of the IR basis functions

In this appendix, we are justifying some of the properties of the IR basis functions discussed
in Sec. 2. Let us do this for the fermionic kernel. For this, we will again write the kernel interms of dimensionless variablesexp − Λ2 x y,K(x, y) = −(A.1)
exp Λ2 y + exp − Λ2 ywhere again x = 2τ/β − 1, y = ω/ωmax , and Λ = βωmax .
R1
R1
This kernel is clearly square integrable, i.e., −1 dx −1 d y K 2 (x, y) < ∞, since its value isbounded. One can show that any such kernel admits a singular value expansion [41]:K(x, y) =∞Xul (x)sl vl∗ ( y) ,

(A.2)

l=0

where again sl are the singular values and ul and vl are the left and right singular functions,
respectively.
First, since the kernel is smooth for any temperature, K ∈ C ∞ [−1, 1]2 , one has that the
singular values decay faster than any power asymptotically [41]. Numerically, one can showthat sl in fact decay faster than exponentially, thus establishing Property 1.
The left and right singular functions in turn are the solution to eigenvalue equations withthe following integral kernels:
Z1
∞
2n
1
1
0
X
Λ2n X Ek ( 2 (x + 1))E2n−k ( 2 (x + 1))
u0
K (x, x ) =d y K(x, y)K(x 0 , y) =, (A.3)4n + 2 k=0
k!(2n − k)!
−1
n=0


Z1
tanh Λ2 y + tanh Λ2 y 0v
00
K ( y, y ) =dx K(x, y)K(x, y ) =,(A.4)
Λ( y + y 0 )−1where En (x) is the n-th Euler polynomial.
These kernels are real and symmetric by construction, which implies that the singular functions can be chosen to be real.
They are also centrosymmetric, i.e., one hasu0u0
K (x, x ) = K (−x, −x ). Since sl are not degenerate (see later), this implies that each ul (x)must be either an odd or an even function in x. The same is true for vl ( y). This establishes
Property 3.
The kernel −K(x, − y) belongs to the “exponential family” of parametrized probability distributions. These kernels satisfy a property called strict total positivity [42]. We briefly statethis property here: The set of k-samples Sk into the space [−1, 1] is defined as
Sk := {(x 1 , . . . , x k ) : −1 ≤ x 1 < . . . < x k ≤ 1}. A kernel K : [−1, 1]2 → R is called strictly
totally positive (STP) if for any k ≥ 0 and (x 1 , . . . , x k ), ( y1 , . . . yk ) ∈ Sk , one hasK(x 1 , y1 ) K(x 1 , y2 ) . . . K(x 1 , yk )
..
..
..
..
 > 0.
det 
(A.5)
.
.
.
.
K(x k , y1 ) K(x k , y2 ) . . . K(x k , yk )22

SciPost Phys. Lect.Notes 63 (2022)

Condition number

cond[F̂F ]104

Chebyshev
Λ = 1e5Λ = 1e6Λ = 1e7∼ N 1/2∼ N 3/2

103

102
102
N

103
102
cond[F ] (in τ )101cond[F̂F ] (in iωn )∼ Λ1/2

103

103

105

107

Λ

Figure 10: Condition number of the IR transformation matrices F̂F and F [12]. Leftpanel shows the condition number of frequency transformation matrices F̂F as a function of basis size N = L, in comparison with the Chebyshev representation. Rightpanel shows the condition number of both τ and iωn transformation matrices withrespect to Λ, where N is chosen to be the maximum number of coefficients with thesame cutoff in singular values Sl , provided in the irbasis library [21].
(This generalizes the notion of positive matrices appearing in, e.g., the theory of finite Markov
chains, to integral kernels.) One can now show that since K is a STP kernel, then the following
holds [42]:1. there exists a countable set of singular values which are all non-degenerate,
2. the n-th left/right singular function has exactly n nodal zeros in (−1, 1) and no zeros ofanother type,
3. the zeros of the n’th and n + 1’st singular function strictly interlace, i.e., exactly one zeroof un falls between two consecutive zeros of un+1 .
This establishes Property 4.

B

Numerical stability of sparse sampling

When the expansion coefficients are numerically evaluated using the sparse sampling “fitting”
procedures, i.e. Eqs. (16) and (17), numerical errors, such as round-off errors from floating
point operations or truncation errors from a finite basis cutoff, may be amplified due to the(pseudo-)inversion process. This error amplification can be quantified by the condition numberof the transformation matrices F and F̂α , defined as the product of the 2-norms of the matrix
and its inverse. In Fig. 10 we show the behaviors of such condition numbers for the IR basis as
a function of the basis size N = L (left panel, compared to the Chebyshev representation), andas a function of Λ (right panel). We can see that up to a significant number of basis functions,the condition numbers are < 104 , which indicates well-conditioned inversion problems. In
addition, the condition numbers show an approximate scaling of O(L 1/2 ), which is slower
that of the Chebyshev representation O(L 3/2 ). Since the values of L and Λ shown in Fig. 10
cover most values used in the calculations reviewed in this paper, the sparse sampling scheme
guarantees stable numerical routines to get accurate results.
23

SciPost Phys. Lect.Notes 63 (2022)

Λ = 10 6

Λ = 10 7
10 1

GW iteration 0
GW iteration 1
GW iteration 2
GW iteration 3
GW iteration 4

10 -1

kGl k

10 -3
10 -5
10 -7
10 -9
0

25

50

l

75

100

125 0

25

50

l

75

100

125

Figure 11: Examples of stable and unstable GW calculations of the Krypton atom
In practical applications of the sparse sampling method, one should take care not to introduce large systematic errors in the basis representation as defined in Eq. (9), such as large
truncation errors ε L due to insufficient basis size L or control parameter Λ. For example, asystematic error at the level of 10−3 in the basis representation, amplified by a condition number of 103 of the “fitting” procedure, may lead to a numerical error greater than the actual
result. As shown in Fig. 11, such a situation could make the simulation unstable. Therefore,
to ensure stable numerical calculations with the sparse sampling method, one should choose
the appropriate basis parameters L and Λ such that the basis representation stays accurate.
Figure 11 shows examples of stable and unstable GW calculation of the Krypton atom in
cc-pVDZ basis with β = 10 000 Eh −1 using sparse sampling. We carried out five GW iterations
with different choices of Λ, and plotted the norms of Gl from each iteration. The left panel ofFig. 11 shows results for Λ = 107 . At each iteration, the Green’s function is well approximatedby the IR basis and the basis truncation error is small, resulting in a stable simulation. In theright panel, a smaller Λ = 106 is used, which is insufficient for this system and introduces alarge systematic error (around 10−3 ), while the condition number of F is ∼ 103 . The systematic
error is amplified by fitting procedure, rendering the GW simulation unstable.

24

SciPost Phys. Lect.Notes 63 (2022)

References[1] H. N. Rojas, R. W. Godby and R. J. Needs, Space-time method for ab initio calculations ofself-energies and dielectric response functions of solids, Phys. Rev. Lett. 74, 1827 (1995),
doi:10.1103/physrevlett.74.1827.
[2] M. Kaltak, J. Klimeš and G. Kresse, Low scaling algorithms for the random phase approximation: Imaginary time and Laplace transformations, J. Chem. Theory Comput. 10, 2498
(2014), doi:10.1021/ct5001268.
[3] N. E. Bickers and D. J. Scalapino, Conserving approximations for strongly fluctuating
electron systems. I. Formalism and calculational approach, Ann. Phys. 193, 206 (1989),
doi:10.1016/0003-4916(89)90359-x.
[4] N. E. Bickers and S. R. White, Conserving approximations for strongly fluctuating electron systems. II. Numerical results and parquet extension, Phys. Rev. B 43, 8044 (1991),
doi:10.1103/physrevb.43.8044.
[5] 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.
[6] E. Gallicchio and B. J. Berne, The absorption spectrum of the solvated electron in fluid helium by maximum entropy inversion of imaginary time correlation functions from path integral Monte Carlo simulations, J. Chem. Phys. 101, 9909 (1994), doi:10.1063/1.467892.
[7] A. A. Kananenka, J. J. Phillips and D. Zgid, Efficient temperature-dependent Green’s functions methods for realistic systems: Compact grids for orthogonal polynomial transforms,
J. Chem. Theory Comput. 12, 564 (2016), doi:10.1021/acs.jctc.5b00884.
[8] M. Schüler and Y. Pavlyukh, Spectral properties from Matsubara Green’s function approach:Application to molecules, Phys. Rev. B 97, 115164 (2018),
doi:10.1103/physrevb.97.115164.
[9] T. Y. Takeshita, W. Dou, D. G. A. Smith, W. A. de Jong, R. Baer, D. Neuhauser and E.
Rabani, Stochastic resolution of identity second-order Matsubara Green’s function theory,
J. Chem. Phys. 151, 044114 (2019), doi:10.1063/1.5108840.
[10] 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.
[11] J. Otsuki, M. Ohzeki, H. Shinaoka and K. Yoshimi, Sparse modeling approach to analytical continuation of imaginary-time quantum Monte Carlo data, Phys. Rev. E 95, 061302(2017), doi:10.1103/physreve.95.061302.
[12] J. Li, M. Wallerberger, N. Chikano, C.-N. Yeh, 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.
[13] T. Nomoto, T. Koretsune and R. Arita, Local force method for the ab initio tight-binding
model: Effect of spin-dependent hopping on exchange interactions, Phys. Rev. B 102,014444 (2020), doi:10.1103/physrevb.102.014444.

25

SciPost Phys. Lect.Notes 63 (2022)[14] T. Nomoto, T. Koretsune and R. Arita, Formation mechanism of the helical Q
structure in Gd-based skyrmion materials, Phys. Rev. Lett. 125, 117204 (2020),
doi:10.1103/physrevlett.125.117204.
[15] Y. Nomura, T. Nomoto, M. Hirayama and R. Arita, Magnetic exchange coupling in cuprate-analog d 9 nickelates, Phys. Rev. Research 2, 043144 (2020),
doi:10.1103/physrevresearch.2.043144.
[16] S. Iskakov, C.-N. Yeh, E. Gull and D. Zgid, Ab initio self-energy embedding forthe photoemission spectra of NiO and MnO, Phys. Rev. B 102, 085105 (2020),
doi:10.1103/physrevb.102.085105.
[17] N. Witt, E. G. C. P. van Loon, T. Nomoto, R. Arita and T. O. Wehling, Efficient
fluctuation-exchange approach to low-temperature spin fluctuations and superconductivity: From the Hubbard model to N a x C oO2 · y H2 O, Phys. Rev. B 103, 205148 (2021),
doi:10.1103/PhysRevB.103.205148.
[18] N. Witt, J. M. Pizarro, T. Nomoto, R. Arita and T. O. Wehling, Doping fingerprints of spin
and lattice fluctuations in Moiré superlattice systems, Phys. Rev. B 105, L241109 (2022),
doi:10.1103/PhysRevB.105.L241109.
[19] 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.
[20] J. Kaye, K. Chen and O. Parcollet, Discrete Lehmann representation of imaginary time
Green’s functions, Phys. Rev. B 105, 235115 (2022), doi:10.1103/PhysRevB.105.235115.
[21] N. Chikano, J. Otsuki and H. Shinaoka, Performance analysis of a physically constructedorthogonal representation of imaginary-time Green’s function, Phys. Rev. B 98, 035104(2018), doi:10.1103/physrevb.98.035104.
[22] L. Boehnke, H. Hafermann, M. Ferrero, F. Lechermann and O. Parcollet, Orthogonal
polynomial representation of imaginary-time Green’s functions, Phys. Rev. B 84, 075145(2011), doi:10.1103/physrevb.84.075145.
[23] M. Jarrell and J. E. Gubernatis, Bayesian inference and the analytic continuation of imaginary-time quantum Monte Carlo data, Phys. Rep. 269, 133 (1996),
doi:10.1016/0370-1573(95)00074-7.
[24] M. Kaltak and G. Kresse, Minimax isometry method: A compressive sensing approachfor Matsubara summation in many-body perturbation theory, Phys. Rev. B 101, 205145(2020), doi:10.1103/physrevb.101.205145.
[25] A. Takatsuka, S. Ten-no and W. Hackbusch, Minimax approximation for the decomposition of energy denominators in Laplace-transformed Møller-Plesset perturbation theories, J.
Chem. Phys. 129, 044112 (2008), doi:10.1063/1.2958921.
[26] M. Kaltak, J. Klimeš and G. Kresse, Cubic scaling algorithm for the random phase approximation: Self-interstitials and vacancies in Si, Phys. Rev. B 90, 054115 (2014),
doi:10.1103/PhysRevB.90.054115.
[27] P. Liu, M. Kaltak, J. Klimeš and G. Kresse, Cubic scaling g w towards fast quasiparticle
calculations, Phys. Rev. B 94, 165109 (2016), doi:10.1103/PhysRevB.94.165109.
[28] H. Shinaoka and Y. Nagai, Sparse modeling of large-scale quantum impurity models withlow symmetries, Phys. Rev. B 103, 045120 (2021), doi:10.1103/physrevb.103.045120.
26

SciPost Phys. Lect.Notes 63 (2022)[29] T. Wang, T. Nomoto, Y. Nomura, H. Shinaoka, J. Otsuki, T. Koretsune and
R. Arita, Efficient ab initio Migdal-Eliashberg calculation considering the retardation effect in phonon-mediated superconductors, Phys. Rev. B 102, 134503 (2020),
doi:10.1103/physrevb.102.134503.
[30] A. I. Liechtenstein, M. I. Katsnelson and V. A. Gubanov, Exchange interactions andspin-wave stiffness in ferromagnetic metals, J. Phys. F: Met. Phys. 14, L125 (1984),
doi:10.1088/0305-4608/14/7/007.
[31] Y. Nagai and H. Shinaoka, Smooth self-energy in the exact-diagonalization-based dynamical
mean-field theory: Intermediate-representation filtering approach, J. Phys. Soc. Jpn. 88,
064004 (2019), doi:10.7566/jpsj.88.064004.
[32] M. Wallerberger et al., sparse-ir: Optimal compression and sparse sampling of many-body
propagators, arXiv:2206.11762v1.
[33] 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.
[34] H. Shinaoka,

https://nbviewer.jupyter.org/github/SpM-lab/ir_
kotaibutsuri/blob/82d34fe0bd7f5313fbdfa87a38983aadba2d85cb/sample.ipynb/.
Github,(2022),[35] H. Shinaoka,

Github, (2022), https://nbviewer.jupyter.org/github/SpM-lab/ir_
kotaibutsuri/blob/82d34fe0bd7f5313fbdfa87a38983aadba2d85cb/second_order_
perturbation.ipynb/.

[36] A. N. Rubtsov, M. I. Katsnelson and A. I. Lichtenstein, Dual fermion approach
to nonlocal correlations in the Hubbard model, Phys. Rev. B 77, 033101 (2008),
doi:10.1103/physrevb.77.033101.
[37] 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.
[38] 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.
[39] H. Shinaoka, D. Geffroy, M. Wallerberger, J. Otsuki, K. Yoshimi, E. Gull and J. Kuneš,
Sparse sampling and tensor network representation of two-particle Green’s functions, SciPost Phys. 8, 012 (2020), doi:10.21468/scipostphys.8.1.012.
[40] M. Wallerberger, H. Shinaoka and A. Kauch, Solving the Bethe-Salpeter equation with exponential convergence, Phys. Rev. Research 3, 033168 (2021),
doi:10.1103/physrevresearch.3.033168.
[41] P. C. Hansen, Discrete inverse problems, Society for Industrial and Applied Mathematics,
Philadelphia, US, ISBN 9780898716962 (2010), doi:10.1137/1.9780898718836.
[42] S. Karlin, Total positivity, Stanford University Press, Stanford, California, US, ISBN
9780804703147 (1968).

27