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, 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

*<EMAIL_ADDRESS>
## 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 the
Migdal–Eliashberg theory. We also describe a numerical library for the IR
basis and the sparse sampling method, sparse-ir, and provide its sample codes.
This lecture note follows the Japanese review article [H. Shinaoka et al.,
Solid State Physics 56(6), 301 (2021)].

###### Contents

  1. 1 Introduction
  2. 2 Intermediate-representation (IR) of Matsubara Green’s functions
  3. 3 Sparse sampling
    1. 3.1 Dyson equation
    2. 3.2 Sparse sampling and Fourier transform
    3. 3.3 Notes on calculations using the sparse sampling method
    4. 3.4 Related approaches
    5. 3.5 Discrete Lehmann representation (DLR)
  4. 4 Review of applications
    1. 4.1 Noble-gas atoms, silicon crystals: Self-consistent $GW$ and second-order perturbation (GF2) calculations
    2. 4.2 First-principles calculations of superconducting transition temperature based on Migdal–Eliashberg equation
    3. 4.3 Calculation of magnetic interaction based on Liechtenstein formula
    4. 4.4 Other applications
  5. 5 sparse-ir library
    1. 5.1 Installation
    2. 5.2 Usage
    3. 5.3 Example: Second-order perturbation theory
  6. 6 Summary and future directions
  7. A Some properties of the IR basis functions
  8. B Numerical stability of sparse sampling

##  1 Introduction

Perturbation and quantum field theories based on Green’s functions are widely
used for ab initio 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 been used for ab initio calculations
and model calculations based on various theories such as the $GW$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 to condensed matter physics but is also used in
quantum chemistry and high energy physics [7, 8, 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 is necessary 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 $\simeq$ 0.1 meV). In such
cases, the number of imaginary (Matsubara) frequencies required in the
calculations increases, and the 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 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,
14, 15, 16, 17, 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 library sparse-ir. We hope that this Lecture Note
help readers start using the IR basis and the sparse sampling method in
various fields.

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

The intermediate representation (IR) is a basis set in which the imaginary
frequency and time dependence 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(\tau)=-\langle T_{\tau}A^{\alpha}(\tau)B^{\alpha}(0)\rangle,$ (1)

where $\beta=1/T$ is the inverse temperature, $\tau$ is the imaginary time,
$T_{\tau}$ is the time ordering operator, and $\alpha$ is either F (fermion)
or B (boson). Here, we take $k_{\mathrm{B}}$ = 1. In addition, $A$ and $B$ are
the operators. In the important special case of a one-particle Green’s
function, $A$ is an annihilation operator, and $B$ is a creation operator.

The spectral representation of the two-point correlation function is given as
follows:

$G(\tau)=-\int_{-{\omega_{\mathrm{max}}}}^{\omega_{\mathrm{max}}}\differential{\omega^{\prime}}K^{\alpha}(\tau,\omega^{\prime})A(\omega^{\prime}),$
(2)

where $A(\omega)$ is a corresponding spectral function and the kernel
$K^{\alpha}(\tau,\omega)$ is defined as

$\displaystyle K^{\alpha}(\tau,\omega)$
$\displaystyle=\begin{cases}\displaystyle\frac{e^{-\tau\omega}}{1+e^{-\beta\omega}}&(\alpha=\mathrm{F})\\\
\displaystyle\frac{e^{-\tau\omega}}{1-e^{-\beta\omega}}&(\alpha=\mathrm{B})\end{cases},$
(3)

for $0<\tau<\beta$. We take ${\omega_{\mathrm{max}}}$ to be sufficiently large
such that the spectral function $A(\omega)$ is nonzero only for
$\omega\in[-{\omega_{\mathrm{max}}},{\omega_{\mathrm{max}}}]$. We can regard
Eq. (2) as an integral equation which connects the spectral function
$A(\omega)$ to the correlation function $G(\tau)$ through the kernel
$K^{\mathrm{\alpha}}(\tau,\omega)$.

To avoid the divergence of the bosonic kernel at $\omega=0$, we reformulate
Eq. (2) as [20]

$G(\tau)=-\int_{-{\omega_{\mathrm{max}}}}^{\omega_{\mathrm{max}}}\differential{\omega^{\prime}}K(\tau,\omega^{\prime})\rho(\omega^{\prime}),$
(4)

where we call $K(\tau,\omega)[=K^{\mathrm{F}}(\tau,\omega)]$ a logistic
kernel, and $\rho(\omega)$ is the modified spectral function defined by

$\displaystyle\rho(\omega)$
$\displaystyle\equiv\begin{cases}A(\omega)&(\alpha=\mathrm{F})\\\
\displaystyle\frac{A(\omega)}{\tanh(\beta\omega/2)}&(\alpha=\mathrm{B})\end{cases}.$
(5)

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].

The task at hand is to find a compact representation of the imaginary-
frequency (Matsubara) Green’s function $G(i\omega)$ that retains full
information of the real-frequency spectral function $\rho(\omega)$. For this,
we perform a singular value expansion of the logistic kernel $K$ for a given
$\beta$ and ${\omega_{\mathrm{max}}}$ [10]:

$K(\tau,\omega^{\prime})=\sum_{l=0}^{\infty}U_{l}(\tau)S_{l}V_{l}(\omega^{\prime}).$
(6)

The Fourier transform of $K(\tau,\omega^{\prime})$ with respect to $\tau$
yields the singular value expansion in the imaginary-frequency domain

$\hat{K}^{\alpha}(\mathrm{i}\omega^{(\alpha)},\omega^{\prime})\equiv-\int_{0}^{\beta}\differential\tau
K(\tau,\omega^{\prime}){e^{\mathrm{i}\omega^{(\alpha)}\tau}}=-\sum_{l=0}^{\infty}\hat{U}_{l}^{\alpha}(\mathrm{i}\omega^{(\alpha)})S_{l}V_{l}(\omega^{\prime}),$
(7)

where

$\displaystyle\hat{U}_{l}^{\alpha}({\mathrm{i}\omega^{(\alpha)}})$
$\displaystyle\equiv\int_{0}^{\beta}\differential\tau
U_{l}(\tau){e^{\mathrm{i}\omega^{(\alpha)}\tau}}.$ (8)

Here, $\mathrm{i}\omega^{(\alpha)}$ is the imaginary frequency corresponding
to the statistics $\alpha$ and we attached a hat $\hat{}$ to the quantities
defined in the imaginary-frequency domain. The minus sign in Eq. (7)
originates from the convention $K(\tau,\omega)>0$. In the following, we will
denote $i\omega^{(\alpha)}$ as $\mathrm{i}\omega$ for simplicity.

