# Robust analytic continuation of Green’s functions via projection, pole
estimation, and semidefinite relaxation

Zhen Huang  Department of Mathematics, University of California, Berkeley,
California 94720 USA  Emanuel Gull  Department of Physics, University of
Michigan, Ann Arbor, Michigan 48109, USA  Lin Lin<EMAIL_ADDRESS>Department of Mathematics, University of California, Berkeley, California
94720 USA  Applied Mathematics and Computational Research Division, Lawrence
Berkeley National Laboratory, Berkeley, CA 94720, USA

(August 28, 2024)

###### Abstract

Green’s functions of fermions are described by matrix-valued Herglotz-
Nevanlinna functions. Since analytic continuation is fundamentally an ill-
posed problem, the causal space described by the matrix-valued Herglotz-
Nevanlinna structure can be instrumental in improving the accuracy and in
enhancing the robustness with respect to noise. We demonstrate a three-pronged
procedure for robust analytic continuation called PES: (1) Projection of data
to the causal space. (2) Estimation of pole locations. (3) Semidefinite
relaxation within the causal space. We compare the performance of PES with the
recently developed Nevanlinna and Carathéodory continuation methods and find
that PES is more robust in the presence of noise and does not require the
usage of extended precision arithmetics. We also demonstrate that a causal
projection improves the performance of the Nevanlinna and Carathéodory
methods. The PES method is generalized to bosonic response functions, for
which the Nevanlinna and Carathéodory continuation methods have not yet been
developed. It is particularly useful for studying spectra with sharp features,
as they occur in the study of molecules and band structures in solids.

##  I Introduction

Green’s functions of quantum many-body systems are central subjects in
condensed matter physics, quantum chemistry, and quantum field theory. Due to
the difficulty of performing correlated finite-temperature calculations on the
real axis, numerical methods typically obtain results on the Matsubara axis
and employ analytic continuation as a postprocessing tool to infer Green’s
functions along the real axis. Analytic continuation is fundamentally an ill-
posed problem, and many numerical methods have been developed, such as the
Padé fit [1, 2], maximum entropy (MaxEnt) [3, 4, 5, 6, 7], stochastic analytic
continuation and its variants [8, 9, 10, 11], sparse modeling [12, 13], and
machine learning approaches [14, 15]. All these methods encounter difficulties
in satisfying causality, recovering sharp features, resolving multiple
features, and/or capturing high frequency information in the spectra.
Furthermore, most methods have been developed only for diagonal entries of the
Green’s function. The analytic continuation of off-diagonal entries of the
Green’s function [16, 17, 18, 19, 20] is often attempted by diagonalizing the
Green’s function at a certain Matsubara frequency, conducting the analytic
continuation only of the diagonal entries with respect to the transformed
basis at other frequencies, and neglecting the remaining non-diagonal entries.

The fermionic Green’s functions is a special type of matrix-valued Herglotz-
Nevanlinna functions. These functions have a wide range of scientific and
engineering applications, and have been (somewhat inconsistently) named after
renowned mathematicians including Carathéodory, Herglotz, Nevanlinna, Pick,
and Riesz. They are also sometimes called $R$-functions [21]. The name matrix-
valued Herglotz-Nevanlinna functions adopted in this paper follows the
suggestion in [22]. This crucial analytic structure has only been taken into
account recently by the Nevanlinna continuation method [23] for diagonal
entries of Green’s functions, and by the Carathéodory method [24] for both
diagonal and off-diagonal entries. In the absence of noise, these methods have
reached unprecedented accuracy in analytic continuation by resolving
complicated spectral functions with multiple features. However, these methods
also have notable drawbacks. (1) They are not numerically stable, and extended
precision arithmetic operations (128 bits of precision or higher) are needed
to carry out such calculations (even in the absence of noise). Therefore the
computational cost of the analytic continuation (especially for the
Carathéodory continuation) can be large. (2) The Nevanlinna/Carathéodory
interpolants exist if and only if Pick’s criterion [25, 26, 23] is satisfied
by the Matsubara data. In practice, noisy data often violates Pick’s
criterion, meaning that the computational result is not guaranteed to be
causal. (3) In their standard form, matrix-valued Herglotz-Nevanlinna
structure is only applicable to fermions, and hence these methods are not
directly applicable to bosonic systems.

The Herglotz-Nevanlinna function and its interpolation have been thoroughly
studied in areas such as the control theory (see e.g., [27] and [28, Chapter
9.2]). The problem of analytic continuation is also intimately related to
problems in signal processing, and in particular fitting with exponential
functions (see e.g., [29, 30, 31]). One advantage of Nevanlinna and
Carathéodory methods is that they require little prior knowledge of the
structure of Green’s functions. With further prior knowledge of the lower and
upper bound of the absolute value of the pole locations, Ying recently showed
that a modified Prony’s method can be used for performing analytic
continuation with noisy data with multiple features using only double
precision arithmetic operations [32, 33].

In this article, we propose a three-pronged Projection-Estimation-Semidefinite
relaxation (PES) method to perform analytic continuation within the causal
space, while avoiding the aforementioned drawbacks. The meaning of the three
steps are as follows: (1) Projection of the noisy data into the causal space;
(2) Estimation of pole locations using the adaptive Antoulas–Anderson (AAA)
algorithm [34]; (3) Semidefinite relaxation (SDR) fitting of Green’s functions
(diagonal and off-diagonal elements) using a bi-level optimization approach
[35]. Each step of the approach aims at resolving certain aspects of the
difficulties in the analytic continuation of noisy Matsubara Green’s
functions. The prior knowledge needed for the PES method is comparable to that
in Nevanlinna and Carathéodory methods. We demonstrate that the PES method can
robustly perform analytical continuation of noisy data using standard double
precision arithmetic operations, and is applicable to both fermionic and
bosonic systems. Table 1 compares the PES method with the Nevanlinna,
Carathéodory, MaxEnt, and Padé method. We emphasize that each of the three
steps of PES can also be useful in improving the robustness of other analytic
continuation methods. For instance, the projection step can be employed as a
data pre-processing step that significantly improves the robustness of the
Nevanlinna and Carathéodory continuation for noisy data. The semidefinite
relaxation can be combined with other continuation methods as a post-
processing refinement.