The singular value expansion has similar properties as the singular value
decomposition of matrices: The singular values $S_{l}$ satisfy
$S_{0}>S_{1}>\ldots>0$. The left singular functions
$\\{U_{0}(\tau),U_{1}(\tau),\ldots\\}$ form an orthonormal set on the
imaginary time axis $\tau\in[0,\beta]$, while the right singular functions
$\\{V_{0}(\omega),V_{1}(\omega),\ldots\\}$ form an orthonormal set on the real
frequency axis $\omega\in[-{\omega_{\mathrm{max}}},{\omega_{\mathrm{max}}}]$.

The so-called IR basis functions are nothing but the left singular functions
$U_{l}(\tau)$ [$\hat{U}_{l}^{\alpha}(\mathrm{i}\omega)$] and the right
singular functions $V_{l}(\omega)$. It should be emphasized that the IR basis
functions only depend on the inverse temperature $\beta$, statistics $\alpha$,
and the energy (frequency) cutoff ${\omega_{\mathrm{max}}}$ and do not depend
on the details of the system.

Figure 1: (a) Singular value $S_{l}$ computed for various values of $\beta$.
(b), (c) IR basis functions $U_{l}(\tau)$ and $V_{l}(\omega)$ computed for
$l=0,1,10$ and $\beta=100$, respectively. ${\omega_{\mathrm{max}}}$ is fixed
at ${\omega_{\mathrm{max}}}=1$.

Figure 1 shows the singular values $S_{l}$ and the IR basis functions
$U_{l}(\tau)$ and $V_{l}(\omega)$. The IR basis functions have the following
interesting properties:

Property 1

The singular values $S_{l}$ are non-degenerate, non-negative, and
monotonically decrease faster than exponentially with respect to $l$ 111For
numerical evidence, refer to Appendix D of [21]..

Property 2

The number of numerically significant singular values (e.g., $S_{l}/S_{0}\geq
10^{-15}$) is determined by the dimensionless quantity
$\Lambda\equiv\beta{\omega_{\mathrm{max}}}$ and increases only logarithmically
with respect to $\Lambda$ [21].

Property 3

$U_{l}(\tau)$ and $V_{l}(\omega)$ can be chosen to be real functions, and they
then become even (odd) functions for even (odd) $l$. In this convention,
$\hat{U}_{l}^{\alpha}(\mathrm{i}\omega)$ is purely imaginary or real.

Property 4

$U_{l}(\tau)$ and $V_{l}(\omega)$ have $l$ roots. In addition, in the limit of
$\beta{\omega_{\mathrm{max}}}\rightarrow 0$ (the high-temperature limit),
$U_{l}(\tau(x))$ and $V_{l}(\omega(x))$ coincide with the Legendre polynomial
$P_{l}(x)$ up to a constant [$\tau(x)\equiv\beta(x+1)/2$, $\omega(x)\equiv
x{\omega_{\mathrm{max}}}$ for $-1<x<1$] 222As a compact basis for imaginary-
time Green’s functions, the Legendre basis [22] has been widely used in the
context of DMFT and QMC. This fact is interesting because it shows that the
Legendre basis corresponds to the high-temperature limit of the IR basis..

We sketch mathematical proofs of Properties 1, 3, and 4 in Appendix A.
Property 2 has been verified only numerically333 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. .

Now, let us expand an “arbitrary” two-point correlation using a finite number
($L$) of the IR basis functions:

$\hat{G}(\mathrm{i}\omega)=\sum_{l=0}^{L-1}\hat{U}_{l}^{\alpha}(\mathrm{i}\omega)G_{l}+\epsilon_{L},$
(9)

where $\epsilon_{L}$ is a truncation error. If the corresponding spectral
function $\rho(\omega)$ is finite only in the interval
$[-{\omega_{\mathrm{max}}},{\omega_{\mathrm{max}}}]$, the expansion
coefficient $G_{l}$ can be evaluated using the orthonormality of the IR
functions as

$\displaystyle G_{l}$
$\displaystyle=-S_{l}\int_{-{\omega_{\mathrm{max}}}}^{\omega_{\mathrm{max}}}\differential\omega\rho(\omega)V_{l}(\omega)\equiv-
S_{l}\rho_{l}.$ (10)

From Eq. (10), it can be seen that $G_{l}$ vanishes at least as fast as
$S_{l}$,444 $G_{l}=\mathcal{O}(S_{l})$ is contingent on the boundedness of
$V_{l}$, i.e., $\exists V_{\mathrm{max}}\in\mathbb{R}\ \forall
l\in\mathbb{N}:||V_{l}||_{\infty}\leq V_{\mathrm{max}}$. This is a conjecture
supported by numerical evidence. i.e., faster than the exponential function
from Property 1.

Figure 2 shows the convergence of the expansion coefficients $\rho_{l}$ and
$G_{l}$ for some simple spectral functions $\rho(\omega)$. The convergence of
$\rho_{l}$ strongly depends on the shape of $\rho(\omega)$. In general,
$\rho_{l}$ decays exponentially in the case of a smooth spectral function, but
does not decay at all in the case of a $\delta$ function spectrum. By
contrast, $G_{l}$ always decays as fast as or even faster than the singular
values, regardless of the spectral functions. We would like to emphasize again
that this convergence is determined by the singular values of the kernel and
does not depend on the details of the system.