Method | $\begin{array}[]{c}\text{Noise}\\\ \text{Robustness}\end{array}$ | ${\begin{array}[]{c}\text{Calculation}\\\ \text{ Precision }\\\ \text{Requirement}\end{array}}$ | $\begin{array}[]{c}\text{Sharp}\\\ \text{Features}\end{array}$ | Causality | Applicability
---|---|---|---|---|---
This work | ✓ | Double | ✓ | ✓ | $\begin{array}[]{c}\text{Fermion}\\\ \text{Boson}\end{array}$
$\begin{array}[]{c}\text{Nevanlinna and}\\\ \text{Carath\'{e}odory}\end{array}$ | ✗ | Extended | ✓ | ${\begin{array}[]{c}\text{{\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\char 51}}\text{ if clean}\\\ \text{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\char 55}}\text{ if noisy}\end{array}}$ | Fermion
MaxEnt | ✓ | Double | ✗ | ✓ | $\begin{array}[]{c}\text{Fermion}\\\ \text{Boson}\end{array}$
Padé | ✗ | Extended | ✓ | ✗ | $\begin{array}[]{c}\text{Fermion}\\\ \text{Boson}\end{array}$
Table 1: Comparison of approaches for analytic continuation. Double means
double precision arithmetic operations, and extended means extended precision
arithmetic operations (requiring 128 bits or larger).

This article is organized as follows. In Section II, we introduce the theory
of the PES method. After giving a brief discussion on the setup of the
analytic continuation problem in Section II.1, we explain in detail each of
the three steps of the PES method: the pre-processing step using semidefinite
projection in Section II.2, the pole-estimation step using AAA algorithm in
Section II.3 and the SDR fitting step using SDR and bi-level optimization in
Section II.4. We summarize the PES method in Section II.5. The numerical
results are presented in Section III. With the Hubbard dimer example, we
compare the results of our methods with the Nevanlinna, Carathéodory and
MaxEnt methods, and show that the PES method is both noise-robust and
efficient for sharp features, and recovers both diagonal and off-diagonal
entries. We also demonstrate the necessity of the pole estimation step using
this model. In Section III.2, we show that the PES method is also applicable
to bosonic response functions and has a much better performance compared to
MaxEnt. In Section III.3, in order to demonstrate the strength of the PES
method, we conduct analytic continuations for the band structure of solids
with hundreds of orbitals.

##  II Theory

###  II.1 Analytic structure of Green’s functions

Let $\hat{c}_{i},\hat{c}_{i}^{\dagger}$ be the annihilation and creation
operator for the $i$-th orbital, $i=1,\cdots,N_{\text{orb}}$, where
$N_{\text{orb}}$ is the number of orbitals. In the Lehmann representation, a
Green’s function can be expressed as follows (see e.g., [36, Chapter 5.2]):

$\mathbb{G}_{ij}(z)=\frac{1}{Z}\sum_{r,s}\frac{\left\langle\Psi_{s}\right|\hat{c}_{i}\left|\Psi_{r}\right\rangle\left\langle\Psi_{r}\right|\hat{c}_{j}^{\dagger}\left|\Psi_{s}\right\rangle}{z+E_{s}-E_{r}}(\mathrm{e}^{-\beta
E_{s}}\mp\mathrm{e}^{-\beta E_{r}}).$ (1)

Here $\beta$ is the inverse temperature. $\left|\Psi_{s}\right\rangle$ is the
$s$-th eigenfunction of the Hamiltonian $\hat{H}$ with energy $E_{s}$, i.e.,

$\hat{H}\left|\Psi_{s}\right\rangle=E_{s}\left|\Psi_{s}\right\rangle,\quad
s=1,2,\cdots,N_{\mathrm{S}},$

where $N_{\mathrm{S}}=2^{N_{\text{orb}}}$, $\hat{H}$ is the many-body
Hamiltonian and $Z=\sum_{s=1}^{N_{\mathrm{S}}}\mathrm{e}^{-\beta E_{s}}$ is
the partition function. The negative sign corresponds to bosons, while the
positive sign corresponds to fermions.

Since $\mathbb{G}_{ij}(z)$ is defined for $i,j\in\\{1,\cdots
N_{\text{orb}}\\}$, the Green’s function $\mathbb{G}$ can be viewed as a
$N_{\text{orb}}\times N_{\text{orb}}$ matrix. In other words, we have
$\mathbb{G}(z)=\sum_{r,s}\dfrac{\mathbb{X}_{r,s}}{z-\lambda_{r,s}}$ where
$\lambda_{r,s}=E_{r}-E_{s}\in\mathbb{R}$ and
$\mathbb{X}_{r,s}=c_{r,s}v_{r,s}v_{r,s}^{\dagger}$ is a rank-1
$N_{\text{orb}}\times N_{\text{orb}}$ matrix, in which

$c_{r,s}=\frac{\mathrm{e}^{-\beta E_{s}}\mp\mathrm{e}^{-\beta E_{r}}}{Z},\quad
v_{r,s}(i)=\left\langle\Psi_{r}\right|\hat{c}_{i}^{\dagger}\left|\Psi_{s}\right\rangle.$
(2)

Here $v_{r,s}(i)$ is the $i$-th component of $v_{r,s}$. Note that for
fermionic systems, $c_{r,s}$ is always nonnegative. For bosonic systems, we
have $c_{r,s}>0$ if $E_{r}>E_{s}$, or equivalently $\lambda_{r,s}>0$; while
$c_{r,s}<0$ if $E_{r}<E_{s}$, or equivalently $\lambda_{r,s}<0$.

For orthogonal orbitals, the creation and annihilation operators satisfy the
canonical relation
$[\hat{c}_{i},\hat{c}_{j}^{\dagger}]_{\mp}=\hat{c}_{i}\hat{c}_{j}^{\dagger}\mp\hat{c}_{j}^{\dagger}\hat{c}_{i}=\delta_{ij}$.
As a result, the matrices $\mathbb{X}_{r,s}$ have the following property:

$\displaystyle\sum_{r,s}\left(\mathbb{X}_{r,s}\right)_{ij}$ (3)
$\displaystyle=$
$\displaystyle\frac{1}{Z}\sum_{r,s}\left\langle\Psi_{s}\right|\hat{c}_{i}\left|\Psi_{r}\right\rangle\left\langle\Psi_{r}\right|\hat{c}_{j}^{\dagger}\left|\Psi_{s}\right\rangle\left(\mathrm{e}^{-\beta
E_{s}}\mp\mathrm{e}^{-\beta E_{r}}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{Z}\sum_{s}\left\langle\Psi_{s}\left|\hat{c}_{i}\left(\sum_{r}\left|\Psi_{r}\right\rangle\left\langle\Psi_{r}\right|\right)\hat{c}_{j}^{\dagger}\right|\Psi_{s}\right\rangle\mathrm{e}^{-\beta
E_{s}}$ $\displaystyle\mp$
$\displaystyle\frac{1}{Z}\sum_{r}\left\langle\Psi_{r}\left|\hat{c}_{j}^{\dagger}\left(\sum_{s}\left|\Psi_{s}\right\rangle\left\langle\Psi_{s}\right|\right)\hat{c}_{i}\right|\Psi_{r}\right\rangle\mathrm{e}^{-\beta
E_{r}}$ $\displaystyle=$
$\displaystyle\frac{1}{Z}\sum_{s}\left\langle\Psi_{s}\left|\hat{c}_{i}\hat{c}_{j}^{\dagger}\mp\hat{c}_{j}^{\dagger}\hat{c}_{i}\right|\Psi_{s}\right\rangle\mathrm{e}^{-\beta
E_{s}}=\delta_{ij}.$

This indicates the sum rule
$\sum_{r,s}\mathbb{X}_{r,s}=\mathbb{I}_{N_{\text{orb}}}$, where
$\mathbb{I}_{N_{\text{orb}}}$ is the $N_{\text{orb}}\times N_{\text{orb}}$
identity matrix. For the rest of this article, let us compress the indices
$r,s$ into $l$, and let $N_{p}$ be the number of matrices $\mathbb{X}_{r,s}$
that are nonzero, i.e., the number of poles $\lambda_{r,s}$ that actually
contribute to the Green’s function. In this way, the Green’s function
$\mathbb{G}(z)$ can be written as

$\mathbb{G}(z)=\sum_{l=1}^{N_{p}}\frac{\mathbb{X}_{l}}{z-\lambda_{l}},$ (4)

where $\mathbb{X}_{l}$ satisfies the rank-1 semidefinite condition:

$\begin{array}[]{cc}\text{For
bosons:}&\begin{array}[]{cc}\operatorname{sign}(\lambda_{l})\cdot\mathbb{X}_{l}\text{
is an }N_{\text{orb}}\times N_{\text{orb}}\\\ \text{rank-1 positive
semidefinite matrix;}\end{array}\end{array}$ (5) $\begin{array}[]{cc}\text{For
fermions:}&\begin{array}[]{cc}\mathbb{X}_{l}\text{ is an }N_{\text{orb}}\times
N_{\text{orb}}\\\ \text{rank-1 positive semidefinite
matrix;}\end{array}\end{array}$ (6)

and the sum rule

$\sum_{l=1}^{N_{p}}\mathbb{X}_{l}=\mathbb{I}_{N_{\text{orb}}}.$ (7)

The function $\mathbb{G}(z)$ is defined on
$\\{0\\}\bigcup\mathbb{C}\backslash\mathbb{R}$, which excludes the real axis.
The Matsubara Green’s function $\mathbb{G}^{\mathrm{M}}(\omega)$ and the
retarded Green’s function $\mathbb{G}^{\mathrm{R}}(\omega)$ share the same
formula $\mathbb{G}(z)$ (Eq. 1 and Eq. 4) in the following way:

  * •

When $z=\mathrm{i}\omega_{n}\in\mathrm{i}\mathbb{R}$ is on the imaginary axis,
$\mathbb{G}^{\mathrm{M}}(\omega_{n})=\mathbb{G}(\mathrm{i}\omega_{n})$ is the
Matsubara (or imaginary frequency) Green’s function, $\omega_{n}$ is called
the Matsubara frequency. For fermionic systems,
$\mathrm{i}\omega_{n}\in\mathrm{i}\frac{2\mathbb{Z}+1}{\beta}$; for bosonic
systems, $\mathrm{i}\omega_{n}\in\mathrm{i}\frac{2\mathbb{Z}}{\beta}$.

  * •

Since $\mathbb{G}(z)$ has poles on the real axis, the real-time (or real-
frequency) Green’s function could only be evaluated at positions
infinitesimally close to the real axis. Let $\eta$ be a positive infinitesimal
number. When $z=\omega+\mathrm{i}\eta$,
$\mathbb{G}^{\mathrm{R}}(\omega)=\mathbb{G}(\omega+\mathrm{i}\eta)$ is the
retarded Green’s function.

The spectral function $A(\omega)$ is defined from
$\mathbb{G}^{\mathrm{R}}(\omega)$ as follows:

$A(\omega)=-\frac{1}{\pi}\operatorname{Im}\left({\operatorname{Tr}\left(\mathbb{G}^{\mathrm{R}}(\omega)\right)}\right)$
(8)

The spectral function contains information of the excitation spectra in
quantum systems. In practice, the retarded Green’s function and the spectral
function are evaluated with $\eta$ a very small positive number.

Green’s functions of fermionic systems are closely related to the matrix-
valued Herglotz-Nevanlinna functions (see e.g., [22]). A matrix-valued
function $\mathbb{N}(z)$ is said to be Herglotz-Nevanlinna if
$\mathbb{N}:\mathbb{C}_{+}\rightarrow\mathbb{C}^{N_{\text{orb}}\times
N_{\text{orb}}}$ is analytic and $\operatorname{Im}(\mathbb{N}(z))$ is a
positive semidefinite matrix for $z\in\mathbb{C}_{+}$. Here $\mathbb{C}_{+}$
is the open complex upper half-plane, $\mathbb{C}^{N_{\text{orb}}\times
N_{\text{orb}}}$ is the set of $N_{\text{orb}}\times N_{\text{orb}}$ matrices
with entries in $\mathbb{C}$, and $\operatorname{Im}(\mathbb{N}(z))$ is the
imaginary part of $\mathbb{N}(z)$. The matrix-valued Herglotz-Nevanlinna
functions admit the following integral representation, (see e.g., [37, Theorem
5.4] for proof):

$\mathbb{N}(z)=Mz+B+\int_{\mathbb{R}}\left(\frac{1}{t-z}-\frac{t}{1+t^{2}}\right)\mathrm{d}\Sigma(t),\text{
}z\in\mathbb{C}_{+},$ (9)

where $M,B\in\mathbb{C}^{N_{\text{orb}}\times N_{\text{orb}}}$, $M$ is
positive semidefinite, $B$ is Hermitian, and $\Sigma(t)$ is a nondecreasing
matrix-valued function on $\mathbb{R}$ such that
$\int_{\mathbb{R}}\frac{u^{\dagger}(\mathrm{d}\Sigma(t))u}{1+t^{2}}<\infty$
for any vector $u\in\mathbb{C}^{N_{\text{orb}}}$. This condition on
$\Sigma(t)$ ensures that the above integration in Eq. 9 is well-defined. For
the detailed mathematical theory of Herglotz-Nevanlinna functions, we refer
the readers to [38, Chapter 2] and [39, Chapter 3].

From Eq. 9, we may set

$\mathrm{d}\Sigma(t)=\sum_{l=1}^{N_{p}}\mathbb{X}_{l}\delta(t-\lambda_{l}),\quad
M=0,\quad B=\int_{\mathbb{R}}\frac{t}{1+t^{2}}\mathrm{d}\Sigma(t),$ (10)

then

$\mathbb{N}(z)=\int_{\mathbb{R}}\frac{\mathrm{d}\Sigma(t)}{t-z}=\sum_{l=1}^{N_{p}}\frac{\mathbb{X}_{l}}{\lambda_{l}-z}=-\mathbb{G}(z),\quad
z\in\mathbb{C}_{+}.$ (11)

In other words, $-\mathbb{G}(z)$ is a matrix-valued Herglotz-Nevanlinna
function. This is the mathematical foundation of Nevanlinna [23] and
Carathéodory continuation [24] methods.

The causal space for Green’s functions of fermionic systems is a subset of the
Herglotz-Nevanlinna function class denoted by $\mathcal{S}_{\mathrm{F}}$:

$\mathcal{S}_{\mathrm{F}}=\left\\{\mathbb{G}\left|\begin{array}[]{c}\mathbb{G}(z)=\sum_{l=1}^{N_{p}}\frac{\mathbb{X}_{l}}{z-\lambda_{l}},\text{
for some }N_{p}\in\mathbb{Z},\\\ \text{ for some
}\lambda_{l}\in\mathbb{R},\text{ and for }\mathbb{X}_{l}\text{ satisfying}\\\
\text{ fermionic rank-1 semidefinite}\text{ conditions}\\\ \text{
\lx@cref{creftype~refnum}{eq:fermionSDconditions} and}\text{ the sum rule
\lx@cref{creftype~refnum}{eq:normcondition}}.\end{array}\right.\right\\}.$
(12)

Similarly, the causal space of Green’s functions of bosonic systems is

$\mathcal{S}_{\mathrm{B}}=\left\\{\mathbb{G}\left|\begin{array}[]{c}\mathbb{G}(z)=\sum_{l=1}^{N_{p}}\frac{\mathbb{X}_{l}}{z-\lambda_{l}},\text{
for some }N_{p}\in\mathbb{Z},\\\ \text{ for some
}\lambda_{l}\in\mathbb{R},\text{ and for }\mathbb{X}_{l}\text{ satisfying}\\\
\text{ fermionic rank-1 semidefinite}\text{ conditions}\\\ \text{
\lx@cref{creftype~refnum}{eq:bosonSDconditions} and}\text{ the sum rule
\lx@cref{creftype~refnum}{eq:normcondition}}.\end{array}\right.\right\\}.$
(13)

In the definition of $S_{\pm}$, the rank-1 condition in $\mathbb{X}_{l}$ can
be dropped, because any semi-definite matrix can be represented as the sum of
several rank-1 semi-definite matrices. Degenerate excitations can be treated
similarly by allowing the values of some $\lambda_{l}$’s to be the same. This
is the mathematical foundation of the semidefinite relaxation fitting [35].

Now we are ready to introduce the setup for analytic continuation problems of
Matsubara data. Given several (possibly noisy) Matsubara data
$\mathbb{G}_{n}\approx\mathbb{G}(\mathrm{i}\omega_{n})$ for
$n=1,2,\cdots,N_{w}$, our goal is to obtain the fitting of these data into the
following analytic form
$\mathbb{G}(z)=\sum_{l=1}^{N_{p}}\frac{\mathbb{X}_{l}}{z-\lambda_{l}}$, where
$\lambda_{l}\in\mathbb{R}$ and $\mathbb{X}_{l}$ is a $N_{\text{orb}}\times
N_{\text{orb}}$ rank-1 semidefinite matrix, and the matrices
$\left\\{\mathbb{X}_{l}\right\\}_{l=1}^{N_{p}}$ satisfies the sum rule
$\sum_{l}\mathbb{X}_{l}=\mathbb{I}_{N_{\text{orb}}}$. For fermionic systems,
$\mathbb{X}_{l}$ is positive semidefinite, while for bosonic cases,
$\operatorname{sign}(\lambda_{l})\cdot\mathbb{X}_{l}$ is positive
semidefinite.

In summary, we would like to obtain the poles $\lambda_{l}$, the semidefinite
matrices $\mathbb{X}_{l}$, and the number of poles $N_{p}$ which we may not
have a priori knowledge. From such information we can infer the spectral
function, as well as other diagonal and off-diagonal entries of the Green’s
function.

Finally, we would like to remark that the Green’s functions are only one type
of response functions. Other response functions commonly considered, such as
the charge, magnetic, or superconducting susceptibilities (see e.g., [36,
Chapter 5.5]) admit the same formula Eq. 4, in which $\mathbb{X}_{l}$ also
satisfies the bosonic/fermionic semidefinite conditions (see Eqs. 5 and 6).
The only difference is that $\mathbb{X}_{l}$ may satisfy a different set of
sum rules. In other words, the analytic structure of Green’s function and
other response functions are similar, therefore the same methodology of
analytic continuation applies.

###  II.2 Step 1: Projection of the noisy data to the causal space

The Matsubara data often contains unphysical noise, i.e.,
$\mathbb{G}^{\mathrm{M}}(\mathrm{i}\omega_{n})$ cannot be expressed as a
matrix-value function in the set $\mathcal{S}_{\mathrm{F}/\mathrm{B}}$.
Therefore the first step of our algorithm is to project the noisy data into
the causal space. This can be achieved efficiently using semidefinite
programming. For simplicity, let us choose a fine uniform grid on the real-
axis

$x_{m}=-\Lambda_{x}+mh_{x},\quad h_{x}=\frac{2C}{M},\quad m=0,\cdots,M,$ (14)

where $\Lambda_{x}$ is the real-axis cutoff, and $h_{x}$ is the grid size. Our
goal is to fit the Matsubara data $\mathbb{G}_{n}$ via the following form:

$\mathbb{G}_{n}\approx\sum_{m=0}^{M}\frac{\mathbb{P}_{m}}{\mathrm{i}\omega_{n}-x_{m}},\quad
n=1,\cdots N_{w},$ (15)

where $\\{\mathbb{P}_{m}\\}$ are semidefinite matrices. The objective function
$\mathcal{E}_{\mathrm{proj}}\left(\\{\mathbb{P}_{m}\\}_{m=0}^{M}\right)$ is
defined as

$\mathcal{E}_{\mathrm{proj}}\left(\\{\mathbb{P}_{m}\\}_{m=0}^{M}\right)=\left(\sum_{n=1}^{N_{w}}\left\|\mathbb{G}_{n}-\sum_{m=0}^{M}\frac{\mathbb{P}_{m}}{\mathrm{i}\omega_{n}-x_{m}}\right\|^{2}_{F}\right)^{1/2}.$
(16)

Here $\left\lVert\cdot\right\rVert_{F}$ is the Frobenius norm. The solutions
of the projection $\mathbb{P}_{m}^{\mathrm{proj}}$ is obtained by solving the
following optimization problem:

$\begin{array}[]{c}\\{\mathbb{P}_{m}^{\mathrm{proj}}\\}=\arg\min_{\\{\mathbb{P}_{m}\\}}\mathcal{E}_{\mathrm{proj}}\left(\\{\mathbb{P}_{m}\\}_{m=1}^{M}\right)\\\
\text{ subject to }\left\\{\begin{array}[]{l}\text{i) semidefinite constraint:
}\\\ \quad\left\\{\begin{array}[]{rr}\mathbb{P}_{m}\succeq
0&\text{(Fermions)},\\\ \operatorname{sign}(x_{m})\cdot\mathbb{P}_{m}\succeq
0&\text{(Bosons case)};\end{array}\right.\\\ \text{ ii) sum rules: }\\\
\quad\sum_{m}\mathbb{P}_{m}=\mathbb{I}_{N_{\text{orb}}}.\end{array}\right.\end{array}$
(17)

This is a convex optimization problem, and can be solved efficiently using
software packages such as CVX [40]. If we are performing analytic continuation
of other correlation functions, $\\{\mathbb{P}_{m}\\}_{m=1}^{M}$ may be
subject to a different sum rule.

For scalar-valued fermionic Green’s function, the positive semidefinite
condition $\mathbb{P}_{m}\succeq 0$ becomes the non-negativity condition
$\mathbb{P}_{m}\geq 0$. Therefore the sum rule $\sum_{m}\mathbb{P}_{m}=1$ can
also be written as $\sum_{m}|\mathbb{P}_{m}|=1$, i.e., an $\ell^{1}$-norm
constraint on the vector $(\mathbb{P}_{1},\cdots,\mathbb{P}_{m})$. The least
squares problem with a $\ell^{1}$-norm constraint is similar to the well-known
Least Absolute Selection and Shrinkage Operator (LASSO) problem [41], which
favors solutions with a sparse structure (see also [42]). This agrees with our
numerical observation that the solution $\mathbb{P}_{m}$ often has relatively
few nonzero entries. For matrix-valued Green’s function, there is a similar
mechanism that induces the sparsity of solution. The natural generalization of
the $\ell^{1}$ norm is the nuclear norm $\|\cdot\|_{*}$, defined as
$\|A\|_{*}=\sqrt{\operatorname{Tr}(A^{\dagger}A)}=\sum_{k}\sigma_{k}(A)$,
where $\sigma_{k}(A)$ is the $k$-th largest singular value for $A$. Note that
each $\mathbb{P}_{m}$ is a positive semidefinite matrix, therefore the $k$-th
largest singular value $\sigma_{k}(\mathbb{P}_{m})$ is equal to the $k$-th
largest eigenvalue $\lambda_{k}(\mathbb{P}_{m})$. Therefore, we have
$\sum_{m}\|\mathbb{P}_{m}\|_{*}=\sum_{m}\sum_{k}\sigma_{k}(\mathbb{P}_{m})=\sum_{m}\sum_{k}\lambda_{k}(\mathbb{P}_{m})=\sum_{m}\operatorname{Tr}(\mathbb{P}_{m})=\operatorname{Tr}(\sum_{m}\mathbb{P}_{m})$.
Combined with the sum rule
$\sum_{m}\mathbb{P}_{m}=\mathbb{I}_{N_{\operatorname{orb}}}$, we can see that
the matrices $\\{\mathbb{P}_{m}\\}$’s are enforced to satisfy the nuclear norm
constraint $\sum_{m}\|\mathbb{P}_{m}\|_{*}=N_{\operatorname{orb}}$. The
nuclear norm is also a sparsity-inducing property (see e.g., [43]).

After obtaining $\mathbb{P}_{m}^{\operatorname{proj}}$, we can construct the
projected Matsubara data $\mathbb{G}_{n}^{\mathrm{proj}}$:

$\mathbb{G}_{n}^{\mathrm{proj}}=\sum_{m=0}^{M}\frac{\mathbb{P}_{m}^{\mathrm{proj}}}{\mathrm{i}\omega_{n}-x_{m}},\quad
n=1,\cdots N_{w},$ (18)

This projection step can be used to improve the robustness of other analytic
continuation methods, and therefore is of independent interest. For instance,
the Nevanlinna analytic continuation method requires the Pick matrix to be
positive semidefinite (see [23] for details). Numerical results indicate that
this criterion can be violated when very small perturbations are applied to
the Matsubara data. Given that noise is inevitable in many Green’s function
solvers (notably, quantum Monte Carlo solvers), the application range of the
Nevanlinna analytic continuation is thus significantly limited by the nature
of the noise. Our numerical results suggest that when the Nevanlinna method is
applied to the projected noisy Matsubara data (for diagonal entries), the
quality of the analytic continuation is significantly improved. (See Section
III.1 and Section III.3 for details.) Similar improvements are observed on
both diagonal and off-diagonal entries of the Green’s functions while applying
Carathéodory continuations on the projected data (also see Section III.1 for
details).

###  II.3 Step 2: Estimation of pole locations using the AAA algorithm

While the projection step can be viewed as an analytic continuation algorithm
by itself, the quality of the continuation is constrained by the resolution of
the grid on the real axis. In practice, the number of grid points (i.e., $M$)
in the uniform grid is often too small to accurately resolve the pole
locations, but is too large to be directly used in the subsequent semidefinite
relaxation (SDR) step to be detailed in Section II.4. Furthermore, the loss
function in the SDR step is highly non-convex, and the optimization with
respect to this loss function requires a proper initial guess. These
considerations lead to the second step of the algorithm for estimating the
locations of a relatively small number of poles.

We use the adaptive Antoulas–Anderson (AAA) algorithm [34] for the pole
estimation, which is available as a subroutine in the Chebfun package [44].
The AAA algorithm is able to obtain the poles of a scalar complex-valued
function $g(z)$. In the context of analytic continuation, the scalar function
$g(z)$ could either be each diagonal entry
$\mathbb{G}_{ii}(z)(i=1,\cdots,N_{\text{orb}})$ of the Green’s function, or
the trace $\operatorname{Tr}\left(\mathbb{G}(z)\right)$. In practice, we find
that performing the pole estimation on each diagonal entry separately often
gives better results.

The AAA algorithm performs a rational approximation to a scalar function
$g(z)$ using barycentric representations [45]:

$g(z)\approx\frac{n(z)}{d(z)}=\frac{\sum_{j=1}^{n}\frac{w_{j}g_{j}}{z-z_{j}}}{\sum_{j=1}^{n}\frac{w_{j}}{z-z_{j}}}.$
(19)

Here $n(z)=\sum_{j=1}^{n}\dfrac{w_{j}g_{j}}{z-z_{j}}$ and
$d(z)=\sum_{j=1}^{n}\dfrac{w_{j}}{z-z_{j}}$ are the numerator and denominator
part of the barycentric representation, respectively. Note that the above
formulation Eq. 19 automatically satisfies that
$\frac{n(z_{j})}{d(z_{j})}=g_{j}$.

We briefly describe the AAA algorithm below, and refer readers to Ref. [34]
for more details. Assume that we are given the values of $g(z)$ on a set
$\mathcal{Z}\subset\mathbb{C}$ of $N$ points (in analytic continuation,
$\mathcal{Z}=\\{\mathrm{i}\omega_{n},n=1,\cdots,N_{w}\\}$), the AAA algorithm
aims at finding a set of support points $\\{z_{j}\\}_{i=1}^{n}$, which allows
us to solve for the weights $\\{w_{j}\\}_{j=1}^{n}$ in Eq. 19. Here
$g_{j}=g(z_{j})$ ($j=1,\cdots,n$). Starting from $n=1$, the AAA algorithm
gradually expands the set $\\{z_{j}\\}_{i=1}^{n}$ following a greedy
algorithm.

At the $(n-1)$-th step ($n=1,2,\cdots$), if the residual
$\mathrm{res}(z)=g(z)-\frac{n(z)}{d(z)}$ is sufficiently small for all
$z\in\mathcal{Z}$, then the algorithm terminates. Otherwise, we pick the next
support point $z_{n}$, by choosing
$z\in\mathcal{Z}\backslash\left\\{z_{1},\ldots,z_{n-1}\right\\}$ where the
residual $\mathrm{res}(z)$ takes its maximum absolute value. Let us write
$\mathcal{Z}\backslash\left\\{z_{1},\cdots,z_{n}\right\\}$ as a column vector
${\bm{Z}}^{(n)}:=\left(Z_{1}^{(n)},\cdots,Z_{N-n}^{(n)}\right)^{T}$ and define
${\bm{g}}^{(n)}=\left(g_{1}^{(n)},\cdots,g_{N-n}^{(n)}\right)^{T}:=\left(g\left(Z_{1}^{(n)}\right),\cdots,g\left(Z_{N-n}^{(n)}\right)\right)^{T}$.
Since $g(z)$ is supposed to be approximated by $g(z)\approx n(z)/d(z)$, i.e.
$g(z)d(z)\approx n(z)$, we choose the normalized weight
$w=(w_{1},\cdots,w_{n})$ by solving the least square problem:

$\displaystyle\min_{\|w\|_{2}=1}$
$\displaystyle\sum_{i=1}^{N-n}\left|g_{i}^{(n)}d\left(Z_{i}^{(n)}\right)-n\left(Z_{i}^{(n)}\right)\right|^{2}$
(20) $\displaystyle=$
$\displaystyle\sum_{i=1}^{N-n}\left|\sum_{j=1}^{n}\frac{w_{j}\left(g_{i}^{(n)}-g_{j}\right)}{Z_{i}^{(n)}-z_{j}}\right|^{2}.$

The minimization problem of Eq. 20 could be viewed as the following linear
algebra problem:

$\min_{\|w\|_{2}=1}\left\|A^{(n)}w\right\|_{2},$ (21)

where $A^{(n)}$ is

$\displaystyle
A^{(n)}=\left(\begin{array}[]{ccc}\frac{g_{1}^{(n)}-g_{1}}{Z_{1}^{(n)}-z_{1}}&\cdots&\frac{g_{1}^{(n)}-g_{n}}{Z_{1}^{(n)}-z_{n}}\\\
\vdots&\ddots&\vdots\\\
\frac{g_{N-n}^{(n)}-g_{1}}{Z_{N-n}^{(n)}-z_{1}}&\cdots&\frac{g_{N-n}^{(n)}-g_{n}}{Z_{N-n}^{(n)}-z_{n}}\end{array}\right).$
(22)

This optimization problem Eq. 21 could be solved by performing a singular
value decomposition (SVD) on the matrix $A^{(n)}=U\Sigma V^{\dagger}$, where
$U$ is a $(N-n)\times r$ orthogonal matrix, $V$ is a $n\times r$ orthogonal
matrix, $\Sigma=\operatorname{diag}(\sigma_{1},\cdots,\sigma_{r})$ is a
diagonal matrix ($\sigma_{1}\geq\cdots\geq\sigma_{r}$), and $r\leq\min(N-n,n)$
is the rank of $A^{(n)}$. The solution $w$ of Eq. 21 should be taken as the
last column of $V$.

After the iteration terminates at the $K$-th step, we have
$\\{z_{1},\cdots,z_{K}\\}$ and its associated weights
$\\{w_{1},\cdots,w_{K}\\}$. We can calculate the zeros of $d(z)$, which serve
as the estimated poles of $g(z)$ that we want to obtain, using the following
generalized eigenvalue problem [34]:

$\left(\begin{array}[]{ccccc}0&w_{1}&w_{2}&\cdots&w_{K}\\\ 1&z_{1}&&&\\\
1&&z_{2}&&\\\ \vdots&&&\ddots&\\\
1&&&&z_{K}\end{array}\right)=\xi\left(\begin{array}[]{ccccc}0&&&&\\\ &1&&&\\\
&&1&&\\\ &&&\ddots&\\\ &&&&1\end{array}\right).$ (23)

At least two of the eigenvalues in this problem are infinite, and the
remaining $(K-1)$ eigenvalues $\\{\xi_{a}\\}_{a=1}^{K-1}$ are the zeros of
$d$, i.e. the complex-valued poles of $g$.

Since all poles of the Green’s function in Eq. 4 are real-valued, after
obtaining $\\{\xi_{a}\\}_{a=1}^{K-1}$, we discard the poles far away from the
real axis. In other words, we only keep poles with
$\left\lvert\operatorname{Im}\xi_{a}\right\rvert\leq\varepsilon_{p}$ for some
$\varepsilon_{p}>0$, and define
$\\{\lambda_{l}^{\operatorname{in}}\\}_{l=1}^{N_{p}}$ to be real parts of the
remaining poles. These $N_{p}$ poles will be used as the initial guess of the
poles in the semidefinite relaxation step below.

###  II.4 Step 3: Semidefinite relaxation

In the final step, we use numerical optimization to obtain an effective
fitting of the Matsubara data in the form of Eq. 4. The fitting error is
defined as

$\displaystyle\text{Err}\left(\\{\lambda_{l}\\}_{l=1}^{N_{p}},\\{\mathbb{X}_{l}\\}_{l=1}^{N_{p}}\right)$
(24) $\displaystyle=$
$\displaystyle\left(\sum_{n=1}^{N}\left\|{\mathbb{G}}_{n}-\mathbb{G}\left(\mathrm{i}\omega_{n};\\{\lambda_{l},\mathbb{X}_{l}\\}_{l=1}^{N_{p}}\right)\right\|_{F}^{2}\right)^{1/2}$
$\displaystyle=$
$\displaystyle\left(\sum_{n=1}^{N}\left\|\mathbb{G}_{n}-\sum_{l=1}^{N_{p}}\frac{\mathbb{X}_{l}}{\mathrm{i}\omega_{n}-\lambda_{l}}\right\|_{F}^{2}\right)^{1/2}.$

The fitting error expressed in Eq. 24 is a highly nonconvex function, and the
minimization of this function can frequently be trapped in local minima, which
strongly depend on the choice of the initial guess.

Recall that in Eq. 4 $\\{\mathbb{X}_{l}\\}$ are required to be rank-1
semidefinite matrix. The semidefinite relaxation drops this rank-1 constraint.
In other words, $\\{\mathbb{X}_{l}\\}$ (or
$\operatorname{sign}(\lambda_{l})\cdot\mathbb{X}_{l}$) are only required to be
positive semidefinite matrices for the analytic continuation of fermionic
(bosonic) Green’s function.

For simplicity, the following discussion focuses on fermionic Green’s
functions, while the same numerical treatment is applicable to bosonic Green’s
functions. The minimization of
$\text{Err}\left(\\{\lambda_{l}\\}_{l=1}^{N_{p}},\\{\mathbb{X}_{l}\succeq
0\\}_{l=1}^{N_{p}}\right)$ can be formulated as a bi-level optimization
problem:

$\begin{aligned} &\min_{\lambda_{l},\mathbb{X}_{l}\succeq
0}\text{Err}\left(\\{\lambda_{l}\\}_{l=1}^{N_{p}},\\{\mathbb{X}_{l}\\}_{l=1}^{N_{p}}\right)\\\
=&\min_{\lambda_{l}\in\mathbb{R}}\min_{\mathbb{X}_{l}\succeq
0}\text{Err}\left(\\{\lambda_{l}\\}_{l=1}^{N_{p}},\\{\mathbb{X}_{l}\\}_{l=1}^{N_{p}}\right)\\\
=&\min_{\lambda_{l}\in\mathbb{R}}\mathcal{E}\left(\\{\lambda_{l}\\}_{l=1}^{N_{p}}\right)\end{aligned},$
(25)

where

$\mathcal{E}\left(\\{\lambda_{l}\\}_{l=1}^{N_{p}}\right)=\min_{\mathbb{X}_{l}\succeq
0}\text{Err}\left(\\{\lambda_{l}\\}_{l=1}^{N_{p}},\\{\mathbb{X}_{l}\succeq
0\\}_{l=1}^{N_{p}}\right).$ (26)

When the poles $\\{\lambda_{l}\\}_{l=1}^{N_{p}}$ are fixed, the optimization
with respect to semidefinite matrices $\\{\mathbb{X}_{l}\succeq 0\\}$ is a
convex optimization problem, and could be efficiently solved using software
packages such as CVX [40]. The optimization of
$\text{Err}\left(\\{\lambda_{l}\\}_{l=1}^{N_{p}},\\{\mathbb{X}_{l}\succeq
0\\}_{l=1}^{N_{p}}\right)$ is transformed into the optimization of
$\mathcal{E}\left(\\{\lambda_{l}\\}_{l=1}^{N_{p}}\right)$, but note that such
an optimization over the pole positions is still a nonconvex problem.

Therefore, we have the following SDR fitting procedure through the bi-level
optimization framework: starting from an initial value for poles
$\\{\lambda_{l}^{\operatorname{in}}\\}_{l=1}^{N_{p}}$ obtained by the previous
pole estimation step, we use numerical optimization technique to minimize
$\mathcal{E}\left(\\{\lambda_{l}\\}_{l=1}^{N_{p}}\right)$, where

  * •

The value of $\mathcal{E}\left(\\{\lambda_{l}\\}_{l=1}^{N_{p}}\right)$ and the
corresponding optimal $\\{\mathbb{X}_{l}^{\text{opt}}\succeq 0\\}$ for fixed
poles $\\{\lambda_{l}\\}_{l=1}^{N_{p}}$ are found by a convex optimization
solver.

  * •

Since
$\mathcal{E}\left(\\{\lambda_{l}\\}_{l=1}^{N_{p}}\right)=\text{Err}\left(\\{\lambda_{l}\\}_{l=1}^{N_{p}},\\{\mathbb{X}_{l}^{\text{opt}}\\}_{l=1}^{N_{p}}\right)$
and the optimal matrices $\\{\mathbb{X}_{l}^{\text{opt}}\\}$ should be
regarded as functions of the poles $\\{\lambda_{l}\\}$, the partial derivative
of $\mathcal{E}\left(\\{\lambda_{l}\\}_{l=1}^{N_{p}}\right)$ could be written
as

$\begin{aligned}
\partial_{m}\mathcal{E}&=\frac{\partial}{\partial\lambda_{m}}\text{Err}\\\
&+\sum_{l=1}^{N_{p}}\sum_{i,j=1}^{N_{\text{orb}}}\left.\frac{\partial\text{
Err}}{\partial\left(\mathbb{X}_{l}\right)_{ij}}\right|_{\mathbb{X}_{l}^{\text{opt}}}\cdot\frac{\partial\left(\mathbb{X}_{l}^{\text{opt}}\right)_{ij}}{\partial\lambda_{m}}\end{aligned},$
(27)

The first order optimality condition for solving Eq. 26 implies
$\left.\frac{\partial\text{
Err}}{\partial\left(\mathbb{X}_{l}\right)_{ij}}\right|_{\mathbb{X}_{l}^{\text{opt}}}=0$.
Such an evaluation of the derivative quantities is similar to the treatment in
the Hellmann-Feynman’s theorem [46]. Therefore we have

$\displaystyle\partial_{m}$
$\displaystyle\mathcal{E}\left(\\{\lambda_{l}\\}_{l=1}^{N_{p}}\right)=\frac{\partial}{\partial\lambda_{m}}\text{Err}\left(\\{\lambda_{l}\\}_{l=1}^{N_{p}},\\{\mathbb{X}_{l}\\}_{l=1}^{N_{p}}\right)\Big{|}_{\mathbb{X}_{l}=\mathbb{X}_{l}^{\text{opt}}}$
(28) $\displaystyle=$
$\displaystyle-\frac{\sum\limits_{n=1}^{N}\sum\limits_{i,j=1}^{N_{\text{orb}}}\operatorname{Re}\left(\left((\mathbb{G}_{n})_{ij}-\sum\limits_{l}\frac{(\mathbb{X}_{l}^{\text{opt}})_{ij}}{\mathrm{i}\omega_{n}-\lambda_{l}}\right)\frac{(\mathbb{X}_{m}^{\text{opt}})_{ij}^{*}}{(\mathrm{i}\omega_{n}+\lambda_{m})^{2}}\right)}{\sqrt{\mathcal{E}\left(\\{\lambda_{l}\\}_{l=1}^{N_{p}}\right)}}.$

Therefore, both the function value and gradient of $\mathcal{E}$ can be
efficiently evaluated. The optimization of
$\mathcal{E}\left(\\{\lambda_{l}\\}_{l=1}^{N_{p}}\right)$ could be conducted
with standard gradient-based optimization solvers, such as the Broyden-
Fletcher-Goldfarb-Shanno (BFGS) method (see e.g. [47]).

Note that the SDR fitting step enforces that the Green’s functions must be in
the causal space, and is thus robust to unphysical noises. Since the
projection step may introduce biases especially when the uniform mesh in Eq.
14 is relatively coarse, we find that the performance of SDR fitting is often
better when applied to the original noisy data, instead of the data obtained
after the projection.

###  II.5 Summary of the PES method

We summarize the PES method in Algorithm 1 below.

Algorithm 1 Projection-Estimation-Semidefinite relaxation (PES) method

1: Parameters: real axis cutoff $\Lambda_{x}$ and mesh size $h_{x}$ in the
projection step; imaginary axis cutoff $\varepsilon_{p}$ in the pole
estimation step.

2: Input: Matsubara data $\left\\{\mathbb{G}_{n}\right\\}_{n=1}^{N_{w}}$,
where $\mathbb{G}_{n}$ is an $N_{\text{orb}}\times N_{\text{orb}}$ matrix.

3: Output: Poles $\left\\{\lambda_{l}\right\\}^{N_{p}}_{l=1}$, weights
$\\{\mathbb{X}_{l}\\}^{N_{p}}_{l=1}$.

4: Projection of the noisy data
$\left\\{\mathbb{G}_{n}\right\\}_{n=1}^{N_{w}}$: form a uniform mesh in Eq.
14, minimize
$\mathcal{E}_{\mathrm{proj}}\left(\\{\mathbb{P}_{m}\\}_{m=1}^{M}\right)$ in
Eq. 17, and then obtain the projected data
$\left\\{\mathbb{G}_{n}^{\text{proj}}\right\\}_{n=1}^{N_{w}}$ using Eq. 18.

5: Estimation of the pole locations using the projected data
$\left\\{\mathbb{G}_{n}^{\text{proj}}\right\\}_{n=1}^{N_{w}}$ via the AAA
algorithm. The output is
$\left\\{\lambda_{l}^{\text{in}}\right\\}_{l=1}^{N_{p}}$.

6: Semidefinite relaxation fitting through bi-level optimization: with the
initial poles $\left\\{\lambda_{l}^{\text{in}}\right\\}_{l=1}^{N_{p}}$, we can
conduct SDR fitting by minimizing
$\mathcal{E}\left(\\{\lambda_{l}\\}_{l=1}^{N_{p}}\right)$ in Eq. 26 to obtain
the final poles $\left\\{\lambda_{l}\right\\}_{l=1}^{N_{p}}$ and weights
$\left\\{\mathbb{X}_{l}\right\\}_{l=1}^{N_{p}}$.

Computational cost: Both the projection and the semidefinite relaxation steps
require the solution of convex optimization problems that can be reformulated
into semidefinite programming (SDP) problems. The cost for solving the convex
optimization problems can depend on the SDP reformulation, and the choice of
algorithms. The default solver in the CVX software package is SDTP3 [48],
which implements the interior point method with a primal-dual path-following
strategy [49] to solve SDP problems. In the primal-dual interior point method,
the cost is dominated by solving the Newton equation to obtain the primal-dual
search direction. The unknown variables are $N_{v}$ positive semidefinite
matrices each of size $N_{d}\times N_{d}$, and the total number of variables
is $N_{v}N_{d}^{2}$. Therefore the cost of solving Newton’s equation can be as
large as $\mathcal{O}(N_{v}^{3}N_{d}^{6})$, and there is an additional cost in
forming the matrices in Newton’s equation (see e.g., [50]). While the cost of
solving Newton’s equation may be reduced using preconditioned and iterative
solvers in certain cases [51], in general, the solution the SDP problem can
become expensive when both $N_{v}$ and $N_{d}$ are large.

In the projection step, $N_{v}=M$ is usually large. Note that the pole
estimation step only uses the diagonal entries of the Green’s function, and it
is sufficient to only conduct projection of these diagonal entries. If $M$ is
very large, we may further reduce the cost by performing projection for each
diagonal entry of the Green’s functions separately. The semidefinite
relaxation fitting requires the solution of $N_{v}=N_{p}$ matrices each of
size $N_{d}=N_{\text{orb}}$. Therefore it is important for the pole estimation
step to obtain a relatively small number of poles. Furthermore, if the non-
diagonal information of the Green’s function is not needed (for example if the
goal is to calculate the spectral function), the cost can be significantly
reduced by conducting the semidefinite relaxation step for the diagonal
entries only. Compared to solving the SDP problems, the cost of the AAA
algorithm for the pole estimation is usually negligible.

Noiseless data: If the data is noiseless or if the Green’s functions are
causal by construction (e.g., Green’s functions are approximated via causal
diagrammatic methods), we can skip the projection step, and conduct estimation
and SDR fitting directly. Furthermore, since the projection is only performed
on a finite grid along the real axis, the pre-processing step may become
detrimental to the accuracy when the grid size is relatively small.

Parameters in AAA: Numerical results suggest that the final result is often
insensitive to the choice of the parameter $\varepsilon_{p}$. For noisy data,
the Lawson refinement step [52] (also implemented in Chebfun), which is based
on an iteratively reweighted least-squares (IRLS) iteration can also be used
to improve the fitting result. Also, the AAA algorithm in Chebfun [44] could
be implemented with or without specifying the number of poles. Prior knowledge
on the on the number of peaks in the physical system can be used to control
the maximal number of poles obtained by the AAA algorithm, and to enhance the
robustness with respect to noise.

##  III Results

In order to demonstrate the strengths of the PES method and give a fair
comparison with other methods, we present various results using the Hubbard
dimer system and the band structure calculation of different materials. When
testing the performances of analytic continuation methods in the presence of
noise, we will manually add multiplicative Gaussian noise to the clean
Matsubara data:

$\mathbb{G}_{\text{noisy }}=\mathbb{G}_{\text{exact }}\cdot\left(1+\sigma
N_{\mathbb{C}}(0,1)\right)$

Here $N_{\mathbb{C}}(0,1)$ is the complex-valued normal Gaussian noise and
$\sigma$ is referred to as the (multiplicative) noise level of the data.

###  III.1 Hubbard dimer

####  III.1.1 Setup of the Hubbard Dimer

The Hubbard dimer example used in [24] is a simple but nontrivial test case,
since it has multiple sharp features while the true value of the Green’s
function could be obtained by the exact diagonalization. It is also a
prototypical system for studying the excitation spectra of molecules. There
are only four orbitals in the Hubbard dimer, namely
$|0\uparrow\rangle,|0\downarrow\rangle,|1\uparrow\rangle,|1\downarrow\rangle$.
The Hamiltonian is

$\hat{H}=\hat{H}_{0}+\hat{H}_{1},$ (29)

where

$\displaystyle\hat{H}_{0}$
$\displaystyle=-\sum_{\sigma\in\\{\uparrow,\downarrow\\}}t\left(\hat{c}_{0\sigma}^{\dagger}\hat{c}_{1\sigma}+\hat{c}_{1\sigma}^{\dagger}\hat{c}_{0\sigma}\right)$
(30) $\displaystyle-$
$\displaystyle\mu\left(\hat{n}_{0,\uparrow}+\hat{n}_{0,\downarrow}+\hat{n}_{1,\uparrow}+\hat{n}_{1,\downarrow}\right),$
$\displaystyle\hat{H}_{1}$
$\displaystyle=(U+U_{a})\hat{n}_{0,\uparrow}\hat{n}_{0,\downarrow}+(U-U_{a})\hat{n}_{1,\uparrow}\hat{n}_{1,\downarrow}$
(31)
$\displaystyle-\frac{{U}}{2}\left(\hat{n}_{0,\uparrow}+\hat{n}_{0,\downarrow}+\hat{n}_{1,\uparrow}+\hat{n}_{1,\downarrow}\right)$
$\displaystyle+h\left(\hat{n}_{0\uparrow}-\hat{n}_{0\downarrow}+\hat{n}_{1\uparrow}-\hat{n}_{1\downarrow}\right)$
$\displaystyle+\mu_{a}\left(\hat{n}_{0\uparrow}+\hat{n}_{0\downarrow}-\hat{n}_{1\uparrow}-\hat{n}_{1\downarrow}\right)$
$\displaystyle+h_{a}\left(\hat{n}_{0\uparrow}-\hat{n}_{0\downarrow}-\hat{n}_{1\uparrow}+\hat{n}_{1\downarrow}\right),$

and

$\hat{n}_{i\sigma}=\hat{c}_{i\sigma}^{\dagger}\hat{c}_{i\sigma},\quad
i=0,1,\quad\sigma=\uparrow,\downarrow.$ (32)

The parameters are chosen the same as in [24], namely
$t=1,\mu=0.7,h=0.3,U_{a}=0.5,\mu_{a}=0.2,h_{a}=0.03,\beta=31.1$ and $U=5$.

####  III.1.2 Fitting results, with comparison to other methods

Figure 1: Spectral functions of the Hubbard dimer system obtained from PES,
ES, Nevanlinna, Projection+Nevanlinna, Carathéodory, Projection+Carathéodory
and MaxEnt methods. Blue line: true value. Red line: calculation results from
analytic continuation.

Let us first demonstrate the efficiency of the PES method. We conduct analytic
continuation for four sets of Matsubara data with different noise levels:
$\sigma=0,\sigma=4\times 10^{-6},\sigma=2.56\times 10^{-4},\sigma=1.6384\times
10^{-2}$. In Fig. 1, we plot the spectral function (evaluated at
$\mathbb{R}+\mathrm{i}\eta$, $\eta=0.01$) obtained from analytic continuation
with comparison to its true values. The methods that we are comparing with PES
are Nevanlinna [23], Carathéodory [24] and MaxEnt [5]. We also implement the
ES method, in which we skip the projection step and conduct the estimation and
the SDR fitting directly. To make it a fair comparison, we also perform
Nevanlinna and Carathéodory continuations on the projected data. Since the
amplitude of MaxEnt spectral functions are very small compared to the delta
peaks in the true spectral functions, we rescale the MaxEnt results for
comparison with true values and results from other methods.

Our observation is summarized as follows:

  * •

When the Matsubara data is clean (noise-level $\sigma=0$), all methods except
the MaxEnt continuation could retrieve all the peaks perfectly. The MaxEnt
method can not deal with sharp features. In all cases, MaxEnt can only result
in two (and up to three) significantly broadened peaks.

  * •

Both the PES method and the ES method are robust to noise. It can retrieve
almost all features of the Green’s function when the noise level is small
($\sigma=4\times 10^{-6}$ and $2.56\times 10^{-4}$), and can retrieve quite a
few features when the noise level is relatively large ( $\sigma=1.6384\times
10^{-2}$). Particularly, the band gap could be accurately resolved even in the
large noise scenario.

  * •

For noisy data, PES behaves better than ES (see noise level $\sigma=4\times
10^{-6}$ and $2.56\times 10^{-4}$), which implies that the projection of noisy
data is indeed helpful. For noiseless data ($\sigma=0$), ES behaves better
than PES.

  * •

Even when the data is only slightly noisy (noise-level $\sigma=4\times
10^{-6}$), the Nevanlinna and Carathéodory continuation breaks down. This is
due to the violation of the Pick criterion [23]. The spectral function is no
longer nonnegative and multiple features are missing.

  * •

If combined with the projection of data into the causal space, the results of
the Nevanlinna and Carathéodory continuation are greatly improved. However,
the quality of Projection+Nevanlinna and Caratheodory is still not as good as
PES. This could be seen in all noisy levels. What’s more, the calculation
results are still not guaranteed to be causal, as shown in the results of
Projection+Carathéodory methods for relatively large noise (i.e.
$\sigma=2.56\times 10^{-4}$ and $\sigma=1.6384\times 10^{-2}$).

Since both PES method and Carathéodory continuation are able to obtain non-
diagonal entries of the Green’s functions, let us show the calculation result
of the real part of $G_{13}$ (also evaluated at $\mathbb{R}+\mathrm{i}\eta$,
$\eta=0.01$) in Fig. 2, in which we compare the performances of PES, ES,
Carathéodory and Projection+Carathéodory.

Figure 2: The non-diagonal element $G_{13}$ obtained by different analytic
continuation methods. Blue line: true value. Red line: calculation results
from analytic continuation.

Similar to what we have learned from the results of spectral functions, we
have the following observations:

  * •

The Projection+Carathéodory does not perform as good as the PES and ES method.
This could be seen at all noise levels. The two high peaks are perfectly
retrieved by the PES and ES method but not by Projection+Carathéodory.

  * •

At noise-level $\sigma=1.6384\times 10^{-2}$, the noise is too large and no
method could recover all peaks perfectly.

With these observations, we conclude that the PES method has the most robust
performance, especially for analytic continuation of the matrix-valued noisy
Matsubara data. The causal projection could help cure the noise-sensitivity
issue of Nevanlinna and Carathéodory continuation, but does not always
guarantee a causal result, and not necessarily perform as well as the PES
method.

####  III.1.3 Importance of the pole estimation

Now we would like to demonstrate the importance of the pole estimation step.
In [35], the SDR step is implemented in the context of hybridization fitting
without the pole estimation. A random initialization of poles could fail to
converge towards the global minimum of the loss function (see [35, Figure 2]).
In the context of analytic continuation, we also find that a random selection
of initial poles will not converge to its true positions. Particularly, the
high frequency features, i.e. the poles away from zero may not be identified
accurately.

To illustrate the difficulty in the SDR step, we fix the poles to be its true
values except for $\lambda_{1}$. The true value of $\lambda_{1}$ is $-4.6361$.
We choose its initial value to be $\lambda_{1}^{\text{in
}}=-5.25,-5.00,-4.75,-4.50,-4.25,-4.00$, and then conduct the SDR fitting
step. The corresponding result of the spectral function at the interval
$[-5,-4]$ is shown in Fig. 3 (top panel). We see that even such a small
deviation from the true value of $\lambda_{1}$ can result in large errors in
analytic continuation.

Figure 3: Top panel: calculation results of spectral functions (red dashed
lines) with comparison to the true value (blue lines) on the interval
$[-5,-4]$. The initial value of $\lambda_{1}$ are chosen as
$\lambda_{1}^{\text{in}}=-5.25,-5.00,-4.75,-4.50,-4.25$ and $-4.00$. Bottom
panel: landscape of $\mathcal{E}\left(\\{\lambda_{l}\\}_{l=1}^{N_{p}}\right)$
with respect to $\lambda_{1}$ on the interval $[-6,-4]$. Other poles are kept
fixed at their true positions: $-5.2479$, $-4.5143$, $-3.1957$, $-2.6075$,
$-2.4857$, $0.8619$, $0.9837$, $1.6181$, $3.4599$, $3.3381$ and $4.0255$.

The reason behind this is as follows. Let us keep other poles fixed at their
true values and let $\lambda_{1}$ vary from $-6$ to $-4$, and plot the
landscape of the loss function $\mathcal{E}(\\{\lambda_{l}\\})$ in Fig. 3
(bottom panel). We can see that the loss function is highly nonconvex. When
$\lambda^{\text{in}}_{1}$ is outside the interval $[-5.25,-4.5]$, the loss
function is essentially flat. In fact, if
$\lambda^{\text{in}}_{1}\notin[-5.25,-4.5]$, the optimized semidefinite matrix
$\mathbb{X}_{1}^{\mathrm{opt}}$ is a zero matrix. This will result in a
vanishing gradient for $\lambda_{1}^{\text{in}}$ (see Eq. 28), and therefore
the optimization could not continue. In such an adversarial scenario, a random
initialization of poles cannot recover the true pole locations.

###  III.2 Bosonic response functions

In order to demonstrate that the PES method is also applicable to bosonic
response functions, let us consider the quantity
$b(\tau)=\langle\hat{n}_{0\uparrow}(\tau)\hat{n}_{0\uparrow}(0)\rangle$ of the
Hubbard dimer model, which indicates the same-spin susceptibility for the
$0$-th site of the system. This function is also used in [2] to test the
performance of analytic continuation methods for bosonic functions. We use the
Hubbard dimer example with the same parameters as in Section III.1 other than
taking $U=0.9$.

Let us use $b(\omega)$ to denote the time-response function in the frequency
domain. We are trying to fit the Matsubara data of $b(\tau)$ into the
following form:

$b(w)=\sum_{l=1}^{N_{p}}\frac{x_{l}}{w-\lambda_{l}},$ (33)

where the quantities $x_{l}$ satisfies the semidefinite constraints Eq. 5 and
the constraint

$\displaystyle\sum_{l}{x}_{l}$
$\displaystyle=\frac{1}{Z}\sum_{r,s}\left|\left\langle\Psi_{s}\left|\hat{n}_{0\uparrow}\right|\Psi_{r}\right\rangle\right|^{2}\left(\mathrm{e}^{-\beta
E_{s}}-\mathrm{e}^{-\beta E_{r}}\right)=0.$ (34)

Here $N_{\text{orb}}=1$ since we are only considering the time response
function with respect to the $0$-th site.

In the current physical system, the spectral function of $b(w)$ has four
peaks, with two on the positive half axis, and two on the negative half axis.
In fact, $b(w)$ is an odd function. Since the susceptibility is mirrored for
negative frequencies, we only choose to plot the positive half axis. Note that
since the Nevanlinna and Carathéodory continuation are not applicable to
bosonic functions, we can only compare the PES method to the MaxEnt
continuation method, for Matsubara data with different noise levels
$\sigma=0,10^{-6},10^{-4}$ and $10^{-2}$. The result is shown in Fig. 4.

Figure 4: Analytic continuation of bosonic response functions.

The MaxEnt method could not resolve the peaks of the bosonic response function
under any noise level. We can see that when the noise level is small
($\sigma=0,10^{-6},10^{-4}$), the PES method could recover both peaks on the
positive half axis. When the noise level is large ($\sigma=10^{-2}$), the PES
method could only recover one peak. This means that the noise-robustness of
PES method for bosonic functions are comparable to that for fermionic
functions.

###  III.3 Band structure

We now applied the PES method to ab initio band structure calculations of
solids. We use two data sets obtained from GW calculations: the matrix-valued
Green’s functions for crystalline Si in the equilibrium geometry, previously
analyzed in Ref. [24], and for AgI obtained with a relativistic exact two-
component formulation in the one-electron approximation (x2c1e) [53]. The data
have been obtained by fully self-consistent GW formulated in Gaussian Bloch
orbitals, as described in detail in Ref. [54]. Both AgI and Si employ the gth-
dzvp-molopt-sr basis [55] with gth-pbe pseudopotential [56], using integrals
generated by the PySCF package [57] and a code based on ALPS [58].
Calculations are performed on a 6 $\times$ 6 $\times$ 6 grid [59] and
interpolated using a Wannier interpolation of Matsubara data [54] on 52
positive IR frequencies [60, 61] on a grid with 200 momentum points along a
high-symmetry path. The broadening parameter of $\eta=0.01\text{Ha}\approx
3157\text{K}$ is larger than the temperature of $T=0.001\text{Ha}\approx
316\text{K}$ and correlation effects.

Figure 5: Band structure of Silicon and AgI, revealed by the analytic
continuation. Data from Ref.[54]. White entries correspond to negative (non-
causal) spectral functions.

The analytic continuation is conducted as follows. Since the number of
orbitals is much larger than that in the Hubbard dimer example
($N_{\text{orb}}\approx 70$ for silicon and $N_{\text{orb}}\approx 230$ for
AgI), the computational cost will be too high if we conduct the analytic
continuation for the entire matrix. Therefore we only implement the analytic
continuation for the diagonal elements. The projection step is conducted on
each diagonal entry separately with a fine grid of $M=12000$ points. After the
projection step, we estimate the poles of each diagonal entry using AAA
algorithm, and then conduct the SDR fitting of each diagonal entry separately.

We conduct analytic continuation for Matsubara data with different noise
levels $\sigma=0,10^{-6},10^{-4}$ and $10^{-2}$. The results of the spectral
functions along a high symmetry path reveal the band structure, and are
plotted in Fig. 5. We place a white marker in the plot if the corresponding
component of the spectral function becomes negative to indicate that the
result is unphysical.

Once again, we can see that Nevanlinna continuations could result in negative
spectral functions even for the clean data, while Nevanlinna continuations
combined with the causal projection can largely resolve the isuse of noise-
sensitivity. PES method, ES method and the Nevanlinna continuation on
projected data routine are all noise-robust. Even at a relatively large noise-
level $\sigma=10^{-2}$, the structure of the valence and conduction bands can
be recovered. This means that our noise-robust analytic continuation methods
could successfully retrieve information such as the band gap from Matsubara
data even under a significant noise level.

##  IV Discussion

In this article, we have proposed a robust analytic continuation method called
PES, consisting of the following three steps: causal projection, pole
estimation and SDR fitting through bi-level optimization. The PES method is
noise-robust, could deal with matrix-valued Green’s functions and is
applicable to both fermionic and bosonic systems. As a by product, we also
find that the causal projection could significantly improve the performance of
noise-sensitive methods such as Nevanlinna and Carathéodory continuation. We
have also demonstrated these properties with extensive numerical experiments.

Our current strategy of pole estimation is to perform AAA algorithm on the
projected data. Note that the AAA algorithm is still an interpolation
procedure can be susceptible to noise. The causal projection alleviates this
drawback, and in practice it often provides a good enough initial guess for
SDR fitting even in the noisy scenario. A better strategy to estimate the pole
locations (e.g., by combining with the Prony method [33]) may still be helpful
in improving the overall performance of the PES method. Moreover, the current
SDP solver (SDPT3) can be computationally expensive when the number of
variables is large. In such a case, more efficient solvers based on first-
order methods may reduce the computational time in both the projection and the
SDR fitting steps.

As we have shown, the PES method is useful for retrieving the excitation
spectra of molecules and band structures in solids. These data feature
multiple sharp peaks modeled by poles on the real axis. It could also be used
to perform bath fitting in the Dynamical Mean Field Theory (DMFT) [35], and it
is applicable in lattice gauge theory (for example, see [62, 63]). Broadened
spectral functions that cannot be resolved into discrete peaks may be modeled
by a relatively small number of poles away from the real axis. This would
require modification of the definition of the causal space. Note that in the
noisy scenario, it may be very difficult to distinguish the spectra described
by a pole in the complex plane and the spectra due to many poles on or near
the real axis. Therefore prior information such as the maximum number of
complex poles and approximate pole locations may become necessary. This
scenario will be investigated in future work.

Acknowledgment

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Advanced Scientific Computing Research and Office
of Basic Energy Sciences, Scientific Discovery through Advanced Computing
(SciDAC) program under Award Number DE-SC0022198 (Z.H.) and DE-SC0022088
(E.G.). The work is partially supported by the Air Force Office of Scientific
Research under award number FA9550-18-1-0095 (L.L.). L.L. is a Simons
Investigator. We thank Jiahao Yao, Lexing Ying, Yang Yu and Dominika Zgid for
helpful discussions.

## References

  * Baker et al. [1996] G. A. Baker, G. A. Baker Jr, P. Graves-Morris, G. Baker, and S. S. Baker, _Pade Approximants: Encyclopedia of Mathematics and It’s Applications, Vol. 59 George A. Baker, Jr., Peter Graves-Morris_ , vol. 59 (Cambridge University Press, 1996).
  * Schött et al. [2016] J. Schött, I. L. Locht, E. Lundin, O. Grånäs, O. Eriksson, and I. Di Marco, Physical Review B 93, 075104 (2016).
  * Jarrell and Gubernatis [1996] M. Jarrell and J. E. Gubernatis, Physics Reports 269, 133 (1996).
  * Asakawa et al. [2001] M. Asakawa, Y. Nakahara, and T. Hatsuda, Progress in Particle and Nuclear Physics 46, 459 (2001).
  * Levy et al. [2017] R. Levy, J. LeBlanc, and E. Gull, Computer Physics Communications 215, 149 (2017).
  * Bryan [1990] R. Bryan, European Biophysics Journal 18, 165 (1990).
  * Rothkopf [2020] A. Rothkopf, Data 5, 85 (2020).
  * Sandvik [1998] A. W. Sandvik, Physical Review B 57, 10287 (1998).
  * Vafayi and Gunnarsson [2007] K. Vafayi and O. Gunnarsson, Physical Review B 76, 035115 (2007).
  * Fuchs et al. [2010] S. Fuchs, T. Pruschke, and M. Jarrell, Physical Review E 81, 056701 (2010).
  * Shao and Sandvik [2022] H. Shao and A. W. Sandvik, _Progress on stochastic analytic continuation of quantum monte carlo data_ (2022), eprint 2202.09870.
  * Yoshimi et al. [2019] K. Yoshimi, J. Otsuki, Y. Motoyama, M. Ohzeki, and H. Shinaoka, Computer Physics Communications 244, 319 (2019).
  * Motoyama et al. [2022] Y. Motoyama, K. Yoshimi, and J. Otsuki, Physical Review B 105, 035139 (2022).
  * Fournier et al. [2020] R. Fournier, L. Wang, O. V. Yazyev, and Q. Wu, Physical Review Letters 124, 056401 (2020).
  * Yoon et al. [2018] H. Yoon, J.-H. Sim, and M. J. Han, Physical Review B 98, 245101 (2018).
  * Tomczak and Biermann [2007] J. M. Tomczak and S. Biermann, Journal of Physics: Condensed Matter 19, 365206 (2007).
  * Dang et al. [2014] H. T. Dang, X. Ai, A. J. Millis, and C. A. Marianetti, Physical Review B 90, 125114 (2014).
  * Gull and Millis [2014] E. Gull and A. J. Millis, Physical Review B 90, 041110 (2014).
  * Kraberger et al. [2017] G. J. Kraberger, R. Triebl, M. Zingl, and M. Aichhorn, Physical Review B 96, 155128 (2017).
  * Reymbaut et al. [2017] A. Reymbaut, A.-M. Gagnon, D. Bergeron, and A.-M. S. Tremblay, Phys. Rev. B 95, 121104 (2017).
  * Kac and Krein [1974] I. S. Kac and M. G. Krein, Amer. Math. Soc. Transl 103, 1–18 (1974).
  * Belyi et al. [2006] S. Belyi, S. Hassi, H. de Snoo, and E. Tsekanovskii, Linear algebra and its applications 419, 331 (2006).
  * Fei et al. [2021a] J. Fei, C.-N. Yeh, and E. Gull, Physical Review Letters 126, 056402 (2021a).
  * Fei et al. [2021b] J. Fei, C.-N. Yeh, D. Zgid, and E. Gull, Physical Review B 104, 165111 (2021b).
  * Pick [1917] G. Pick, Mathematische annalen 78, 270 (1917).
  * Khargonekar and Tannenbaum [1985] P. Khargonekar and A. Tannenbaum, IEEE Transactions on Automatic Control 30, 1005 (1985).
  * Tannenbaum [1987] A. Tannenbaum, in _26th IEEE Conference on Decision and Control_ (IEEE, 1987), vol. 26, pp. 1635–1636.
  * Doyle et al. [2013] J. C. Doyle, B. A. Francis, and A. R. Tannenbaum, _Feedback control theory_ (Courier Corporation, 2013).
  * Osborne and Smyth [1995] M. R. Osborne and G. K. Smyth, SIAM J. Sci. Comput. 16, 119 (1995).
  * Golub and Pereyra [2003] G. Golub and V. Pereyra, Inverse problems 19, R1 (2003).
  * Beylkin and Monzón [2005] G. Beylkin and L. Monzón, Applied and Computational Harmonic Analysis 19, 17 (2005).
  * Ying [2022a] L. Ying, arXiv preprint arXiv:2202.09719 (2022a).
  * Ying [2022b] L. Ying, arXiv preprint arXiv:2202.02670 (2022b).
  * Nakatsukasa et al. [2018] Y. Nakatsukasa, O. Sète, and L. N. Trefethen, SIAM Journal on Scientific Computing 40, A1494 (2018).
  * Mejuto-Zaera et al. [2020] C. Mejuto-Zaera, L. Zepeda-Núñez, M. Lindsey, N. Tubman, B. Whaley, and L. Lin, Physical Review B 101, 035143 (2020).
  * Negele and Orland [1988] J. W. Negele and H. Orland, _Quantum many-particle systems_ (Westview, 1988).
  * Gesztesy and Tsekanovskii [2000] F. Gesztesy and E. Tsekanovskii, Mathematische Nachrichten 218, 61 (2000).
  * Shohat and Tamarkin [1950] J. A. Shohat and J. D. Tamarkin, _The problem of moments_ , vol. 1 (American Mathematical Society (RI), 1950).
  * Akhiezer [2020] N. I. Akhiezer, _The classical moment problem and some related questions in analysis_ (SIAM, 2020).
  * Grant and Boyd [2014] M. Grant and S. Boyd, _CVX: Matlab software for disciplined convex programming, version 2.1_ (2014).
  * Tibshirani [1996] R. Tibshirani, J. Royal Stat. Soc. B 58, 267 (1996).
  * Candes et al. [2006] E. J. Candes, J. K. Romberg, and T. Tao, Communications on Pure and Applied Mathematics 59, 1207 (2006).
  * Bach, Francis and Jenatton, Rodolph and Mairal, Julien and Obozinski, Guillaume [2012] Bach, Francis and Jenatton, Rodolph and Mairal, Julien and Obozinski, Guillaume (2012).
  * Driscoll et al. [2014] T. A. Driscoll, N. Hale, and L. N. Trefethen, _Chebfun guide_ (2014).
  * Berrut and Trefethen [2004] J.-P. Berrut and L. N. Trefethen, SIAM Rev. 46, 501 (2004).
  * Feynman [1939] R. P. Feynman, Physical review 56, 340 (1939).
  * Nocedal and Wright [1999] J. Nocedal and S. J. Wright, _Numerical optimization_ (Springer Verlag, 1999).
  * Toh et al. [2012] K.-C. Toh, M. J. Todd, and R. H. Tütüncü, in _Handbook on semidefinite, conic and polynomial optimization_ (Springer, 2012), pp. 715–754.
  * Mehrotra [1992] S. Mehrotra, SIAM Journal on optimization 2, 575 (1992).
  * Vandenberghe and Boyd [1996] L. Vandenberghe and S. Boyd, SIAM Rev. 38, 49 (1996).
  * Jiang, Kaifeng and Sun, Defeng and Toh, Kim-Chuan [2012] Jiang, Kaifeng and Sun, Defeng and Toh, Kim-Chuan, SIAM J. Optim. 22, 1042 (2012).
  * Nakatsukasa and Trefethen [2020] Y. Nakatsukasa and L. N. Trefethen, SIAM Journal on Scientific Computing 42, A3157 (2020).
  * Yeh et al. [2022a] C.-N. Yeh, A. Shee, Q. Sun, E. Gull, and D. Zgid (2022a).
  * Yeh et al. [2022b] C.-N. Yeh, S. Iskakov, D. Zgid, and E. Gull, _Fully self-consistent finite-temperature $gw$ in gaussian bloch orbitals for solids_ (2022b).
  * VandeVondele and Hutter [2007] J. VandeVondele and J. Hutter, The Journal of Chemical Physics 127, 114105 (2007), eprint https://doi.org/10.1063/1.2770708.
  * Goedecker et al. [1996] S. Goedecker, M. Teter, and J. Hutter, Phys. Rev. B 54, 1703 (1996).
  * Sun et al. [2020] Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, et al., The Journal of Chemical Physics 153, 024109 (2020), eprint https://doi.org/10.1063/5.0006074.
  * Gaenko et al. [2017] A. Gaenko, A. Antipov, G. Carcassi, T. Chen, X. Chen, Q. Dong, L. Gamper, J. Gukelberger, R. Igarashi, S. Iskakov, et al., Computer Physics Communications 213, 235 (2017), ISSN 0010-4655.
  * Iskakov et al. [2020] S. Iskakov, C.-N. Yeh, E. Gull, and D. Zgid, Phys. Rev. B 102, 085105 (2020).
  * Shinaoka et al. [2017] H. Shinaoka, J. Otsuki, M. Ohzeki, and K. Yoshimi, Phys. Rev. B 96, 035147 (2017).
  * Li et al. [2020] J. Li, M. Wallerberger, N. Chikano, C.-N. Yeh, E. Gull, and H. Shinaoka, Phys. Rev. B 101, 035144 (2020).
  * Burnier and Rothkopf [2013] Y. Burnier and A. Rothkopf, Phys. Rev. Lett. 111, 182003 (2013).
  * Tripolt et al. [2019] R.-A. Tripolt, P. Gubler, M. Ulybyshev, and L. von Smekal, Computer Physics Communications 237, 129 (2019), ISSN 0010-4655.