In Eq. (10), we can see that when we transform $\rho_{l}$ into $G_{l}$, the
contribution of $\rho_{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 $\rho_{l}=-\quantity(S_{l})^{-1}G_{l}$, the error of $G_{l}$ 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
$\hat{U}_{l}^{\alpha}(\mathrm{i}\omega)$ and $U_{l}(\tau)$ efficiently
represent the remaining information in imaginary-time quantities.

Finally, we present a numerical example that demonstrates the advantages of
the IR basis at low temperatures. Figure 3 shows the $\beta$-dependence of the
number of expansion coefficients, $N$, required to reconstruct the Green’s
function within a given accuracy for the model in 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 the
imaginary-frequency representation and the Legendre basis [22], $N$ scales as
$O(\beta)$ and $O(\sqrt{\beta})$, respectively. The advantage of the IR basis
becomes more pronounced at lower temperatures.

Figure 2:  IR expansion coefficients $G_{l}$ (Green’s function) and $\rho_{l}$
(spectral function) for (a) a semicircular spectral function and (b) a delta
function. We take $\beta=100$, ${\omega_{\mathrm{max}}}=1$, and
$\alpha=\mathrm{F}$. In (a), $G_{l}$ and $\rho_{l}$ for an odd $l$ are not
plotted because they are values of zero.  Figure 3:  $\beta$–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 in Fig. 2(b). We
counted the minimum number of basis functions required to represent
$G(\tau=0)$ with an accuracy of $10^{-8}$. For the the imaginary-frequency
representation, we count the number of nonnegative imaginary frequencies. A
second-order high-frequency expansion is applied in the Fourier transform into
$G(\tau)$.

##  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 as

$\displaystyle\hat{G}(\mathrm{i}\omega)$
$\displaystyle=\frac{1}{\mathrm{i}\omega+\mu-\mathcal{H}-\hat{\Sigma}(\mathrm{i}\omega)}.$
(11)

The self-energy $\hat{\Sigma}(\mathrm{i}\omega)$ can also be compactly
expanded with the IR basis:

$\displaystyle\hat{\Sigma}(\mathrm{i}\omega)$
$\displaystyle\simeq\sum_{l=0}^{L-1}\Sigma_{l}\hat{U}_{l}^{\mathrm{F}}(\mathrm{i}\omega).$
(12)

Here, we excluded the Hartree term for simplicity. The question is whether the
expansion coefficients of Eq. (11), $G_{l}$, can be computed efficiently from
the given $\Sigma_{l}$. In the following demonstration, for simplicity, we set
$\mathcal{H}=0$ and $\mu=0$.

Solving the Dyson equation directly in the IR basis is not the most
computationally efficient method. To see this, we transform the Dyson equation
(11) into a linear equation for $\hat{G}(\mathrm{i}\omega)$ as

$\displaystyle\sum_{\mathrm{i}\omega^{\prime}}A_{\mathrm{i}\omega,\mathrm{i}\omega^{\prime}}\hat{G}(\mathrm{i}\omega^{\prime})$
$\displaystyle=1.$ (13)

Here, the linear operator $A$ is represented as a diagonal matrix defined by
$A_{\mathrm{i}\omega,\mathrm{i}\omega^{\prime}}\equiv\delta_{\mathrm{i}\omega,\mathrm{i}\omega^{\prime}}(\mathrm{i}\omega-\hat{\Sigma}(\mathrm{i}\omega))$.
In principle, we can transform
$A_{\mathrm{i}\omega,\mathrm{i}\omega^{\prime}}$ and
$\hat{G}(\mathrm{i}\omega^{\prime})$ 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\times L$. Therefore, the
computational cost of the Dyson equation scales as
$O(N_{\mathrm{band}}^{3}L^{3})=O(N_{\mathrm{band}}^{3}(\log(\beta{\omega_{\mathrm{max}}}))^{3})$
for multiband systems. This computational complexity grows more slowly than
any power of $\beta{\omega_{\mathrm{max}}}$. Although this is asymptotically
faster than the conventional approach, this method is not advantageous for
typical parameters, $N_{\mathrm{band}}=100$ and $L\simeq 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

Figure 4:  (Upper panel) The basis function at $l=39$ and the sampling points
(red crosses) 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 identical to the reconstructed value. (a) Imaginary time domain
and (b) imaginary frequency domain. We used the semicircular model shown in
Fig. 2(a) ($\beta=100$, ${\omega_{\mathrm{max}}}=1$). Revised version of Fig.
2 in Ref. [12].  Figure 5: Overview of the sparse sampling method. One can
efficiently transform data between sparse meshes in the imaginary frequency
and time domains through the IR basis. The data can be evaluated at any
frequency and time once the IR expansion coefficients are known.

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

$\displaystyle G_{l}$
$\displaystyle=\sum_{\mathrm{i}\omega=-\mathrm{i}\infty}^{\mathrm{i}\infty}\hat{U}_{l}^{\mathrm{F}}(\mathrm{i}\omega)^{*}\hat{G}(\mathrm{i}\omega).$
(14)

Here, we ignore $\epsilon_{L}$ in Eq. (9). In principle, to determine $G_{l}$
from this equation, we need to know $\hat{G}(\mathrm{i}\omega)$ in the whole
imaginary-frequency domain. However, we do not need to evaluate this equation
in practice because the frequency dependence of $\hat{G}(\mathrm{i}\omega)$
has only $L$ degrees of freedom at most, as shown in Fig. 2. Therefore, if we
choose the appropriate $L$ imaginary frequencies $\mathrm{i}\omega$ and know
$\hat{G}(\mathrm{i}\omega)$ at those points, we can determine the values of
$G_{l}$ for all $l=0,\cdots,L-1$. We call these frequencies “sampling
frequencies.”

The selection of sampling frequencies is not unique. In particular, we use the
fact that $\hat{U}_{l}^{\alpha}(\mathrm{i}\omega)$ has a sign structure
similar to a polynomial, and choose the maximum and minimum neighborhoods of
$\hat{U}_{L-1}^{\alpha}(\mathrm{i}\omega)$ (see Fig. 4):

$\mathcal{W}^{\alpha}=\\{\bar{\mathrm{i}\omega}_{1}^{\alpha},\bar{\mathrm{i}\omega}_{2}^{\alpha},\ldots,\bar{\mathrm{i}\omega}_{L}^{\alpha}\\}.$
(15)

Given the values of the Green’s function on the sampling frequencies, we can
easily calculate $G^{\alpha}_{l}$ using the least-squares fitting 555Note that
the $L_{2}$ norm is not only the choice in the fitting. It is claimed that
using the infinite norm leads to smaller errors [25]. as

$\displaystyle G_{l}$
$\displaystyle=\underset{G_{l}}{\mathrm{argmin}}\sum_{\mathrm{i}\omega\in\mathcal{W}^{\alpha}}\bigg{|}\hat{G}(\mathrm{i}\omega)-\sum_{l=0}^{L-1}\hat{U}_{l}^{\alpha}(\mathrm{i}\omega)G_{l}\bigg{|}^{2}$
$\displaystyle=\pqty{\hat{\mathbf{F}}_{\alpha}^{+}\hat{\bm{g}}}_{l}.$ (16)

Here, $\hat{\mathbf{F}}_{\alpha}^{+}$ is the Moore-Penrose pseudo-inverse of
$\hat{\mathbf{F}}_{\alpha}$, where $\hat{\mathbf{F}}_{\alpha}$ is a matrix
consisting of the values of the basis functions at the sampling frequencies:
$(\hat{\mathbf{F}}_{\alpha})_{kl}=\hat{U}^{\alpha}_{l}(\bar{\mathrm{i}\omega}^{\alpha}_{k})$.
By contrast, $\hat{\bm{g}}$ is a vector consisting of the values of the
Green’s function to be fitted:
$(\hat{\bm{g}})_{k}=\hat{G}^{\alpha}(\bar{\mathrm{i}\omega}^{\alpha}_{k})$.
Because $\hat{\mathbf{F}}_{\alpha}$ is a small matrix of at most $100\times
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 $\hat{\mathbf{F}}_{\alpha}^{+}$ is computed, Eq. (16) is a matrix-vector
product: The computational complexity is $O(L^{2})$, 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 $U_{l}(\tau)$. For the convention in Eq.
(4), one can use the same sampling points for fermions and bosons:

$\displaystyle G_{l}$
$\displaystyle=\underset{G_{l}}{\mathrm{argmin}}\sum_{k}\bigg{|}\hat{G}(\bar{\tau}_{k})-\sum_{l=0}^{L-1}\hat{U}_{l}(\bar{\tau}_{k})G_{l}\bigg{|}^{2}$
$\displaystyle=\pqty{{\mathbf{F}}^{+}\bm{g}}_{l},$ (17)

where $\bar{\tau}_{k}$ are the sampling points and $({\mathbf{F}})_{kl}\equiv
U_{l}(\bar{\tau}_{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 ${\omega_{\mathrm{max}}}=1$ and $\beta=100$. It can be
clearly seen that the basis functions have sign changes. The imaginary-time
sampling points are densely distributed near $\tau=0$ and $\beta$, where
$G(\tau)$ 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 essentially has 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 $G_{l}$. The reconstructed Green’s
functions agree with the exact values with an accuracy of approximately 15
digits, indicating the accuracy and numerical stability of the sparse sampling
method.

Some may wonder why we can determine the expansion coefficient $G_{l}$ of
$l<L-1$ accurately using the sampling points determined from the structure of
$U_{L-1}(\tau)$ and $U_{L-1}^{\alpha}(\mathrm{i}\omega)$. 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)666More
precisely, for $l^{\prime}<l$, between two adjacent roots of
$U_{l^{\prime}}(\tau)$, $U_{l}(\tau)$ contains one or more roots. In the case
of $l=l^{\prime}+1$, this is mathematically proven. In a general case, no
counterexample was found in the numerical experiments.. Therefore, the
sampling points determined from $U_{L-1}(\tau)$ and
$U_{L-1}^{\alpha}(\mathrm{i}\omega)$ 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 of size $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”777The minimax isometry method [24] uses a similar sparse
mesh.. The sparse sampling method is summarized in Fig. 5. We can efficiently
transform the data on the sparse mesh (either in imaginary time or imaginary
frequency) to the IR basis, and vice versa. Once the expansion coefficients
$G_{l}$ 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
physically motivated 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.

| Memory cost | Truncation error | Runtime cost
---|---|---|---
Sparse sampling | $L=O(\log({\omega_{\mathrm{max}}}\beta))$ | $O(e^{-\alpha L})$ | $O(L^{2})$
Conventional method (FFT) | $N_{\omega}=O({\omega_{\mathrm{max}}}\beta)$ | $O(1/N_{\omega}^{\gamma})$ | $O(N_{\omega}\log N_{\omega})$
Table 1: Comparison between the sparse sampling method and the conventional
method (FFT). The power $\gamma$ ($\geq 1$) depends on the details of the
high-frequency expansion.

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 the conventional 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 $\beta$ and
bandwidth ${\omega_{\mathrm{max}}}$. Note that the sparse sampling method can
also be used in combination with conventional polynomial bases [12].

###  3.3 Notes on calculations using the sparse sampling method

In practical applications using the sparse sampling method, one should make
sure if ${\omega_{\mathrm{max}}}$ and $L$ are large enough. If
${\omega_{\mathrm{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 $\epsilon_{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 $\Lambda$ such that the basis
representation stays accurate. The choice of the basis parameters ($L$,
$\Lambda$) can be verified by checking if expansion coefficients $G_{l}$ decay
as fast as the singular values down to the noise level. Otherwise, one should
increase the value of $\Lambda$, e.g., by a factor of 10 or increase $L$. The
basis size $L$ must be chosen so that the $S_{L-1}/S_{0}$ 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 temperature. 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 their sampling 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 $\delta$ functions in the real-frequency domain as well as a
systematic approach to choose sampling frequencies and times. Their “discrete
Lehmann representation (DLR)” is generated by essentially the same idea 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
$\beta{\omega_{\mathrm{max}}}$. IR is unique for a given kernel and $\beta$,
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 that difference to matter, which is however to be confirmed in a future
study.

We first model the spectral function as

$\displaystyle\rho(\omega)$
$\displaystyle=\sum_{p=1}^{L}c_{p}\delta(\omega-\bar{\omega}_{p}),$ (18)

where sampling frequencies $\\{\bar{\omega}_{1},\cdots,\bar{\omega}_{L}\\}$
are chosen to be the extrema of $V^{\prime}_{L-1}(\omega)$. This choice is
heuristic but allows us a numerically stable transform between $\rho_{l}$ and
$c_{p}$ through the relation

$\displaystyle\rho_{l}$ $\displaystyle=\sum_{p=1}^{L}\bm{V}_{lp}c_{p},$ (19)

where the matrix $\bm{V}_{lp}~{}[\equiv V_{l}(\bar{\omega}_{p})]$ is well-
conditioned. Figure 6(a) shows the sampling frequencies generated for
$\beta=100$ and ${\omega_{\mathrm{max}}}=1$ ($L=40$).

It is clear that for any given $\\{g_{l}\\}$ ($l=0,\cdots,L-1$), one can find
$\\{\rho_{l}\\}$ ($l=0,\cdots,L-1$) that approximates the Green’s function
within a desired accuracy by minimizing
$\sum_{l=0}^{L-1}|g_{l}+S_{l}\rho_{l}|^{2}$. Then, the pole coefficients
$\\{c_{p}\\}$ are obtained from $\\{\rho_{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 $G_{l}$. The DLR
expansion coefficients $c_{p}$ exhibits a dip around $\omega=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.

Figure 6: Discrete Lehmann representation (DLR). (a) $V_{L-1}(\omega)$
computed for $\beta=100$ and ${\omega_{\mathrm{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 coefficients reconstructed 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 can be calculated with high accuracy from the values of
Green’s functions for a small number of sampling points. This makes it
possible to efficiently transform data between the imaginary-time 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 self-consistent $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 numerical calculations. 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 the original paper [12] for the results for silicon crystals.

We consider the Hamiltonian

$H=\sum_{ij\sigma}h_{ij}c^{\dagger}_{i\sigma}c_{j\sigma}+\frac{1}{2}\sum_{ijkl}\sum_{\sigma\sigma^{\prime}}V_{ijkl}c^{\dagger}_{i\sigma}c^{\dagger}_{k\sigma^{\prime}}c_{l\sigma^{\prime}}c_{j\sigma},$
(20)

where $h_{ij}$ is the non-interacting Hamiltonian; $V_{ijkl}$ is the Coulomb
integral; $i$, $j$, $k$, and $l$ are orbital indices; and $\sigma$ and
$\sigma^{\prime}$ are spin indices. In this case, the Dyson equation reads:

$\displaystyle\hat{G}(i\omega^{\mathrm{F}}_{n})$
$\displaystyle=[(i\omega^{\mathrm{F}}_{n}+\mu)I-F-\hat{\tilde{\Sigma}}(i\omega^{\mathrm{F}}_{n})]^{-1},$
(21)

where the Fock matrix is given by $F=h+\Sigma^{\mathrm{HF}}$. The Hartree-Fock
term can be computed as

$\Sigma^{\mathrm{HF}}_{ij}=(2V_{ijkl}-V_{ilkj})\rho_{kl},$ (22)

where we used Einstein’s notation. We use different
$\hat{\tilde{\Sigma}}(i\omega^{\mathrm{F}}_{n})$ for different approximations.
Based on the GF2 theory, we use

$\displaystyle\tilde{\Sigma}_{ij}(\tau)$
$\displaystyle\simeq\Sigma_{ij}^{(2)}(\tau)=-G_{kl}(\tau)G_{qm}(\tau)G_{np}(-\tau)$
$\displaystyle\times V_{ikpq}(2V_{ljmn}-V_{mjln}).$ (23)

In the $GW$approximation, we take

$\tilde{\Sigma}_{ij}(\tau)\simeq\Sigma_{ij}^{{GW}{}}(\tau)=-G_{lk}(\tau)\tilde{W}_{ilkj}(\tau).$
(24)

Here, $\tilde{W}=W-V$ is the dynamic part of the screened interaction $W$:

$\hat{W}_{ijkl}(i\omega^{\mathrm{B}}_{n})=V_{ijkl}+V_{ijpq}\hat{P}_{qpsr}(i\omega^{\mathrm{B}}_{n})\hat{W}_{rskl}(i\omega^{\mathrm{B}}_{n}).$
(25)

The bare polarization function is given by

$P_{ijkl}(\tau)=-G_{il}(\tau)G_{jk}(-\tau).$ (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-hand side of each of are known, the values of the left-
hand side can be calculated quickly on the sampling 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 ($\beta=$1000\text{\,}\mathrm{\text{
$E_{\textup{ h }}$ }}^{-1}$,T\simeq$316\text{\,}\mathrm{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
$\Lambda=\beta{\omega_{\mathrm{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.

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
$E_{\mathrm{h}}$), and the lower panel shows the relative error. From Ref.
[12].

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

The transition temperature ($T_{c}$) 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 treatment of a phonon-mediated attraction and Coulomb
repulsion, the energy scales of which are two or three orders of magnitude
apart. In this section, we introduce a recent study[29], which tackles this
problem using the IR basis.

In the framework of the Migdal-Eliashberg theory, $T_{c}$ is obtained by
solving the following linearized gap equation:

$\tilde{\lambda}\Delta_{i}(i\omega_{n}^{\mathrm{F}})=-T\sum_{j,\omega_{n^{\prime}}^{\mathrm{F}}}K_{ij}(i\omega_{n}^{\mathrm{F}}-i\omega_{n^{\prime}}^{\mathrm{F}})F_{j}(i\omega_{n^{\prime}}^{\mathrm{F}}).$
(27)

Here, subscript $i$ denotes a combined index of the momentum and the energy
band, $\Delta_{i}$ is the superconducting gap function, and $K_{ij}=K^{\rm el-
ph}_{ij}+K^{\rm el-el}_{ij}$ denotes the electron-phonon and (screened)
Coulomb interaction. Because the anomalous Green’s function $F_{i}$ is given
by

$F_{i}(i\omega_{n}^{\mathrm{F}})=|G_{i}(i\omega_{n}^{\mathrm{F}})|^{2}\Delta_{i}(i\omega_{n}^{\mathrm{F}}).$
(28)

Equation (27) is an eigenvalue problem for a matrix, and $\tilde{\lambda}$ is
its eigenvalue. The general framework of the theory aims to solve Eq. (27)
using $G_{i}$, including the self-energy correction by $K^{\rm el-ph}_{ij}$,
to find $T_{c}$ as the temperature at which the maximum value of
$\tilde{\lambda}$ 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 $k$-points 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.

Figure 8:  The convergence of the maximum eigenvalue
$\tilde{\lambda}_{\mathrm{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 obtained using
the IR basis. This plot is essentially the same as Fig. 2 of Ref. [29]. The
minor differences in the values of $\tilde{\lambda}_{\mathrm{max}}$ originate
from differences in implementation details.

The actual eigenvalues of the linearized gap equation $\tilde{\lambda}$ for
BCC-Nb and cubic LaH10 under high pressure (250 GPa) are shown in Fig. 8. The
calculated temperatures are 19.7 and 271.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_{\omega}=512$ for LaH10, but still does not converge at $N_{\omega}=2048$
for Nb. As a rough estimate, if we assume that the upper limit of the band is
$\omega_{\mathrm{max}}\sim 10$ eV, which leads to
$\beta{\omega_{\mathrm{max}}}\simeq 10^{4}$, we will need a higher
$N_{\omega}$ 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
$\Lambda=10^{4}$. Comparing the two calculations with $N_{\omega}=4096$ (shown
in Ref. [29]) and $\Lambda=10^{5}$ 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
$\Lambda<10^{2}$ (LaH10) and $\Lambda<10^{4}$ (Nb), which may be because
$\Lambda$ 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 can be predicted; however,
the transition temperature cannot be calculated. Although there are several
approaches to this problem, the most convenient and commonly used approach is
to derive an effective spin model from the results of the electronic structure
calculation, and to obtain the transition temperature by a calculation using
the spin model (e.g., mean-field approximation 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 efficiency of the
calculation by using the IR basis + sparse sampling method, which is briefly
introduced in this subsection [13, 14, 15].

In this theory, we compute the interaction parameters $J_{ij}$ in the
Heisenberg-type spin model $H=-2\sum_{\langle
i,j\rangle}J_{ij}{\bm{S}}_{i}\cdot{\bm{S}}_{j}$ as follows:

$J_{ij}=-\frac{T}{2}\sum_{\omega_{n}}\Tr[G_{ij}(i\omega_{n}^{\mathrm{F}})M_{j}G_{ji}(i\omega_{n}^{\mathrm{F}})M_{i}].$
(29)

Here, $G_{ij}$ is Green’s function for an effective tight-binding model
derived from DFT calculations, and $i,j$ are site indices. When the number of
orbitals belonging to site $i$ is $N_{\mathrm{orb}}^{i}$, $G_{ij}$ is an
$N_{\mathrm{orb}}^{i}\times N_{\mathrm{orb}}^{j}$ matrix. In addition,
although $M_{i}$ is a perturbation term to the Hamiltonian owing to the
infinitesimal rotation of the spins, here it should simply be regarded as a
constant matrix of $N_{\mathrm{orb}}^{i}\times N_{\mathrm{orb}}^{i}$, which
does not depend on $\omega_{n}$. Moreover, $\mathrm{Tr}$ in Eq. (29) denotes
the trace regarding the orbital component. Considering the computational cost
required to evaluate this term, 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_{\omega})$ using the usual Matsubara frequency mesh. Here,
the coefficient factor $c$ is approximately $c=O((N_{k}\log
N_{k})^{3}N_{\mathrm{orb}}^{2})$ when the number of $k$ meshes in one
dimension is $N_{k}$. By contrast, if Green’s function calculation is limited
to the sparse meshes, the computational cost becomes approximately $O{(cL)}$,
which is much lower at low temperatures. In the actual calculation, the trace
of Eq. (29) is calculated on the sparse mesh and finally converted into the IR
basis. Subsequently, Eq. (29) is obtained by recovering the values at the
imaginary time $\tau=0$ from the expansion coefficients in the IR basis. In
this method, the loss in the basis conversion is almost negligible, and the
computation time can be reduced by a factor of $N_{\omega}/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 Hamiltonian $H=\sum_{12}A_{12}c_{1}^{\dagger}c_{2}$, with
the derivative $\delta_{A_{12}}F,\;\delta_{A_{12}}\delta_{A_{34}}F,\;\cdots$
with respect to $A$. Although this computation is not extremely heavy, it is a
good example showing that the IR basis can be 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 can be
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 field of ab initio calculations, such as the calculation
of the superconducting transition temperature using 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 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 and
sampling 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 values of $\Lambda$, sparse-ir allows on-the-fly
computation of the IR basis functions at an arbitrary value of $\Lambda$. 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:

    $ python3 -m pip install sparse-ir xprec

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 $\beta$ and
${\omega_{\mathrm{max}}}$. Below is an example code that construct a basis for
$\Lambda=1000$ (fermion) and $\beta=100$ then evaluates the basis functions at
certain values of $\tau$, $\omega$, and imaginary frequency $n$.

import sparse_ir import numpy as np

# Compute IR basis for fermions and β= 100 and ω_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.55110810e+00 1.42891296e+00 1.05883628e+00 8.46945531e-01

6.03088545e-01 4.42562468e-01 3.10786283e-01 2.18949094e-01

1.51512956e-01 1.04326660e-01 7.11284259e-02 4.81825788e-02

3.24024355e-02 2.16548403e-02 1.43828941e-02 9.49804870e-03

6.23739033e-03 4.07434932e-03 2.64776559e-03 1.71215756e-03

1.10183651e-03 7.05766389e-04 4.50018387e-04 2.85677201e-04

1.80569039e-04 1.13651753e-04 7.12383254e-05 4.44726207e-05

2.76533293e-05 1.71281232e-05 1.05684116e-05 6.49643881e-06

3.97862594e-06 2.42777207e-06 1.47612553e-06 8.94337640e-07

5.39962581e-07 3.24885087e-07 1.94813421e-07 1.16425754e-07

6.93485829e-08 4.11719016e-08 2.43643475e-08 1.43719004e-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.00013308167309003305

n= 1 -0.15952790996681684

n= 2 2.745512426092119e-05

n= 2 0.24340701602860684

(44, 10)

###  5.3 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 [34].

First, for $\beta=10^{3}$ and ${\omega_{\mathrm{max}}}=10^{2}$
($\Lambda=10^{5}$), we construct a basis object (with 138 basis functions).
The generation of sampling points and the computation of
$\hat{\mathbf{F}}_{\alpha}$ and ${\mathbf{F}}$ as well as their pseudo-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.MatsubaraSampling(b)

We now compute the Green’s function for a band dispersion:

$\displaystyle\epsilon(\bm{k})$ $\displaystyle=-2(\cos{k_{1}}+\cos{k_{2}}).$
(30)

Note that the Green’ function

$\displaystyle G(\bm{k},\mathrm{i}\omega)$
$\displaystyle=\frac{1}{\mathrm{i}\omega-\epsilon(\bm{k})}$ (31)

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

# 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
$\times$ the number of $\bm{k}$-points. In imaginary time and real spaces, the
second-order self-energy is given by

$\displaystyle\Sigma(\tau,r)$
$\displaystyle=U^{2}G^{2}(\tau,r)G(\beta-\tau,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 $\tau$ 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)

# Sigma(iw, r): (nw, nk) srf = smpl_matsu.evaluate(srl)

At this point, the computed self-energy $\Sigma(\tau,r)$ is stored in a two-
dimensional array srt. Finally, we transform the self-energy into sampling
frequencies (through the IR basis) and momentum 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**2

Figure 9 plots the computed self-energy. The self-consistent calculations can
be performed on sampling frequencies/imaginary times without using a dense
mesh.

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.

##  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 the Migdal-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
already shifted towards two-particle quantities. There are various theories at
the two-particle level: diagrammatic calculations with vertex corrections,
dynamical susceptibility calculations based on dynamical mean field theory,
and nonlocal extensions of dynamical mean field theory that can handle
unconventional superconductivity (such as the dual fermion method [35] and the
dynamical vertex approximation approach [36]). These two-particle quantities
have not only three augments for imaginary frequencies but also arguments for
orbitals, spins, and wavenumbers. To carry out calculations using these
theories, we need an efficient numerical treatment of two-particle quantities,
which remains a challenging issue. The IR basis and the sparse sampling method
have been extended to the two-particle quantities to address such problems.
The interested readers are referred to recent articles [37, 38, 39].

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.

#### 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 (J-Physics),
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.

##  Appendix 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 in terms of dimensionless variables, cf.
Sec. LABEL:sec:dimensionless:

$K(x,y)=-\frac{\exp(-\frac{\Lambda}{2}xy)}{\exp(\frac{\Lambda}{2}y)+\exp(-\frac{\Lambda}{2}y)},$
(33)

where again $x=2\tau/\beta-1$, $y=\omega/{\omega_{\mathrm{max}}}$, and
$\Lambda=\beta{\omega_{\mathrm{max}}}$.

This kernel is clearly square integrable,
i.e.,$\int_{-1}^{1}\differential{x}\int_{-1}^{1}\differential{y}K^{2}(x,y)<\infty$,
since its value is bounded. One can show that any such kernel admits a
singular value expansion [40]:

$K(x,y)=\sum_{l=0}^{\infty}u_{l}(x)s_{l}v^{*}_{l}(y),$ (34)

where again $s_{l}$ are the singular values and $u_{l}$ and $v_{l}$ are the
left and right singular functions, respectively.

First, since the kernel is smooth for any temperature, $K\in
C^{\infty}[-1,1]^{2}$, one has that the singular values decay faster than any
power asymptotically [40]. Numerically, one can show that $s_{l}$ in fact
decay faster than exponentially, thus establishing Property 1.

The left and right singular functions in turn are the solution to eigenvalue
equations with the following integral kernels:

$\displaystyle K^{u}(x,x^{\prime})$ $\displaystyle=\int_{-1}^{1}\differential
yK(x,y)K(x^{\prime},y)=\sum_{n=0}^{\infty}\frac{\Lambda^{2n}}{4n+2}\sum_{k=0}^{2n}\frac{E_{k}(\frac{1}{2}(x+1))E_{2n-k}(\frac{1}{2}(x^{\prime}+1))}{k!(2n-k)!},$
(35) $\displaystyle K^{v}(y,y^{\prime})$
$\displaystyle=\int_{-1}^{1}\differential
xK(x,y)K(x,y^{\prime})=\frac{\tanh(\tfrac{\Lambda}{2}y)+\tanh(\tfrac{\Lambda}{2}y^{\prime})}{\Lambda(y+y^{\prime})},$
(36)

where $E_{n}(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 has $K^{u}(x,x^{\prime})=K^{u}(-x,-x^{\prime})$. Since $s_{l}$ are
not degenerate (see later), this implies that each $u_{l}(x)$ must be either
an odd or an even function in $x$. The same is true for $v_{l}(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 [41]. We briefly state this property here: The set of
$k$-samples $\mathcal{S}_{k}$ into the space $[-1,1]$ is defined as
$\mathcal{S}_{k}:=\\{(x_{1},\ldots,x_{k}):-1\leq x_{1}<\ldots<x_{k}\leq 1\\}$.
A kernel $K:[-1,1]^{2}\to\mathbb{R}$ is called strictly totally positive (STP)
if for any $k\geq 0$ and $(x_{1},\ldots,x_{k}),(y_{1},\ldots
y_{k})\in\mathcal{S}_{k}$, one has

$\det\left(\begin{array}[]{cccc}K(x_{1},y_{1})&K(x_{1},y_{2})&\ldots&K(x_{1},y_{k})\\\
\vdots&\vdots&\ddots&\vdots\\\
K(x_{k},y_{1})&K(x_{k},y_{2})&\ldots&K(x_{k},y_{k})\\\ \end{array}\right)>0.$
(37)

(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 [41]:

  1. 1.

there exists a countable set of singular values which are all non-degenerate,

  2. 2.

the $n$-th left/right singular function has exactly $n$ nodal zeros in
$(-1,1)$ and no zeros of another type,

  3. 3.

the zeros of the $n$’th and $n+1$’st singular function strictly interlace,
i.e., exactly one zero of $u_{n}$ falls between two consecutive zeros of
$u_{n+1}$.

This establishes Property 4.

##  Appendix 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 number of the
transformation matrices ${\mathbf{F}}$ and $\hat{\mathbf{F}}_{\alpha}$,
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), and as a function of $\Lambda$ (right panel). We can see that
up to a significant number of basis functions, the condition numbers are
$<10^{4}$, 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 $\Lambda$ 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.

Figure 10: Condition number of the IR transformation matrices
$\hat{\mathbf{F}}_{\mathrm{F}}$ and ${\mathbf{F}}$ [12]. Left panel shows the
condition number of frequency transformation matrices
$\hat{\mathbf{F}}_{\mathrm{F}}$ as a function of basis size $N=L$, in
comparison with the Chebyshev representation. Right panel shows the condition
number of both $\tau$ and $i\omega_{n}$ transformation matrices with respect
to $\Lambda$, where $N$ is chosen to be the maximum number of coefficients
with the same cutoff in singular values $S_{l}$, provided in the irbasis
library [21].

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 $\epsilon_{L}$ due to
insufficient basis size $L$ or control parameter $\Lambda$. For example, a
systematic error at the level of $10^{-3}$ in the basis representation,
amplified by a condition number of $10^{3}$ 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 $\Lambda$ 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 $\beta=$10\,000\text{\,}\mathrm{\text{
$E_{\textup{ h }}$ }}^{-1}$$ using sparse sampling. We carried out five $GW$
iterations with different choices of $\Lambda$, and plotted the norms of
$G_{l}$ from each iteration. The left panel of Fig. 11 shows results for
$\Lambda=10^{7}$. At each iteration, the Green’s function is well approximated
by the IR basis and the basis truncation error is small, resulting in a stable
simulation. In the right panel, a smaller $\Lambda=10^{6}$ is used, which is
insufficient for this system and introduces a large systematic error (around
$10^{-3}$), while the condition number of ${\mathbf{F}}$ is $\sim 10^{3}$. The
systematic error is amplified by fitting procedure, rendering the $GW$
simulation unstable.

Figure 11:  Examples of stable and unstable $GW$ calculations of the Krypton
atom

## References

  * [1] H. N. Rojas, R. W. Godby and R. J. Needs,  _Space-time method for ab initio calculations of self-energies and dielectric response functions of solids_ ,  Physical Review Letters 74(10), 1827 (1995),  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_ ,  Journal of Chemical Theory and Computation 10(6), 2498 (2014),  10.1021/ct5001268.
  * [3] N. E. Bickers and D. J. Scalapino,  _Conserving approximations for strongly fluctuating electron systems. I. Formalism and calculational approach_ ,  Annals of Physics 193(1), 206 (1989),  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_ ,  Physical Review B 43(10), 8044 (1991),  10.1103/physrevb.43.8044.
  * [5] A. Georges, G. Kotliar, W. Krauth and M. J. Rozenberg,  _Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions_ ,  Reviews of Modern Physics 68(1), 13 (1996-01),  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_ ,  The Journal of Chemical Physics 101(11), 9909 (1994),  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_ ,  Journal of Chemical Theory and Computation 12(2), 564 (2016),  10.1021/acs.jctc.5b00884.
  * [8] M. Schüler and Y. Pavlyukh,  _Spectral properties from Matsubara Green’s function approach: Application to molecules_ ,  Physical Review B 97(11), 115164 (2018),  10.1103/physrevb.97.115164,  1710.08660.
  * [9] T. Y. Takeshita, W. Dou, D. G. A. Smith, W. A. d. Jong, R. Baer, D. Neuhauser and E. Rabani,  _Stochastic resolution of identity second-order Matsubara Green’s function theory_ ,  The Journal of Chemical Physics 151(4), 044114 (2019),  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_ ,  Physical Review B 96(3), 035147 (2017),  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_ ,  Physical Review E 95(6), 061302(R) (2017),  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_ ,  Physical Review B 101(3), 035144 (2020),  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_ ,  Physical Review B 102(1), 014444 (2020),  10.1103/physrevb.102.014444,  2003.11162.
  * [14] T. Nomoto, T. Koretsune and R. Arita,  _Formation Mechanism of the Helical Q Structure in Gd-Based Skyrmion Materials_ ,  Physical Review Letters 125(11), 117204 (2020),  10.1103/physrevlett.125.117204,  2003.13167.
  * [15] Y. Nomura, T. Nomoto, M. Hirayama and R. Arita,  _Magnetic exchange coupling in cuprate-analog d9 nickelates_ ,  Physical Review Research 2(4), 043144 (2020),  10.1103/physrevresearch.2.043144,  2006.16943.
  * [16] S. Iskakov, C.-N. Yeh, E. Gull and D. Zgid,  _Ab initio self-energy embedding for the photoemission spectra of NiO and MnO_ ,  Physical Review B 102(8), 085105 (2020),  10.1103/physrevb.102.085105,  2003.04440.
  * [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 ${\mathrm{na}}_{x}{\mathrm{coo}}_{2}\cdot{}{y\mathrm{H}}_{2}\mathrm{O}$_,  Phys. Rev. B 103, 205148 (2021),  10.1103/PhysRevB.103.205148,  2012.04562.
  * [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_ (2021), https://arxiv.org/abs/2108.01121.
  * [19] J. Otsuki, M. Ohzeki, H. Shinaoka and K. Yoshimi,  _Sparse Modeling in Quantum Many-Body Problems_ ,  Journal of the Physical Society of Japan 89(1), 012001 (2020),  10.7566/jpsj.89.012001.
  * [20] J. Kaye, K. Chen and O. Parcollet,  _Discrete lehmann representation of imaginary time green’s functions_ (2021), 2107.13094.
  * [21] N. Chikano, J. Otsuki and H. Shinaoka,  _Performance analysis of a physically constructed orthogonal representation of imaginary-time Green’s function_ ,  Physical Review B 98(3), 035104 (2018),  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_ ,  Physical Review B 84(7), 075145 (2011-08),  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_ ,  Physics Reports 269(3), 133 (1996),  10.1016/0370-1573(95)00074-7.
  * [24] M. Kaltak and G. Kresse,  _Minimax isometry method: A compressive sensing approach for Matsubara summation in many-body perturbation theory_ ,  Physical Review B 101(20), 205145 (2020),  10.1103/physrevb.101.205145,  1909.01740.
  * [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_ ,  The Journal of Chemical Physics 129(4), 044112 (2008),  10.1063/1.2958921,  https://doi.org/10.1063/1.2958921.
  * [26] M. Kaltak, J. c. v. Klimeš and G. Kresse,  _Cubic scaling algorithm for the random phase approximation: Self-interstitials and vacancies in si_ ,  Phys. Rev. B 90, 054115 (2014),  10.1103/PhysRevB.90.054115.
  * [27] P. Liu, M. Kaltak, J. c. v. Klimeš and G. Kresse,  _Cubic scaling $gw$: Towards fast quasiparticle calculations_,  Phys. Rev. B 94, 165109 (2016),  10.1103/PhysRevB.94.165109.
  * [28] H. Shinaoka and Y. Nagai,  _Sparse modeling of large-scale quantum impurity models with low symmetries_ ,  Physical Review B 103(4), 045120 (2021),  10.1103/physrevb.103.045120.
  * [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_ ,  Physical Review B 102(13), 134503 (2020),  10.1103/physrevb.102.134503.
  * [30] A. I. Liechtenstein, M. I. Katsnelson and V. A. Gubanov,  _Exchange interactions and spin-wave stiffness in ferromagnetic metals_ ,  Journal of Physics F: Metal Physics 14(7), L125 (1984),  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_ ,  Journal of the Physical Society of Japan 88(6), 064004 (2019),  10.7566/jpsj.88.064004.
  * [32] M. Wallerberger, S. Badr, S. Hoshino, F. Kakizawa, T. Koretsune, Y. Nagai, K. Nogaki, T. Nomoto, H. Mori, J. Otsuki, S. Ozaki, R. Sakurai _et al._ ,  _sparse-ir: optimal compression and sparse sampling of many-body propagators_ ,  10.48550/arxiv.2206.11762 (2022).
  * [33] N. Chikano, K. Yoshimi, J. Otsuki and H. Shinaoka,  _irbasis: Open-source database and software for intermediate-representation basis functions of imaginary-time Green’s function_ ,  Computer Physics Communications 240, 181 (2018),  10.1016/j.cpc.2019.02.006,  1807.05237.
  * [34] https://nbviewer.jupyter.org/github/SpM-lab/ir_kotaibutsuri/blob/82d34fe0bd7f5313fbdfa87a38983aadba2d85cb/second_order_perturbation.ipynb/.
  * [35] A. N. Rubtsov, M. I. Katsnelson and A. I. Lichtenstein,  _Dual fermion approach to nonlocal correlations in the Hubbard model_ ,  Physical Review B 77(3), 033101 (2008),  10.1103/physrevb.77.033101.
  * [36] A. Toschi, A. Katanin and K. Held,  _Dynamical vertex approximation: A step beyond dynamical mean-field theory_ ,  Physical Review B 75(4), 045118 (2007-01),  10.1103/physrevb.75.045118.
  * [37] H. Shinaoka, J. Otsuki, K. Haule, M. Wallerberger, E. Gull, K. Yoshimi and M. Ohzeki,  _Overcomplete compact representation of two-particle Green’s functions_ ,  Physical Review B 97(20), 205111 (2018),  10.1103/physrevb.97.205111.
  * [38] 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 Physics 8(1), 012 (2020),  10.21468/scipostphys.8.1.012.
  * [39] M. Wallerberger, H. Shinaoka and A. Kauch,  _Solving the Bethe-Salpeter equation with exponential convergence_ ,  Physical Review Research 3(3), 033168 (2021),  10.1103/physrevresearch.3.033168.
  * [40] P. C. Hansen,  _Discrete Inverse Problems: Insights and Algorithms_ ,  SIAM,  ISBN 978-0-898716-96-2,  10.1137/1.9780898718836 (2010).
  * [41] S. Karlin,  _Total Positivity_ ,  Stanford University Press,  ISBN 978-0804703147 (1968).