# Analytical Continuation of Matrix-Valued Functions: Carathéodory Formalism

Jiani Fei  Department of Physics, University of Michigan, Ann Arbor, Michigan
48109, USA  Chia-Nan Yeh  Department of Physics, University of Michigan, Ann
Arbor, Michigan 48109, USA  Dominika Zgid  Department of Chemistry, University
of Michigan, Ann Arbor, Michigan 48109, USA  Department of Physics, University
of Michigan, Ann Arbor, Michigan 48109, USA  Emanuel Gull  Department of
Physics, University of Michigan, Ann Arbor, Michigan 48109, USA

###### Abstract

Finite-temperature quantum field theories are formulated in terms of Green’s
functions and self-energies on the Matsubara axis. In multi-orbital systems,
these quantities are related to positive semidefinite matrix-valued functions
of the Carathéodory and Schur class. Analysis, interpretation and evaluation
of derived quantities such as real-frequency response functions requires
analytic continuation of the off-diagonal elements to the real axis. We derive
the criteria under which such functions exist for given Matsubara data and
present an interpolation algorithm that intrinsically respects their
mathematical properties. For small systems with precise Matsubara data, we
find that the continuation exactly recovers all off-diagonal and diagonal
elements. In real-materials systems, we show that the precision of the
continuation is sufficient for the analytic continuation to commute with the
Dyson equation, and we show that the commonly used truncation of off-diagonal
self-energy elements leads to considerable approximation artifacts. Our method
paves the way for the systematic evaluation of Matsubara data with equations
of many-body theory on the real-frequency axis.

##  I Introduction

The central object of finite-temperature field theory is the Matsubara Green’s
function $\mathcal{G}_{ij}(i\omega_{n})$. Finite temperature simulations
ranging from perturbative calculations [1, 2, 3] to lattice [4] and
continuous-time [5] quantum Monte Carlo and lattice QCD [6, 7, 8] simulations
obtain this quantity. In post-processing, an analytical continuation step to
the retarded real axis Green’s function $G^{R}_{ij}(\omega)$ is performed with
the purpose of obtaining single-particle excitation spectra, self-energy
information, and approximations to susceptibilities that can then directly be
related to experiment.

As small imprecisions in the Matsubara data result in large deviations of the
continued quantities, the direct solution of the continuation problem is
typically avoided. Instead, methods that aim to fit the Matsubara data with a
physically reasonable (i.e. smooth, positive and normalized) spectral
function, such as the maximum entropy analytic continuation [9, 10, 11, 12,
13, 14, 15, 16, 17, 18], the stochastic analytic continuation (SAC) and
variants [19, 20, 21, 22, 23, 24], the sparse modeling method [25, 26], or
machine learning approaches [27] are used, from which the real and imaginary
parts of the retarded Green’s functions, self-energies, and susceptibilities
are extracted. These methods work well for single-orbital systems and the
diagonal components of Green’s functions, especially in the presence of data
with statistical uncertainties.

However, Green’s functions for all but the simplest quantum systems are
matrix-valued objects with both diagonal and off-diagonal entries. While
certain quantities such as the total spectral function only depend on the
diagonal entries, the analysis of the self-energy, computation of the
susceptibilities, or evaluation of equations such as the Dyson equation on the
real axis requires knowledge of both diagonal and off-diagonal entries. As the
off-diagonal entries may change sign as a function of frequency, standard
continuation methods that rely on positivity [11] fail.

So far, no reliable and general algorithms for the continuation of matrix-
valued Green’s functions exist, even though several approaches have been
explored. For instance, one may perform a Padé [28, 29, 13, 30, 31, 32, 33]
continuous fraction interpolation of both the diagonal and off-diagonal terms.
However, the continued results typically exhibit continuation artifacts such
as negative spectral functions even for the diagonal part. One may instead
transform to a basis that diagonalizes the Green’s function or the self-energy
for a given Matsubara frequency [34, 35, 36], neglect the remaining off-
diagonal elements, and employ a continuation method for diagonal Green’s
functions. Unless symmetry dictates that all off-diagonal elements for all
frequencies must be zero, this is an uncontrolled approximation and results
will depend on the basis chosen. Finally, one may generalize the maximum
entropy method to matrix-valued functions and off-diagonal terms [37, 38],
using either a positive-negative or a maximum ‘quantum’ entropy approach.
These methods enforce the positive semidefiniteness of the Green’s function
but eliminate sharp and high-energy features, such as the band structure
contained in Green’s function data, due to the intrinsic limitations of the
fitting procedure.

As we show in this paper, the application of Nevalinna theory [39] to matrix-
valued functions overcomes these limitations, leading to basis-indepedent
interpolations that intrinsically respect the analytic structure of matrix-
valued Green’s functions. Using a generalization of Nevanlinna-Pick
interpolation to the class of Carathéodory functions, we demonstrate that the
analytic continuation of Green’s functions, self-energies, and cumulants [40]
leads to matrix-valued functions that are indistinguishable from the exact
results for a simple model system. In the context of a real-materials multi-
orbital simulation we then demonstrate that the analytic continuation step
commutes with the Dyson equation, such that self-energy, cumulant, and Green’s
function continuations yield consistent results. Finally, we assess the
quality of the frequently used approximation of neglecting off-diagonal self-
energy components.

##  II Theory

The main mathematical objects considered in this paper are matrix-valued
Green’s functions, self-energies, and cumulants. As we will show below, the
analytical properties of these functions ensure that they lie, up to a factor
of $i=\sqrt{-1}$, in the class of Carathéodory functions [41].

In the following, we will present the mathematical framework by introducing
Carathéodory functions and their properties, along with the related class of
Schur functions [42] and a mapping between the two classes. We will then
reformulate the problem of analytical continuation as an interpolation problem
in the class of Schur or Carathéodory functions, in analogy to Ref. [39].
After detailing the conditions under which a Carathéodory interpolant through
given Matsubara points exists, we will derive an interpolation algorithm.
Finally, we will show that Green’s functions, self-energies, and cumulants are
Carathéodory functions.

###  II.1 Carathéodory and Schur functions

Consider an open subset $\mathcal{B}$ of the complex plane, such as the unit
disk $\mathcal{D}=\\{z\mathrel{\mathop{\mathchar 58\relax}}\lvert
z\rvert<1\\}$ or the upper half of the complex plane
$\mathcal{C}^{+}=\\{z\mathrel{\mathop{\mathchar 58\relax}}\Im{z}>0\\}$. A
matrix-valued function $F(z)\mathrel{\mathop{\mathchar 58\relax}}\
\mathcal{B}\rightarrow\mathbb{C}^{m\times m}$, holomorphic on $\mathcal{B}$,
belongs to the class of Carathéodory functions $\mathfrak{C}$ if, for any
$z\in\mathcal{B}$, the Hermitian matrix $(F(z)+F^{\dagger}(z))/2$ is positive
semi-definite (PSD) [41, 43, 44].

Analogously, a matrix-valued function $\Psi(z),\
\mathcal{D}\rightarrow\mathbb{C}^{m\times m}$, holomorphic on $\mathcal{D}$,
belongs to the Schur class $\mathcal{S}$ if $\lVert\Psi(z)\rVert\leq 1$ for
any $\lvert z\rvert<1$, where $\lVert\Psi\rVert$ denotes norm of the matrix
$\Psi$ [42, 45]. We use the spectral norm as the matrix norm in this paper,
i.e. the largest eigenvalue of the constant matrix
$\Psi_{s}\mathrel{\mathop{\mathchar 58\relax}}=[\Psi\Psi^{\dagger}]^{1/2}$,
with $(\cdot)^{1/2}$ indicating the Hermitian square root.

While the domain of interest for many-body objects such as the Green’s
function is the upper half of the complex plane $\mathcal{C}^{+}$, the
traditional mathematical literature mostly considers functions on the open
unit disk $\mathcal{D}$. A Möbius transform $h\mathrel{\mathop{\mathchar
58\relax}}\mathcal{C}^{+}\rightarrow\mathcal{D},h(z)=\frac{z-i}{z+i}$ maps
$\mathcal{C}^{+}$ to $\mathcal{D}$, and its inverse
$h^{-1}(z)\mathrel{\mathop{\mathchar
58\relax}}\mathcal{D}\rightarrow\mathcal{C}^{+}$ maps the unit disk back to
the upper half of the complex plane.

As shown in detail in [43, 45], every Carathéodory function
$F(z)\in\mathfrak{C}$ restricted to $\mathcal{D}$ can be mapped to a
corresponding Schur function $\Psi(z)\in\mathcal{S}$ with the Cayley transform

$\displaystyle\Psi(z)=[I-F(z)][I+F(z)]^{-1},$ (1)

and its inverse

$\displaystyle F(z)=[I+\Psi(z)]^{-1}[I-\Psi(z)].$ (2)

Figure 1: Mapping of the input Carathéodory problem to a Schur interpolation
problem. $x_{j}$ denote Matsubara frequencies, $Y_{j}$ Matsubara values,
$z_{j}$ transformed Matsubara points, and $J_{j}$ transformed Matsubara values

As illustrated in Fig. 1, the construction of a Carathéodory interpolant on
the upper half of the complex plane is therefore equivalent to the
construction of a Carathéodory interpolant on the unit disk, which in turn is
equivalent to the construction of the corresponding Schur interpolant on the
unit disk.

###  II.2 Pick Criterion

We aim to interpolate matrix-valued Green’s functions, self-energies, and
cumulants obtained in Matsubara frequencies. Thus, after mapping to the unit
disk, the interpolation problem is specified by a set of $n$ Carathéodory
$m\times m$ matrices $Y_{j}(z_{j})$ at $n$ points $z_{j}\in\mathcal{D}$ (see
Fig. 1). We first consider the conditions under which such an interpolation
problem has solutions.

Generalizing the Pick criterion for scalar functions [46, 47, 48] to the
matrix-valued case, Refs. [43, 49] derived an existence criterion for
Carathéodory interpolants directly based on input data. Solutions for the
interpolation problem exist if and only if the Pick matrix defined for the
Carathéodory function on the unit disk $\mathcal{D}$,

$\displaystyle
P_{C}=[\frac{Y_{k}+Y_{l}^{*}}{1-z_{k}^{*}z_{l}}]_{(mn)\times(mn)}$ (3)

or alternatively, the Pick matrix defined for the transformed Schur function
on $\mathcal{D}$,

$\displaystyle
P_{S}=[\frac{I-J_{k}^{*}J_{l}}{1-z_{k}^{*}z_{l}}]_{(mn)\times(mn)}$ (4)

is positive semi-definite; and a unique solution only exists if it is
singular.

We note that this criterion is very restrictive in practice. Numerical noise
in Monte Carlo typically leads to negative eigenvalues such that a
Carathéodory (or Nevanlinna [39]) solution of the interpolation problem does
not exist. Round-off and convergence issues in semi-analytical calculations
such as GW simulations of real materials (see Sec. III) lead to eigenvalues
that are small and negative, such that the evaluation just above the real axis
still leads to PSD spectral functions, while analytical solutions such as
those for the Hubbard dimer presented below always satisfy the Pick criterion.

###  II.3 Interpolation of Schur functions

Matrix-valued functions in the Carathéodory class have continued fraction
expansions [43]. An algorithm to obtain such interpolation is given by matrix
extensions of the classical Schur algorithm [43, 49, 50], which we use as
interpolation method in this paper.

Given input data $F(x_{j})=Y_{j}$ ($j=0,1,\dots,n-1$;
$x_{j}\in\mathcal{C}^{+}$; $Y_{j}\in\mathbb{C}^{m\times m}$) for a
Carathéodory matrix-valued function $F(z)$, we Möbius transform the domain and
conformally map the function value according to Eq. 1 (see Fig. 1), and
reformulate the problem as a problem of finding a matrix-valued function
$\Psi(z)$ in the Schur class, such that

$\displaystyle\Psi(z_{j})$
$\displaystyle=J_{j}=\Psi(\frac{x_{j}-i}{x_{j}+i})=[I-Y_{j}][I+Y_{j}]^{-1},$
$\displaystyle j$ $\displaystyle=0,1,\dots,n-1.$ (5)

We will proceed as follows. Given a interplation problem $\Psi=\Psi_{0}$ with
$n$ nodes $\Psi(z_{j})=J_{j}$, $j=0,1,\dots,n-1$, we will find a function that
interpolates the first node $z_{0}$ and express the remaining problem as an
interpolation problem $\Psi_{1}$ through $n-1$ nodes. This interpolation
problem will then be expressed as a function that interpolates the second node
$z_{1}$ and a remaining problem $\Psi_{2}$ through $n-2$ nodes. The procedure
will be repeated until all nodes are interpolated and only a free Schur
function $\Psi_{n}$ remains.

We first show how to reduce the Schur function $\Psi_{i}(z)\in\mathcal{S}$ to
$\Psi_{i+1}(z)\in\mathcal{S}$, while releasing the node constraint of
$\Psi_{i}(z)$ at $z_{i}$. The reduction step is based on the theories of
J-contractive transformations [51] and as follows [43]. Assuming
$\Psi_{i}(z_{i})=W_{i}$, define the matrix-valued function $L_{i}(z)$ by

$\displaystyle y_{i}L_{i}(z)=$
$\displaystyle[I-W_{i}W_{i}^{\dagger}]^{-1/2}[\Psi_{i}(z)-W_{i}]\cdot$
$\displaystyle[I-W_{i}^{\dagger}\Psi_{i}(z)]^{-1}[I-W_{i}^{\dagger}W_{i}]^{1/2}$
(6)

where $y_{i}=\lvert z_{i}\rvert(z_{i}-z)/(z_{i}(1-z_{i}^{*}z))$.
$L_{i}(z)\in\mathcal{S}$ by the Schwarz lemma [51]. Define $\Psi_{i+1}(z)$ as

$\displaystyle\Psi_{i+1}(z)=$
$\displaystyle[I-K_{i}K_{i}^{\dagger}]^{-1/2}[L_{i}(z)-K_{i}]\cdot$
$\displaystyle[I-K_{i}^{\dagger}L_{i}(z)]^{-1}[I-K_{i}^{\dagger}K_{i}]^{1/2}$
(7)

where $K_{i}$ is an arbitrary matrix such that $\lVert K_{i}\rVert<1$. Since
for any $z_{k},z_{l}\in\mathcal{D}$, the matrices
$\left[I-L(z_{k})^{\dagger}L(z_{l})\right]$ and
$\left[I-\Psi_{i+1}(z_{k})^{\dagger}\Psi_{i+1}(z_{l})\right]$ are equivalent
under similarity transformation, $\Psi_{i+1}(z)\in\mathcal{S}$ as well.

Notice that in Eq. 6, left- and right-hand-side are both zero at node $z_{i}$
so that the value of $L_{i}(z_{i})$, and in turn of $\Psi_{i+1}(z_{i})$, is
free; and that for all remaining nodes $z_{j},j>i$, $\Psi_{i+1}(z_{j})$ is
completely determined by $L_{i}(z_{j})$, in turn by $\Psi_{i}(z_{j})$, up to a
freedom in the choice of $K_{i}$. We therefore have a new interpolation
problem $\Psi_{i+1}\in\mathcal{S}$ with one less node constraint. Iterating
the complete algorithm backwards, an arbitrary Schur function
$\Psi_{n}(z)\in\mathcal{S}$ will yield $\Psi_{0}(z)\in\mathcal{S}$ that hits
all interpolation nodes.

$K_{i}$, $i=0,1,\dots,n-1$ and $\Psi_{n}$ are free parameters that can be used
to enforce additional conditions, such as smoothness [39], and cover all
possible interpolants in the class of Schur functions. For convenience, we
choose $\Psi_{n}$ to be the identity matrix and $K_{i}$ to be zero, for all
$i$.

Denoting $\Psi_{i}(z_{j})=W_{j}^{i}$, $i,j=0,1,\dots,{n-1}$ and
$W_{i}^{i}=W_{i}$, consistent with the above notation, the first stage of the
algorithm consists of computing all $W_{i}$ and storing them. By Eq. (6) and
(7), we have

$\displaystyle W_{j}^{i+1}=$
$\displaystyle\frac{z_{i}(1-z_{i}^{*}z_{j})}{\lvert
z_{i}\rvert(z_{i}-z_{j})}[1-W_{i}W_{i}^{\dagger}]^{-1/2}[W_{j}^{i}-W_{i}]\times$
$\displaystyle[1-W_{i}^{\dagger}W_{j}^{i}]^{-1}[1-W_{i}^{\dagger}W_{i}]^{1/2}\quad
j\geq i+1$ (8)

Iterating through $i=0,\dots,n-2$, we obtain all $W_{i}$.

The second stage of the algorithm consists of solving
$\Psi_{n}(z)\rightarrow\Psi_{n-1}(z)\dots\rightarrow\Psi_{0}(z)=\Psi(z)$ where
$\Psi(z),z\in\mathcal{D}$ is the desired Schur interpolant. This is done by
transcribing Eq. (6) and (7) as

$\displaystyle V_{i}(z)=$ $\displaystyle\frac{\lvert
z_{i}\rvert(z_{i}-z)}{z_{i}(1-z_{i}^{*}z)}$ (9)
$\displaystyle\times[I-W_{i}W_{i}^{\dagger}]^{1/2}\Psi_{i+1}(z)[I-W_{i}^{\dagger}W_{i}]^{-1/2}$
$\displaystyle\Psi_{i}(z)=$
$\displaystyle[I+V_{i}(z)W_{i}^{\dagger}]^{-1}[V_{i}(z)+W_{i}]$ (10)

with $i$ iterating from $n-1$ to $0$.

Lastly, we map the Schur class fuctions back to the Carathéodory space by
constructing the Carathéodory interpolant $F(x)$ via

$\displaystyle F(x)=[I+\Psi(\frac{x-i}{x+i})]^{-1}[I-\Psi(\frac{x-i}{x+i})]$
(11)

The first stage is computed only once, while the second stage depends on $z$
and is repeated by Eq. 11 for each $x\in\mathcal{C}^{+}$ that needs to be
evaluated.

We note that, as in most Padé codes, our numerical interpolation is highly
sensitive to numerical roundoff and is therefore performed in $512$ or
$1024$-bit precision.

Figure 2: Schematic of the transformations between $G(i\omega_{n})$,
$G(\omega)$, $\Sigma(i\omega_{n})$, $\Sigma(\omega)$, $M(i\omega_{n})$ and
$M(\omega)$. AC stands for analytical continuation. BC stands for back
continuation.

###  II.4 Carathéodory Functions in Many-Body Theory

We now show that the Green’s functions, self-energies, and cumulants commonly
encountered in many-body theory are Carathéodory functions (up to factors of
$i$) and can therefore be interpolated with the interpolation method described
above.

####  II.4.1 $iG(z)$ is a Carathéodory function

We assume a Hamiltonian system with eigenvectors $\lvert m\rangle$ and eigen
energies $\textit{E}_{m}$. $c_{\bm{i}}^{\dagger}$ ($c_{\bm{i}}$) is the
creation (annihilation) operator for the single-particle orbital $\bm{i}$,
$Z=\sum_{m}e^{-\beta\textit{E}_{m}}$ the partition function and $\beta$ the
inverse temperature. In the Lehmann representation, the matrix elements of the
fermionic Green’s function $G(z)$ in the upper half complex plane
$\mathcal{C}^{+}$ are

$\displaystyle G_{\bm{i}\bm{j}}(z)$
$\displaystyle=\frac{1}{Z}\sum_{m,n}\frac{\langle n\lvert c_{\bm{i}}\rvert
m\rangle\langle m\lvert c^{\dagger}_{\bm{j}}\rvert
n\rangle}{z+\textit{E}_{n}-\textit{E}_{m}}(e^{-\beta\textit{E}_{n}}+e^{-\beta\textit{E}_{m}}).$
(12)

We aim to prove for $z\in\mathcal{C}^{+}$ that $iG(z)+(iG(z))^{\dagger}$ is a
PSD matrix, i.e., $iG(z)$ is a Carathéodory function on $\mathcal{C}^{+}$.
This follows directly from the definition of a PSD matrix. For any complex
vector $\rvert x\rangle$,

$\displaystyle\langle x\lvert iG(z)$ $\displaystyle+(iG(z))^{\dagger}\rvert
x\rangle$ (13) $\displaystyle=\frac{1}{Z}\sum_{mn\bm{i}\bm{j}}$
$\displaystyle\left[\frac{i(e^{-\beta\textit{E}_{m}}+e^{-\beta\textit{E}_{n}})}{z+\textit{E}_{n}-\textit{E}_{m}}-\frac{i(e^{-\beta\textit{E}_{m}}+e^{-\beta\textit{E}_{n}})}{z^{*}+\textit{E}_{n}-\textit{E}_{m}}\right]$
$\displaystyle\langle n\lvert c_{\bm{i}}x_{\bm{i}}^{*}\rvert m\rangle\langle
m\lvert c^{\dagger}_{\bm{j}}x_{\bm{j}}\rvert n\rangle$ (14)
$\displaystyle=\frac{1}{Z}\sum_{mn}$
$\displaystyle\left[\frac{i(e^{-\beta\textit{E}_{m}}+e^{-\beta\textit{E}_{n}})}{z+\textit{E}_{n}-\textit{E}_{m}}-\frac{i(e^{-\beta\textit{E}_{m}}+e^{-\beta\textit{E}_{n}})}{z^{*}+\textit{E}_{n}-\textit{E}_{m}}\right]$
$\displaystyle\langle n\lvert\sum_{\bm{i}}c_{\bm{i}}x_{\bm{i}}^{*}\rvert
m\rangle^{2}$ (15) $\displaystyle=\frac{1}{Z}\sum_{mn}$
$\displaystyle\frac{2\text{Im}\\{z\\}\
(e^{-\beta\textit{E}_{m}}+e^{-\beta\textit{E}_{n}})}{\text{Im}\\{z\\}^{2}+(\text{Re}\\{z\\}+\textit{E}_{n}-\textit{E}_{m})^{2}}$
$\displaystyle\langle n\lvert\sum_{\bm{i}}c_{\bm{i}}x_{\bm{i}}^{*}\rvert
m\rangle^{2}\geq 0$ (16)

$iG(z)$ is therefore a Carathéodory function on $\mathcal{C}^{+}$.

####  II.4.2 $i\Sigma(z)$ is a Carathéodory function

In order to show that self-energies are Carathéodory functions, we make use of
the Lehmann representation of the self-energy $\Sigma(z)$ proposed in Ref.
[52, 53].

The Lehmann representation of $\Sigma(z)$ is constructed with respect to a
general fermionic Hamiltonian

$\displaystyle H(t)=$
$\displaystyle\sum_{\bm{i}\bm{j}}[T_{\bm{i}\bm{j}}(t)-\mu\delta_{\bm{i}\bm{j}}]c_{\bm{i}}^{\dagger}(t)c_{\bm{j}}(t)$
$\displaystyle+\frac{1}{2}\sum_{\bm{i}\bm{j}\bm{i}^{\prime}\bm{j}^{\prime}}U_{\bm{i}\bm{i}^{\prime}\bm{j}\bm{j}^{\prime}}(t)c_{\bm{j}}(t)c_{\bm{i}^{\prime}}^{\dagger}(t)c_{\bm{j}^{\prime}}(t)c_{\bm{i}}(t),$
(17)

where $c_{\bm{i}}(t)=U^{\dagger}(t,0)c_{\bm{i}}U(t,0)$,
$U(t,t^{\prime})=\mathcal{T}_{C}\exp(-i\int_{t}^{t^{\prime}}H(t_{1})dt_{1})$
is the system’s time-evolution operator and $\mathcal{T}_{C}$ is the time-
ordering operator along the Keldysh-Matsubara contour
($t=0\rightarrow\infty\rightarrow 0\rightarrow-i\beta$; note that we restrict
ourselves here to time-translation invariant systems).

For the self-energy, additional bath degrees of freedom (denoted $\bm{s}$) are
added to the physical orbitals (denoted $\bm{i}$, $\bm{j}$; all orbitals
denoted $\bm{x}$, $\bm{y}$) in order to emulate the retardation effect (see
Ref. [53] for details),

$\displaystyle
H_{\text{eff}}(t)=\sum_{\bm{x}\bm{y}}h_{\bm{x}\bm{y}}(t)c_{\bm{i}}^{\dagger}c_{\bm{j}}$
(18)

$H_{\text{eff}}$ is determined uniquely by making the $\bm{s}\times\bm{s}$
virtual sector diagonal.

The explicit construction of the Lehmann representation of $\Sigma(z)$ for
$H(t)$ is as follows [53],

$\displaystyle\Sigma_{\bm{i}\bm{j}}(t,t^{\prime})$
$\displaystyle=\delta_{C}(t,t^{\prime})\Sigma_{\bm{i}\bm{j}}^{\text{HF}}(t)+\Sigma_{\bm{i}\bm{j}}^{C}(t,t^{\prime})$
(19) $\displaystyle\Sigma_{\bm{i}\bm{j}}^{\text{HF}}(t)$ $\displaystyle\equiv
2\sum_{\bm{i}^{\prime}\bm{j}^{\prime}}U_{\bm{i}\bm{i}^{\prime}\bm{j}\bm{j}^{\prime}}(t)\langle\mathcal{T}_{C}\
c_{\bm{i}^{\prime}}^{\dagger}(t)c_{\bm{j}^{\prime}}(t)\rangle_{H_{\text{eff}}}$
(20) $\displaystyle\Sigma_{\bm{i}\bm{j}}^{C}(t,t^{\prime})$
$\displaystyle\equiv\sum_{\bm{s}}h_{\bm{i}\bm{s}}(t)g(h_{\bm{s}\bm{s}};t,t^{\prime})h_{\bm{j}\bm{s}}^{*}(t^{\prime}).$
(21)

where the correlated $\Sigma^{C}$ term is the self-energy of the effective
model;
$g(\epsilon;t,t^{\prime})=i[1/(e^{\beta\epsilon}+1)-\Theta_{C}(t,t^{\prime})]e^{i\epsilon(t-t^{\prime})}$
[$\Theta_{C}(t,t^{\prime})=1$ for $t\geqslant_{C}t^{\prime}$,
$\Theta_{C}(t,t^{\prime})=0$ otherwise] is the non-interacting Green’s
function of an isolated one-particle mode ($h_{\text{mode}}=\epsilon
c^{\dagger}c$) with excitation energy $\epsilon$.

We aim to prove that for $z\in\mathcal{C}^{+}$,
$i\Sigma(z)+(i\Sigma(z))^{\dagger}$ is a PSD matrix, i.e., $i\Sigma(z)$ is a
Carathéodory function on $\mathcal{C}^{+}$. Fourier transforming $\Sigma^{C}$
from time to frequency yields

$\displaystyle\Sigma^{C}_{\bm{i}\bm{j}}(i\omega_{n})$
$\displaystyle=\int_{0}^{\beta}-i\Sigma^{C}_{\bm{i}\bm{j}}(-i\tau,0)e^{i\omega_{n}\tau}d\tau$
(22)
$\displaystyle=\int_{0}^{\beta}\sum_{\bm{s}}h_{\bm{i}\bm{s}}(0)h^{*}_{\bm{j}\bm{s}}(0)\frac{-e^{(\beta+\tau)h_{\bm{s}\bm{s}}}}{e^{\beta
h_{\bm{s}\bm{s}}}+1}e^{i\omega_{n}\tau}d\tau$ (23)
$\displaystyle=\sum_{\bm{s}}\frac{h_{\bm{i}\bm{s}}(0)h^{*}_{\bm{j}\bm{s}}(0)(e^{i\omega_{n}\beta}-e^{h_{\bm{s}\bm{s}}\beta})}{(1+e^{h_{\bm{s}\bm{s}}\beta})(h_{\bm{s}\bm{s}}-i\omega_{n})}$
(24)
$\displaystyle=\sum_{\bm{s}}\frac{h_{\bm{i}\bm{s}}(0)h^{*}_{\bm{j}\bm{s}}(0)}{i\omega_{n}-h_{\bm{s}\bm{s}}}$
(25) $\displaystyle\Sigma^{C}_{\bm{i}\bm{j}}(z)$
$\displaystyle=\sum_{\bm{s}}\frac{h_{\bm{i}\bm{s}}(0)h^{*}_{\bm{j}\bm{s}}(0)}{z-h_{\bm{s}\bm{s}}}$
(26)

where $\omega_{n}=\frac{(2n+1)\pi}{\beta}$ are the fermionic Matsubara
frequencies so that $e^{i\omega_{n}\beta}=-1$. Observe that
$\Sigma^{\text{HF}}(z)$ is Hermitian and $z$-independent and that self-energy
has the property $\Sigma(x+yi)=(\Sigma(x-yi))^{\dagger}$ for $x,y>0$, for any
complex vector $\rvert x\rangle$,

$\displaystyle\langle x\lvert i\Sigma(z)+(i\Sigma(z))^{\dagger}\rvert
x\rangle$ (27) $\displaystyle=$ $\displaystyle\langle x\lvert
i\Sigma^{C}(x+yi)-i\Sigma^{C}(x-yi)\rvert x\rangle$ (28) $\displaystyle=$
$\displaystyle\sum_{\bm{s}}\frac{2y\sum_{\bm{i}\bm{j}}x_{\bm{i}}x^{*}_{\bm{j}}h_{\bm{i}\bm{s}}(0)h^{*}_{\bm{j}\bm{s}}(0)}{(h_{\bm{s}\bm{s}}-x)^{2}+y^{2}}$
(29) $\displaystyle=$
$\displaystyle\sum_{\bm{s}}\frac{2y(\sum_{\bm{i}}x_{\bm{i}}h_{\bm{i}\bm{s}}(0))^{2}}{(h_{\bm{s}\bm{s}}-x)^{2}+y^{2}}\geqslant
0$ (30)

$i\Sigma(z)$ is therefore a Carathéodory function on $\mathcal{C}^{+}$.

####  II.4.3 $iM(z)$ is a Carathéodory function

The Cumulant $M$ is a derived object that has the properties of a Green’s
function but lacks the ‘Fock’ contribution of the one-body Hamiltonian [40,
54]. It has the definition of

$\displaystyle M^{-1}(z)$ $\displaystyle=G^{-1}(z)+F$ (31)
$\displaystyle=(z+\mu)I-\Sigma(z)$ (32)

where $F$ is the Fock matrix. It is related to the Green’s function as

$\displaystyle-iG^{-1}(z)+i(G^{-1}(z))^{\dagger}$ (33) $\displaystyle=$
$\displaystyle-iM^{-1}(z)+iF+i(M^{-1}(z))^{\dagger}-iF^{\dagger}$ (34)
$\displaystyle=$ $\displaystyle-iM^{-1}(z)+i(M^{-1}(z))^{\dagger}$ (35)

Since $F$ is a constant Hermitian matrix, and the inverse of a Carathéodory
function (if it exists) is a Carathéodory function [55], $iM(z)$ is a matrix-
valued Carathéodory function on $\mathcal{C}^{+}$.

Figure 3: Total spectral function of the Hubbard dimer, obtained in the site
basis and in a randomly rotated basis, illustrating the basis independence of
the continuation procedure. Figure 4:  Diagonal and off-diagonal entries of
Green’s functions, self-energies, and cumulants of a Hubbard dimer (see text
for parameters). Dashed lines: exact data evaluated at $\omega+i\eta$ for
$\eta=0.01$. Solid lines: analytically continued data evaluated at
$\eta=0.01$. Inset and red line: Maximum entropy continuation of $G_{22}$. Top
Row: Green’s function. Middle Row: Self-energy. Bottom Row: Cumulant. First
and second column: diagonal and off-diagonal elements obtained in site-basis.
Third column: Off-diagonal entry of data rotated by a random unitary matrix.

##  III Results

###  III.1 Hubbard Dimer

We now show results for analytical continuations. As a first example we choose
the Hubbard dimer, which consists of two fermionic spin $1/2$ sites (‘0’ and
‘1’) which may be occupied by a particle with either spin $\uparrow$ or
$\downarrow$. The Fock space of the system contains sixteen configurations,
and exact solutions for Green’s functions or self-energies can readily be
obtained.

The Hamiltonian of the system is

$\displaystyle H=H_{0}+H_{V}+H_{\bm{H}}+H_{SB},$ (36)

where
$H_{0}=-\sum_{\sigma}t(c^{\dagger}_{0\sigma}c_{1\sigma}+h.c.)-\sum_{i,\sigma}\mu
n_{i\sigma}$ describes the usual quadratic spin-diagonal hopping and chemical
potential terms and
$H_{V}=\sum_{i}UD_{i}-\sum_{i,\sigma}\frac{U}{2}n_{i\sigma}$ a Hubbard-type
[56] on-site interaction, with $D_{i}$ denoting the double occupancy on site
$i$. In addition to these terms, we add a magnetic field and a symmetry-
breaking term $H_{\bm{H}}=\sum_{i}\bm{H}(n_{i\uparrow}-n_{i\downarrow})$ and
$H_{SB}=U_{a}(D_{0}-D_{1})+\mu_{a}(n_{0\uparrow}+n_{0\downarrow}-n_{1\uparrow}-n_{1\downarrow})+\bm{H}_{a}(n_{0\uparrow}-n_{0\downarrow}-n_{1\uparrow}+n_{1\downarrow})$,
with the aim of breaking degeneracies and thereby producing a spectral
function with additional features. For concreteness, we choose $\beta=10$,
$t=1$, $U=5$, $\mu=0.7$, $\bm{H}=0.3$, $U_{a}=0.5$, $\mu_{a}=0.2$,
$\bm{H}_{a}=0.03$ and the summations are over $i=0,1$ and
$\sigma=\uparrow,\downarrow$. All Green’s function results are evaluated on
the Matsubara axis and stored in double precision as an input to the
interpolation algorithm.

Fig. 3 shows the spectral function of the system. Fig. 4 shows the real
(green) and imaginary (purple) parts of the Green’s function (top row), self-
energy (middle row), and cumulant (bottom row) as a function of frequency.
Shown are a sample diagonal (22) element (left column), a sample offdiagonal
(31) element (middle column), and an off-diagonal element that has been
obtained after rotating the fermion operators by a random unitary rotation in
orbital space (right column).

All dashed lines were obtained exactly by diagonalizing the system in Fock
space, obtaining all eigenvalues and eigenvectors, and subsequently computing
the real-frequency Green’s function via the Lehmann representation. Once the
Green’s function was known, the process was repeated for the non-interacting
system and self-energies and cumulants were obtained by inverting the Dyson
equation in real frequencies.

All colored (green, purple) lines were obtained using the algorithm described
in Sec. II. For the Hubbard dimer the results both for the imaginary and the
real parts match the exact results precisely, showing the success of the
Carathéodory continuation method. Data was evaluated on the Matsubara axis on
a non-equidistant intermediate representation [57, 58] grid with $36$ positive
Matsubara frequency points, interpolated, and evaluated just above the real
frequency axis at $\omega+i\eta$ for $\eta=0.01$.

The two top left panels illustrate existing capabilities: diagonal entries of
the Green’s function are strictly negative, as Green’s functions are
Nevanlinna [39]. Standard continuation methods such as the Maximum Entropy
method [11, 16] obtain a broadened version of these spectra (red line in the
inset of Fig. 4), from which the real parts can be obtained via a Kramers
Kronig relation. If only spectral functions are desired, diagonal components
of the Green’s functions are sufficient.

The top middle panel exhibits the new capability developed in this work: the
analytic continuation of off-diagonal matrix elements. Shown is $G_{13}$,
which is one of many off-diagonal components. As is evident, the imaginary
part contains both positive and negative entries, and such Green’s functions
can therefore not be obtained using standard Maximum Entropy [11]. As with the
diagonal part, real and imaginary parts are related by a Kramers Kronig
transform.

Figure 5: Total spectral function of Si evaluated along a high-symmetry path
obtained within three approaches, evaluated just above real axis
($\omega+i\eta$ for $\eta=0.01$). Left: continuation of the Green’s function
$G$. Middle: continuation of the self-energy $\Sigma$ followed by Dyson
equation at $\omega+i\eta$. Right: continuation of the cumulant $M$ and
application of the Dyson equation at $\omega+i\eta$.

The top right panel shows a particularly difficult case: an orbital rotation
with a random unitary matrix that mixes all Green’s function components. While
individual elements of the Green’s functions are mixed, the rotation preserves
the positive definiteness of the Green’s function and therefore its
Carathéodory character. The total spectral function, which is basis
independent in theory, remains basis independent in practice, giving a first
indication that the precision of the continuation and its mathematical
properties may be sufficient to perform algebraic operations on the real axis.

To explicitly show that a random orbital rotation leaves the spectral function
invariant, Fig. 3 shows the total spectral function of this example, both
obtained in the original site basis and in the randomly rotated basis. Peak
locations as well as peak heights of the continuations agree after the
rotation.

We now turn our attention to the middle row of Fig. 4. Shown are continuations
(real and imaginary parts) of the self-energy. The middle left panel
illustrates the present capabilities: similar to Green’s functions, diagonal
entries of the self-energy are Carathéodory, and can therefore be continued
with established methods [59]. Real parts can then be obtained via a Kramers
Kronig transform. However, knowledge of the diagonal parts is not sufficient
for obtaining the spectral function, as the inversion in the Dyson equation
mixes diagonal and off-diagonal elements. Similar to the continuation of
Green’s function, no deviation between the exact answer and the continued
spectral function is visible.

The middle panel shows off-diagonal self-energies, which exhibit both positive
and negative contributions. The right panel contains the self-energies in a
randomly rotated basis chosen to maximize off-diagonal contributions. The fact
that off-diagonal entries are accessible now allows one to perform arithmetic
operations such as the inversion of the Dyson equation on the real axis.

Finally, the bottom three panels of Fig. 4 show the cumulant [40] for this
case. The cumulant has the same units as the Green’s function and can be
thought of as the Green’s function of an interacting system ‘without its band
structure contribution’. Given the knowledge of the Fock matrix and the
matrix-valued cumulant on the real axis, the spectral functions and self-
energies can then be obtained via inversion. As in the case of the self-
energies and Green’s functions, all structure of the system is recovered.

###  III.2 Realistic example: Silicon

Next, we present results for a realistic ab-initio simulation. As an example
we show data for crystalline Si, which is a weakly correlated system for which
the band structure is well known. The system is solved within the fully self-
consistent GW approximation in Gaussian orbitals (gth-dzvp-molopt-sr basis
[60] and gth-pbe pseudopotential [61]), using integrals generated by the pySCF
package [62]. Calculations are performed on a $6\times 6\times 6$ grid [63]
with $26$ orbitals per unit cell and interpolated using a Wannier
interpolation of Matsubara data [64] on 52 IR positive frequencies [57, 58] on
a grid with $200$ $k$-points along a high-symmetry path. The total spectral
function of the system, obtained by taking the trace of the $26$ elements of
the matrix-valued Carathéodory interpolation of Sec. II, is shown in the left
panel of Fig. 5. The broadening parameter of $\eta=0.01$ Ha $\sim 3157K$ used
for all continuations here is larger than the temperature of $T=0.001$ Ha
$\sim 316K$ and correlation effects. Due to the precision of the continuation,
individual bands, band gap, and degeneracies at high-symmetry points and
directions are clearly visible. (The indirect band gap of $\sim 2$ eV differs
from the experimental band gap of 1.1eV due to the self-consistent GW
approximation, basis set effects, and finite size effects). High- and low-
frequency values near $\pm$ 10 eV show continuation artifacts. We note that
continuations for the spectral function can be obtained directly from the
diagonal elements of the Green’s functions via Nevanlinna continuation [39],
and do not require continuation of the off-diagonal elements.

The middle panel shows the spectral function obtained from a continuation of
the self-energy with all off-diagonal elements, followed by the Dyson equation
on the real axis (see Fig. 2). No deviation from the Green’s function data is
visible, illustrating the precision of the self-energy continuation. As in the
left panel, bands, degeneracies, and band gaps are clearly visible.

Finally, we show the spectral function obtained from the cumulant in the right
panel. In realistic systems, cumulants are more convenient to work with than
self-energies and Green’s functions, as they have the same units as Green’s
functions, and as (due to the elimination of the Fock term) most of their
structure lies at low energies, thereby increasing the precision of arithmetic
operations and analytical continuations. No deviation from the Green’s
function data is visible, illustrating that the numerical solution of the
Dyson equation is possible on the real axis.

Figure 6: Eigenvalues of the Pick matrices, Eqs. 3 (blue) and 4 (red). Top
panel: scalar (Nevanlinna [39]) continuation of the Green’s function of Si.
Bottom panel: Eigenvalues of the generalized Pick matrix of the matrix-valued
problem. Figure 7: Effect of the commonly used diagonal approximation of the
self-energy on the spectral function. Left panel: Spectral function obtained
by continuing the self-energy including all off-diagonal parts; Middle panel:
Spectral function obtained by omitting the dynamical self-energy altogether;
Right panel: Spectral function obtained by continuing only the diagonal parts
of the dynamical self-energy using the method of Ref. [39].

Carathéodory continuations are sensitive to the precision of the input data.
In practice, for large systems, we find that some eigenvalues of the Pick
matrix dip below zero due to finite double precision and roundoff errors, and
that therefore no Carathéodory interpolant exists on the unit disk. This is
illustrated by Fig. 6, which shows the Pick matrix eigenvalues for the Green’s
function continuations of Fig. 5. However, even in the case where the Pick
criterion is violated, interpolations can still be performed, and the
evaluation just above the real axis (in our case shifted by $\eta=0.01$) often
leads to positive definite spectral functions. Negative contributions, where
they occur, lie very close to the real axis and in frequency regions above the
$10$ eV shown here.

###  III.3 Diagonal Approximation of the Self-energy

It is common practice to perform continuations only of the diagonal parts of
the self-energy or cumulant, truncating all off-diagonal contributions. In
systems where the self-energy can be chosen to be diagonal for all frequencies
due to symmetry, such as in rotationally invariant three-orbital systems [65,
66, 67] and certain model systems, this is exact. However, self-energies are
in general off-diagonal, and can only be diagonalized for one Matsubara
frequency at a time. Performing a diagonal truncation of the self-energy is an
uncontrolled approximation.

To illustrate the effect of the off-diagonal entries of the self-energy on the
spectral function, we show in Fig. 7 results from continuing the self-energy
with all entries along with results obtained by only continuing the diagonal
part of the dynamical self-energy and results that completely omit the
dynamical part of the self-energy. While overall features remain robust, the
resulting bands of the diagonal and mean field continuation are shifted and
distorted with respect to the full matrix continuation, indicating that a
correct treatment of the off-diagonal matrix elements of the self-energy is
essential for precise results [38], even in systems that are only weakly
correlated.

##  IV Conclusions

In conclusion, we have shown that matrix-valued Green’s functions, self-
energies, and cumulants are related to Carathéodory and Schur functions. With
this analytic structure, we were able to demonstrate the conditions under
which a causal positive semi-definite interpolation exists, to construct an
interpolation algorithm, and to parametrize all causal interpolations. An
application to a simple benchmark problem showed that diagonal and off-
diagonal parts of the self-energies, cumulants, and Green’s functions could be
recovered to high precision. A demonstration for a real materials multi-
orbital problem showed that the analytic continuation commutes with the
application of the Dyson equation, and that a diagonal approximation of the
self-energy did not produce satisfactory results.

Knowledge of the off-diagonal elements of the real-frequency Green’s
functions, self-energies, and cumulants is a prerequisite for the evaluation
of the equations of many-body theory on the real axis including the real-
frequency calculation of susceptibilities and other response functions, and
the analysis and inversion of self-energies.

In systems where precise input data is available, our method provides accurate
continuations of off-diagonal terms. In system where substantial noise is
present, such as in those solved by multi-orbital QMC impurity solvers [5], a
positive definite interpolation does not exist. Causal projections or
generalizations of the Maximum entropy fitting procedure [37, 38] may then be
used instead. The derivation of such methods, in addition to the development
of methods for the continuation of bosonic quantities, is an interesting
problem for future investigation.

###### Acknowledgements.

EG and FJ were supported by the Simons Foundation via the Simons Collaboration
on the Many-Electron Problem. DZ and CNY were supported by the U.S. Department
of Energy under Award No. DE-SC0019374.

## References

  * Hedin [1965] L. Hedin, Phys. Rev. 139, A796 (1965).
  * Dahlen and van Leeuwen [2005] N. E. Dahlen and R. van Leeuwen, The Journal of Chemical Physics 122, 164102 (2005).
  * Phillips and Zgid [2014] J. J. Phillips and D. Zgid, The Journal of Chemical Physics 140, 241101 (2014).
  * Blankenbecler  _et al._ [1981] R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Phys. Rev. D 24, 2278 (1981).
  * Gull  _et al._ [2011] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Rev. Mod. Phys. 83, 349 (2011).
  * Asakawa  _et al._ [2001] M. Asakawa, Y. Nakahara, and T. Hatsuda, Progress in Particle and Nuclear Physics 46, 459 (2001).
  * Tripolt  _et al._ [2019] R.-A. Tripolt, P. Gubler, M. Ulybyshev, and L. von Smekal, Computer Physics Communications 237, 129 (2019).
  * Rothkopf [2020] A. Rothkopf, Data 5, 10.3390/data5030085 (2020).
  * Bryan [1990] R. K. Bryan, European Biophysics Journal 18, 165 (1990).
  * Creffield  _et al._ [1995] C. E. Creffield, E. G. Klepfish, E. R. Pike, and S. Sarkar, Phys. Rev. Lett. 75, 517 (1995).
  * Jarrell and Gubernatis [1996] M. Jarrell and J. Gubernatis, Physics Reports 269, 133 (1996).
  * Beach [2004] K. S. D. Beach, Identifying the maximum entropy method as a special limit of stochastic analytic continuation (2004), arXiv:cond-mat/0403055 [cond-mat.str-el] .
  * Gunnarsson  _et al._ [2010a] O. Gunnarsson, M. W. Haverkort, and G. Sangiovanni, Phys. Rev. B 82, 165125 (2010a).
  * Gunnarsson  _et al._ [2010b] O. Gunnarsson, M. W. Haverkort, and G. Sangiovanni, Phys. Rev. B 81, 155107 (2010b).
  * Bergeron and Tremblay [2016] D. Bergeron and A.-M. S. Tremblay, Phys. Rev. E 94, 023303 (2016).
  * Levy  _et al._ [2017] R. Levy, J. LeBlanc, and E. Gull, Computer Physics Communications 215, 149 (2017).
  * Gaenko  _et al._ [2017] A. Gaenko, A. Antipov, G. Carcassi, T. Chen, X. Chen, Q. Dong, L. Gamper, J. Gukelberger, R. Igarashi, S. Iskakov, M. Könz, J. LeBlanc, R. Levy, P. Ma, J. Paki, H. Shinaoka, S. Todo, M. Troyer, and E. Gull, Computer Physics Communications 213, 235 (2017).
  * Rumetshofer  _et al._ [2019] M. Rumetshofer, D. Bauernfeind, and W. von der Linden, Phys. Rev. B 100, 075137 (2019).
  * Sandvik [1998] A. W. Sandvik, Phys. Rev. B 57, 10287 (1998).
  * Mishchenko  _et al._ [2000] A. S. Mishchenko, N. V. Prokof’ev, A. Sakamoto, and B. V. Svistunov, Phys. Rev. B 62, 6317 (2000).
  * Vafayi and Gunnarsson [2007] K. Vafayi and O. Gunnarsson, Phys. Rev. B 76, 035115 (2007).
  * Fuchs  _et al._ [2010] S. Fuchs, M. Jarrell, and T. Pruschke, Journal of Physics: Conference Series 200, 012041 (2010).
  * Goulko  _et al._ [2017] O. Goulko, A. S. Mishchenko, L. Pollet, N. Prokof’ev, and B. Svistunov, Phys. Rev. B 95, 014102 (2017).
  * Krivenko and Harland [2019] I. Krivenko and M. Harland, Computer Physics Communications 239, 166 (2019).
  * Otsuki  _et al._ [2017] J. Otsuki, M. Ohzeki, H. Shinaoka, and K. Yoshimi, Phys. Rev. E 95, 061302 (2017).
  * Otsuki  _et al._ [2020] J. Otsuki, M. Ohzeki, H. Shinaoka, and K. Yoshimi, Journal of the Physical Society of Japan 89, 012001 (2020).
  * Yoon  _et al._ [2018] H. Yoon, J.-H. Sim, and M. J. Han, Phys. Rev. B 98, 245101 (2018).
  * Vidberg and Serene [1977] H. J. Vidberg and J. W. Serene, Journal of Low Temperature Physics 29, 179 (1977).
  * Beach  _et al._ [2000] K. S. D. Beach, R. J. Gooding, and F. Marsiglio, Phys. Rev. B 61, 5147 (2000).
  * Östlin  _et al._ [2012] A. Östlin, L. Chioncel, and L. Vitos, Phys. Rev. B 86, 235107 (2012).
  * Osolin, Žiga and Žitko [2013] Osolin, Žiga and R. Žitko, Phys. Rev. B 87, 245135 (2013).
  * Schött  _et al._ [2016] J. Schött, I. L. M. Locht, E. Lundin, O. Grånäs, O. Eriksson, and I. Di Marco, Phys. Rev. B 93, 075104 (2016).
  * Han  _et al._ [2017] X.-J. Han, H.-J. Liao, H.-D. Xie, R.-Z. Huang, Z.-Y. Meng, and T. Xiang, Chinese Physics Letters 34, 077102 (2017).
  * 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, Phys. Rev. B 90, 125114 (2014).
  * Gull and Millis [2014] E. Gull and A. J. Millis, Phys. Rev. B 90, 041110 (2014).
  * Kraberger  _et al._ [2017] G. J. Kraberger, R. Triebl, M. Zingl, and M. Aichhorn, Phys. Rev. B 96, 155128 (2017).
  * Sim and Han [2018] J.-H. Sim and M. J. Han, Phys. Rev. B 98, 205102 (2018).
  * Fei  _et al._ [2021] J. Fei, C.-N. Yeh, and E. Gull, Phys. Rev. Lett. 126, 056402 (2021).
  * Stanescu and Kotliar [2006] T. D. Stanescu and G. Kotliar, Phys. Rev. B 74, 125110 (2006).
  * Carathéodory [1907] C. Carathéodory, Mathematische Annalen 64, 95 (1907).
  * Schur [1918] J. Schur, Journal für die reine und angewandte Mathematik 148, 122 (1918).
  * Delsarte  _et al._ [1979] P. Delsarte, Y. Genin, and Y. Kamp, SIAM Journal on Applied Mathematics 36, 47 (1979).
  * Arov and Rozhenko [2008] D. Arov and N. Rozhenko, St Petersburg Mathematical Journal 19, 375 (2008).
  * Fritzsche  _et al._ [2012] B. Fritzsche, B. Kirstein, A. Lasarow, and A. Rahn, On reciprocal sequences of matricial carathéodory sequences and associated matrix functions, in  _Interpolation, Schur Functions and Moment Problems II_, edited by D. Alpay and B. Kirstein (Springer Basel, Basel, 2012) pp. 57–115.
  * Pick [1917] G. Pick, Mathematische annalen 78, 270 (1917).
  * Akhiezer [1965] N. Akhiezer,  _The Classical Moment Problem and Some Related Questions in Analysis_ (Oliver and Boyd., London, 1965).
  * Khargonekar and Tannenbaum [1917] P. Khargonekar and A. Tannenbaum, IEEE Transactions on Automatic Control 30, 1005 (1917).
  * Chen and Çetin Kaya Koç [1994] G. Chen and Çetin Kaya Koç, Linear Algebra and its Applications 203-204, 253 (1994).
  * Yazici and Sevindir [2013] C. Yazici and H. K. Sevindir, AIP Conference Proceedings 1558, 2474 (2013).
  * Potapov [1955] V. Potapov,  _The multiplicative structure of J-contractive matrix functions_ (GITTL, Moscow, 1955) pp. 125 – 236.
  * Balzer and Eckstein [2014] K. Balzer and M. Eckstein, Phys. Rev. B 89, 035148 (2014).
  * Gramsch and Potthoff [2015] C. Gramsch and M. Potthoff, Phys. Rev. B 92, 235135 (2015).
  * Sakai  _et al._ [2012] S. Sakai, G. Sangiovanni, M. Civelli, Y. Motome, K. Held, and M. Imada, Phys. Rev. B 85, 035102 (2012).
  * Johnson [1972] C. R. Johnson,  _Matrices whose hermitian part is positive definite_ , Ph.D. thesis (1972).
  * Qin  _et al._ [2021] M. Qin, T. Schäfer, S. Andergassen, P. Corboz, and E. Gull, The hubbard model: A computational perspective (2021), arXiv:2104.00064 [cond-mat.str-el] .
  * 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).
  * Wang  _et al._ [2009] X. Wang, E. Gull, L. de’ Medici, M. Capone, and A. J. Millis, Phys. Rev. B 80, 045101 (2009).
  * VandeVondele and Hutter [2007] J. VandeVondele and J. Hutter, The Journal of Chemical Physics 127, 114105 (2007).
  * 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, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov, and G. K.-L. Chan, The Journal of Chemical Physics 153, 024109 (2020).
  * Iskakov  _et al._ [2020] S. Iskakov, C.-N. Yeh, E. Gull, and D. Zgid, Phys. Rev. B 102, 085105 (2020).
  * [64] C. Yeh, Y. Chen, J. Fei, E. Gull, and D. Zgid, Wannier interpolation of the self-energy in real materials systems, in preparation.
  * Werner  _et al._ [2009] P. Werner, E. Gull, and A. J. Millis, Phys. Rev. B 79, 115119 (2009).
  * Georges  _et al._ [2013] A. Georges, L. d. Medici, and J. Mravlje, Annual Review of Condensed Matter Physics 4, 137 (2013).
  * Yin  _et al._ [2011] Z. P. Yin, K. Haule, and G. Kotliar, Nature Materials 10, 932 (2011).


# Score-informed syllable segmentation for a cappella singing voice with
convolutional neural networks

###### Abstract

This paper introduces a new score-informed method for the segmentation of
jingju a cappella singing phrase into syllables. The proposed method estimates
the most likely sequence of syllable boundaries given the estimated syllable
onset detection function (ODF) and its score. Throughout the paper, we first
examine the jingju syllables structure and propose a definition of the term
“syllable onset”. Then, we identify which are the challenges that jingju a
cappella singing poses. Further, we investigate how to improve the syllable
ODF estimation with convolutional neural networks (CNNs). We propose a novel
CNN architecture that allows to efficiently capture different time-frequency
scales for estimating syllable onsets. In addition, we propose using a score-
informed Viterbi algorithm –instead of thresholding the onset function–,
because the available musical knowledge we have (the score) can be used to
inform the Viterbi algorithm in order to overcome the identified challenges.
The proposed method outperforms the state-of-the-art in syllable segmentation
for jingju a cappella singing. We further provide an analysis of the
segmentation errors which points possible research directions.

##  1 Introduction

The ultimate goal of our research project is to automatically evaluate the
jingju a cappella singing of a student in the scenario of jingju singing
education – see Figure 1. Jingju, a traditional Chinese performing art form
also known as Peking or Beijing opera, is extremely demanding in the clear
pronunciation and accurate intonation for each syllabic or phonetic singing
unit. To this end, during the initial learning stages, students are required
to completely imitate tutor’s singing. Therefore, the automatic jingju singing
evaluation tool we envision is based on this training principle and measures
the intonation and pronunciation similarities between the student’s and the
tutor’s singings. Before measuring the similarities, the singing phrase should
be automatically segmented into syllabic or phonetic units in order to capture
the temporal details.

In this paper we tackle the problem of score-informed automatic syllable
segmentation for a cappella singing (bold rectangle in Figure 1). Jingju music
scores, which contain the duration information for each singing syllable, will
be a helpful hint for the segmentation method.

Figure 1: Framework of the entire research project. The module with bold
border is addressed in this paper. Figure 2: A dan role-type singing phrase
with its last syllable prolongated. Top: log-mel bands spectrogram. Bottom:
syllable ODF estimated with Schlüter’s method [16] (blue), and ground truth
syllable onsets (red).

The syllable segmentation task consists on determining the time positions of
the syllable boundaries – onset and offset. In this research, we consider the
onset of the subsequent syllable to be the offset of the current one.
Therefore, we treat the segmentation task as an onset detection problem.

Most segmentation methods rely on first estimating an ODF. For example,
Klapuri [7] utilizes a band-wise processing principle inspired by
psychoacoustics. He divides the audio signal into 21 non-overlapping bands,
then detects onset components in each band ODF and finally combines them to
yield onsets. Obin et al. [10] obtain a syllable ODF by fusing mel-frequency
intensity profiles and voicing profiles. One shortcoming of these band-wise
methods has been already pointed out by Klapuri [7]: “they are unable to deal
with strong amplitude modulations” – what is very common in singing voice
recordings.

On the other hand, some methods are based on features and supervised learning.
For example, Toh et al. [17] use two Gaussian Mixture Models (GMMs) to
classify singing audio frames into onset or non-onset classes. Three timbral
features (MFCCs, LPCCs and ERB-bands) are chosen as input. Their results show
that this supervised method is superior to many band-wise ODF-based methods.

Neural networks have also been successfully explored for the task of musical
onset detection. Eyben et al. [4] trained a bidirectional long-short term
memory neural network on mel-scale magnitude spectrograms to estimate an ODF.
Schlüter et al. [16] proposed to use CNNs to estimate the ODF, which defines
the current state-of-the-art for onset detection. The advantage of applying
deep learning methods in musical onset detection is that one does not need to
design rules or handcraft features to capture the relevant facets for
detecting onset/non-onset frames – which are very difficult to design.
Besides, if context (more than one frame) is input into the network, spectro-
temporal features can also be learned – what might be useful to learn slow-
transient features defining some onsets.

Thresholding or peak-picking operations are normally executed over the ODF in
order to determine the final onsets [7, 17, 4, 16]. However, we will argue
that these operations are not suitable for selecting jingju singing syllable
onsets. Probabilistic models –which allow incorporating prior domain knowledge
for decision making– usually result in a performance gain compared to simple
thresholding and peak-picking methods. For example, Böck et al. [1] tracked
the beat/downbeat by using a dynamic bayesian network (DBN) observing a beat
ODF estimated by a recurrent neural network (RNN). Or Obin et al. [10] applied
a segmental Viterbi algorithm over a syllable ODF to detect speech syllable
onsets. These two approaches are closely related to ours: former one relates
with our work because it uses deep learning to estimate the ODF, and latter
one because we also consider using the Viterbi algorithm for estimating the
onset candidates.

This paper introduces a new score-informed method for the segmentation of
jingju singing into syllables. We first define what a “syllable” and “syllable
onset” is in the context of jingju music. By doing so, in section 2 we
introduce the challenges we aim to address. The audio dataset we use is
described in section 3. Section 4 explains the CNN architecture used for
estimating the syllable ODF, and the Viterbi decoding algorithm that exploits
the prior syllable duration information extracted from the score. Evaluation
and error analysis are conducted in section 5, and section 6 concludes this
work.

##  2 Background and Challenges

Jingju singing is the most precise articulated rendition of the spoken
Mandarin language. Although certain special pronunciations in jingju
theatrical language differ from their normal Mandarin pronunciations –due to:
firstly, the adoption of certain regional dialects; and secondly, the ease or
variety in pronunciation and projection of sound– the mono-syllabic
pronouncing structure of the standard Mandarin doesn’t change [18].

A syllable/character of jingju singing is composed of three distinct parts in
most of the cases: the “head” (tou), the “belly” (fu) and the “tail” (wei).
The “head” consists of an initial consonant or semi-vowel – and a medial
vowel, if the syllable includes one. The “head” is not normally prolonged in
its pronunciation except when there is a medial vowel. The “belly” follows the
“head” and consists of the central vowel and it is normally prolonged. The
“belly” is the most sonorous part of a jingju singing syllable and can be
analogous to the nuclei of a speech syllable. The “tail” is composed of the
terminal vowel or consonant [18]. However, there are syllables in jingju
singing where “head” is absent – only “belly” is an obligatory element. To
avoid ambiguity, we define the syllable onset as: “the start of the initial
consonant if the syllable includes one or the start of the central vowel
otherwise”. This definition is slightly different from that agreed by most of
the phonological theories [5]: “any consonants that precede the nuclear
element (the vowel)”. It is also worth stressing that our notion of singing
syllable onset is different from the singing onset defined in Toh et al.’s
paper [17]: “the start of a new human-perceived note, taking into account
contextual cues”, of which the latter emphasizes on the intonational aspect
instead of the phonological aspect of the singing voice.

###  2.1 Challenges

Figure 2 shows an example of a dan role-type singing phrase in which the last
syllable lasts approximately 20s. This singing method is more common in dan
singing than in laosheng singing. However, both role-types use it as a way of
improving artistic expression and showing off their singing skills. These
skills include breathing techniques, intonational techniques and dynamic
control techniques, among others. The syllable ODF (blue curve in Figure 2) is
generated using a CNN model based on Schlüter et al.’s work [16], which is
considered the state-of-the-art. This CNN model is trained with the jingju
dataset presented in section 3. The resulting syllable ODF is quite robust to
pitch variations and continuous dynamic change since no prominent peaks are
observed in these regions. However, numerous peaks can be found in the start
and end positions of each articulation throughout this long syllable
($\approx$ 20s). Note that thresholding or peak-picking would choose these
peaks (false syllable onsets) since they have similar amplitude level than
real onsets. For that reason, we propose using a score-informed decoding
method as an alternative to naive thresholding.

Stealing breath (tou qi) is one of the major methods of taking breath in
jingju singing. A breath is stolen when a sound is too long to be delivered in
one breath – and no vocal pauses are desired [18]. However, this is not the
only technique which can lead to pauses within a syllable. Another singing
technique (zu yin – literal translation: block sound), provokes also pauses
without occurring exhalation or inhalation. This kind of pause can be very
short in duration and can be easily found in jingju singing syllables, see
Figure 3 for an example. These silences can be another source of false
positives if thresholding or peak-picking is used. However, these can also be
avoided by using an score-informed decoding method.

A priori syllable duration information is often easy to obtain from the score
and this is an advantage which we exploit to avoid previously described
issues. The repertoire of jingju includes around 1400 plays [18], and most
used teaching pieces are transcribed into scores. We use the score to guide
the Viterbi decoding throughout the onset detection process.

Figure 3: Another dan role-type singing phrase. A zu yin (block sound) is
located by the arrow within a normal length syllable. Schlüter’s syllable ODF
(blue curve). Ground truth syllable onset positions (red vertical lines).

##  3 Dataset

The a cappella singing audio dataset111https://goo.gl/y0P7BL used for this
study focuses on two most important jingju role-types [14]: dan (female) and
laosheng (old man). It contains 39 interpretations of 31 unique arias sung by
11 jingju singers. Audio is sampled at 44.1 kHz and is is pre-segmented into
phrases. The syllable onset ground truth is manually annotated in Praat [2] –
with 291 phrases and 2641 syllables (including padding written-characters
[18]). Two Mandarin native speakers and a jingju musicologist have been
devoted to this annotation work. Syllable durations are manually transcribed
from music scores considering as unit duration the quarter note. The whole
dataset is randomly split into training, validation and test sets (60%, 20%
and 20%, respectively). The percentage of presence of each role-type in a
split is kept constant throughout sets.

In order to highlight the differences between jingju music and popular western
music, we compare the vowel’s duration between Kruspe’s dataset [8] (a
cappella singing of commercial pop songs) and ours. In Kruspe’s dataset, the
standard deviation duration of voiced phonemes is of 0.31s, whereas this
duration is more than doubled in our dataset: 0.75s. Since voiced phonemes are
the main component of a syllable, it thus becomes clear that jingju singing
syllable durations show huge variations. Therefore, it is impossible to model
syllable durations with a single distribution. For that reason the proposed
model must be able to handle different syllable lengths depending on the case.
To this end, syllable durations extracted from scores will be a valuable
information.

##  4 Approach

A score-informed syllable boundary detection approach is explained in this
section. We first propose some improvements to the current state-of-the-art
CNN onset detection model proposed by Schlüter et al. The syllable ODF output
from the CNN serves as observation probability for the syllable boundary
decoding process. Then, an a priori syllable duration model based on the score
is proposed. This model guides the Viterbi algorithm by informing it in which
time-positions is likely to occur a syllable boundary. Therefore, the syllable
boundaries sequence is decoded by taking advantage of a score-informed Viterbi
algorithm.

###  4.1 CNN syllable onset detection function

Studied CNNs are inspired by Schlüter et al.’s [16] and Pons et al.’s work
[11, 12]. Schlüter et al. have shown that CNNs fed with spectrograms can
achieve state-of-the-art performance for the onset detection task [16]. On the
other hand, Pons et al. [11] have recently proposed a novel design strategy
for spectrograms-based CNNs. They propose using different filter shapes in the
first layer so that local stationarities in spectrograms (present at different
time/frequency-scales) can be efficiently captured.

Explored models are based on Schlüter et al.’s architecture [16], which
consists of: two convolutional layers, and a dense layer of 256 units
connected to an output softmax layer – with two output units standing for
onset and non-onset. First CNN layer uses 20x $3{\times}7$ filters222First
number denote the number of filters (ie. 20x). Second and third respectively
denote the frequency and temporal size of the filter (ie. $3{\times}7$). and
is followed by a $3{\times}1$ max pool layer. Second CNN layer has 20x
$3{\times}3$ filters and is followed by a $3{\times}1$ max-pool layer. Input
is set to be a log-mel spectrogram333Original Schlüter et al.’s model inputs
three channels with different resolution spectrograms to the network. However,
preliminary results showed that a single channel was performing better than
three. of size $80{\times}21$ – note that the network takes a decision for
every frame given its context: $\pm$10ms, 21 frames in total.

Proposed architectures follow Pons et al. [11, 12] design strategy and several
filter shapes in the first layer are used. This approach allows to efficiently
capture different time-frequency scales, what might be interesting for the
task at hand because onsets can exhibit different time-frequency patterns. For
example, Figure 4 depicts two onsets: middle onset is expressed as an abrupt
time-frequency change corresponding to an onset starting with a consonant, and
right onset is expressed as a gradual timbral change that corresponds to an
onset starting with a central vowel. Note that these examples are
representative of the proposed definition of syllable onset. Moreover, using
different filter shapes in the first layer promotes a much richer
representation out of the first layer. In the following, several CNNs are
designed to efficiently capture the relevant time-frequency contexts and
scales for onset detection.

Figure 4: Log-mel spectrograms (80 bands, 21 frames). left: Non-onset
spectrogram. middle and right: Onset spectrograms. Red line denotes the onset
frame.

Schlüter et al.’s architecture and proposed ones only differ in the way the
first convolutional layer and its following max-pool layer are set up. Input
and remaining layers are kept intact – unless it is explicitly stated:

$\rightarrow$ Temporal architecture. Note that the small filters proposed by
Schlüter et al. might have difficulties on learning the longer time-scales
required to model the gradual changes that some onsets express – see Figure 4,
right. Several filter shapes are designed to efficiently capture different and
longer time-scales [11] than Schlüter et al.:

$\cdot$ 12x $1{\times}7$, 6x $3{\times}7$ and 3x $5{\times}7$

$\cdot$ 12x $1{\times}12$, 6x $3{\times}12$ and 3x $5{\times}12$

These CNNs are followed by a $3{\times}5$ max-pool layer.

$\rightarrow$ Timbral architecture. Since our syllable onset definition
considers phonological aspects, we consider the timbre to be an important
feature for our task. Filter shapes are designed to learn timbral
representations, note that these filters span throughout the frequency domain
[12]:

$\cdot$ 12x $50{\times}1$, 6x $50{\times}5$ and 3x $50{\times}10$

$\cdot$ 12x $70{\times}1$, 6x $70{\times}5$ and 3x $70{\times}10$

These CNNs are followed by a $5{\times}3$ max-pool layer.

We also explore combining temporal and timbral architectures. We propose a
late-fusion approach where we multiply the estimated output probabilities from
temporal and timbral (independently trained) architectures444Using temporal
and timbral filters in first layer yield to worse results..

In addition, we also explore increasing the number of filters in the second
layer from 20x to 32x.

For concatenating several feature maps resulting of CNNs with different filter
shapes, it is required to use zero padding. We apply same padding in the first
CNN layer, so that all resulting feature maps have the same length and these
can be concatenated. STFT was performed using a 25ms window (2048 samples with
zero-padding) with a hop size of 10ms. The 80 log-mel bands energies are
calculated on frequencies between 0Hz and 11000Hz and spectrograms are
standardized to have zero mean and unit variance. We use L2 weight decay
regularization, ELU activation functions [3] and 30% dropout for each layer.
The model parameters are learned with mini-batch training (batch size 128)
using the ADAM update rule [6] and early stopping – if validation loss
(categorical cross-entropy) was not decreasing after 10 epochs. In order to
allow a fair comparison between Schlüter’s and Pons’ architectures, the
original Schlüter architecture is modified to meet above described hyper-
parameters.

Finally, syllable ODFs estimated with CNNs are smoothed by convolving those
with a 5 frames Hanning window [16] – since estimated syllable ODFs are
typically very spiky.

###  4.2 A priori duration model

The a priori duration model is shaped with a Gaussian function
$\mathcal{N}(x;{\mu}_{l},\sigma_{l}^{2})$ whose mean ${\mu}_{l}$ represents
the $l$-th relative syllable duration – according to the score. Its standard
deviation $\sigma_{l}$ is proportional to $\mu_{l}$:
$\sigma_{l}=\gamma\mu_{l}$ and $\gamma$ is heuristically set to 0.35.

$\mathcal{N}(x;\mu_{l},\sigma_{l}^{2})=\frac{1}{\sqrt{2\pi}\sigma_{l}}\exp\left(-\frac{(x-\mu_{l})^{2}}{2\sigma_{l}^{2}}\right).$
(1)

Figure 5 provides an intuitive example of how the a priori duration model
works. Observe that it provides the prior likelihood of an onset to occur
according to the duration in the score.

Figure 5: A priori relative duration distributions (bottom) of the syllables
of a singing phrase.

The relative duration of each note is measured considering a quarter note
length as a unit, so an eighth note has a duration of 0.5. We only keep the
relative duration and discard the tempo information of the score. By
normalizing the summation of the notes’ relative durations to the incoming
audio recording’s duration, we obtain the absolute notes’ durations. Then, a
sequence of syllable absolute durations $M=\mu_{1}\mu_{2}\cdots\mu_{L}$ is
deduced by summing the notes’ absolute durations corresponded with each
syllable, where $L$ is the total syllable number of the score. The a priori
duration model distributions will be incorporated into the Viterbi algorithm
as state transition probabilities to inform the algorithm where syllable
boundary are likely to occur.

###  4.3 Decoding of the syllable boundaries

To decode the syllable boundaries, we construct a hidden markov model
characterized by the following:

  1. 1.

The hidden state space is a set of $N$ candidate onset positions
$S_{1},S_{2},\cdots,S_{N}$ discretized by the hop size, where $S_{N}$ is the
offset position of the last syllable.

  2. 2.

The state transition probability at the time instant $l$ associated with state
changes is defined by a priori duration distribution
$\mathcal{N}(d_{ij};\mu_{l},\sigma_{l}^{2})$, where $d_{ij}$ is the time
distance between states $S_{i}$ and $S_{j}$ ($j>i$). The length of the decoded
state sequence is equal to the total syllable number $L$ written in the score.

  3. 3.

The observation probability for the state $S_{j}$ is represented by its
corresponding value in the syllable ODF $p$, which is denoted as $p_{j}$.

The goal is to find the best onset position state sequence
$Q={q_{1}q_{2}\cdots q_{L}}$ for a given a priori duration sequence $M$, where
$q_{i}$ denotes the onset of the $i+1$th decoding syllable. $q_{0}$ and
$q_{L}$ are fixed as $S_{1}$ and $S_{N}$ as we expect that the onset of the
first syllable to be located at the beginning of the incoming audio and the
offset of the last syllable is located at the ending of the audio. One can
fulfill this assumption by truncating the silences at the beginning and at the
end of the incoming audio. According to the logarithmic form of Viterbi
algorithm [13], we define:

$\delta_{l}(i)=\max_{q_{1},q_{2},\cdots,q_{l}}{\log P[q_{1}q_{2}\cdots
q_{l},\,\mu_{1}\mu_{2}\cdots\mu_{l}]}$

with the initial step as follows:

$\displaystyle\delta_{1}(i)$
$\displaystyle=\log(\mathcal{N}(d_{1i};\mu_{1},\sigma_{1}^{2}))+\log(p_{i})$
$\displaystyle\psi_{1}(i)$ $\displaystyle=S_{1}$

with the following recursive step:

$\displaystyle\delta_{l}(j)$ $\displaystyle=\max_{1\leqslant
i<j}[\delta_{l-1}(i)+\log(\mathcal{N}(d_{ij};\mu_{l},\sigma_{l}^{2}))]+\log(p_{j})$
$\displaystyle\psi_{l}(j)$ $\displaystyle=\arg\max_{1\leqslant
i<j}[\delta_{l-1}(i)+\log(\mathcal{N}(d_{ij};\mu_{l},\sigma_{l}^{2}))]$

and the following termination step:

$\displaystyle\log P^{*}$ $\displaystyle=\max_{1\leqslant
i<N}[\delta_{L-1}(i)+\log(\mathcal{N}(d_{iN};\mu_{L},\sigma_{L}^{2}))]$
$\displaystyle q_{L}^{*}$ $\displaystyle=\arg\max_{1\leqslant
i<N}[\delta_{L-1}(i)+\log(\mathcal{N}(d_{iN};\mu_{L},\sigma_{L}^{2}))]$

Finally, the best offset position state sequence $Q$ is obtained by the
backtracking step. An example of the best boundary position decoding path is
showed in figure 6.

Figure 6: Illustration of an example of the best boundary position decoding
path (black) with $N=558$ and $L=8$. The intermediate states are omitted in
the vertical direction for a clear visualization. Blue curve: syllable ODF;
red horizontal lines: decoded boundary positions.

##  5 Experiments and results

###  5.1 Performance metrics

The syllable segmentation task consists on determining the time positions of
syllable boundaries. The proposed evaluation consists in comparing the
detected syllable onsets/offsets to their reference ones. We report F-measure
results in this paper. The definition of a correct segmented syllable is
borrowed from the note transcription literature [9]. For syllable onsets, we
choose an evaluation tolerance of $\pm\tau$ ms. For offsets, we choose an
evaluation tolerance of either (a) $\pm$20% of the syllable’s duration
annotation, or (b) $\pm\tau$ ms – whichever is larger. If both the onset and
the offset of a syllable lie within the tolerance of their annotated
counterparts and the syllable is correctly labeled, we consider that it’s
correctly segmented. We report the results for a tolerance of $\tau=0.05$
(seconds).

###  5.2 Results and discussion

Two state-of-the-art methods are set as baselines: Obin et al. [10] as a
traditional approach, and Schlüter et al. [16] as a deep learning method –
both already introduced. We explore Pons et al.’s CNNs design strategy [11,
12] as a way to improve our results. Code is available
online555https://github.com/ronggong/jingjuSyllabicSegmentaion/tree/v0.1.0.
All syllable ODFs are decoded by using the same approach (described in section
4.3). F-measure results are reported in Table 1. The evaluation can not be
performed by using peak-picking because the syllabic label can not be
attached. However, in section 2 we already discussed the potential problems of
peak-picking which are explicitly addressed with the proposed score-informed
method.

Table 1: Syllable segmentation results using different methods for estimating the syllable ODF. Numbers in parenthesis indicate the #filters in the second CNN layer. Syllable ODFs | #params | F-measure (%)
---|---|---
Late-fusion (32) | - | 86.37
Temporal (32) | 210,403 | 83.28
Timbral (32) | 185,420 | 84.83
Late-fusion (20) | - | 84.05
Temporal (20) | 132,127 | 83.86
Timbral (20) | 119,432 | 82.89
Schlüter et al. [16] | 535,290 | 81.93
Obin et al. [10] | - | 40.61
Figure 7: Three syllable ODFs. Red vertical line: decoded syllable onset time
positions. Black arrows: ground truth.

Results show that proposed architectures improve the state-of-the-art in all
cases. These results support the idea that using different filter shapes in
the first CNN layer is beneficial. Moreover, observe that #params is reduced
at least to the half – interpret #params (number of parameters of the CNN) as
a measure for the representational capacity of the network. This denotes how
efficient can be CNNs if these are designed to capture the relevant features
for the task at hand. Also note that if the #params of a network is reduced,
the over-fitting risk is reduced as well.

We also observe that late-fusing the predictions helps improving the results.
The interesting idea behind this approach is that individual networks still
need to be able to solve the task by their own means – and we propose fusing
models that are tailored towards learning complementary time-frequency scales:
temporal and timbral architectures.

We also see a big performance gap between Obin et al. and CNN-based methods.
For example in Figure 7 numerous noisy peaks can be observed for Obin et al.’s
method, what is significantly decreasing its performance. This result denotes
the difficulty of designing rules or handcrafting features for estimating the
syllable ODF.

###  5.3 Error analysis

We conduct an error analysis to study the causes of the segmentation errors
produced by our best performing method, what points out future research
directions for further improvement. We study falsely detected onsets in the
test dataset by comparatively observing the detected syllable onset positions,
the ground truth positions, the spectrogram and the corresponding scores.
Errors out of 0.05s tolerance are analyzed and 70 syllable onsets out of 512
are identified as false detections. Errors are then classified according to
their causes. Two main error classes are found: score-performance mismatching
and ambiguous syllable transitions.

The majority of the errors (71.42%) are caused by performance deviations from
the score: syllable insertion, syllable deletion or performing different
syllable durations than the score. The proposed syllable boundary decoding
algorithm is not able to preserve the correctness when the difference between
the score and its performance becomes too large. Furthermore, since the length
of the decoded state sequence of the algorithm must be equal to the number of
syllables in the score, singing syllable insertions and deletions not
expressed in the score can not be handled properly. One possible solution is
to incorporate more domain knowledge of jingju singing into the decoding
process – such as considering the probability of a padding-character to be
inserted/deleted, or the expectation of prolonging the last syllable in a
phrase.

The second source of errors include ambiguous syllable transitions (15.61%) –
such as transitions from vowel to vowel, from vowel or to semi-vowel, etc.
These errors are very difficult to be corrected because no prominent spectral
changes can be discerned within such transitions [15]. In future
investigations we shall detect “onset regions” rather than “onset time
positions” since this kind of syllable transitions usually manifest themselves
as gradual spectral changes.

##  6 Conclusions

This paper introduces a new score-informed method for the segmentation of
jingju singing into syllables. Two main contributions are presented in this
paper: (i) improvements for estimating the syllable ODF with CNNs, and (ii) a
method for incorporating score information into Viterbi’s algorithm for
estimating syllable boundaries.

The improvements to the CNN architecture consisted on using different filter
shapes in the first layer and late-fusing the predictions of two models –
designed to learn complementary representations (temporal and timbral). By
doing so, we increased the expressiveness of the first layer and enabled the
networks to efficiently capture different time-frequency scales useful for
detecting syllable onsets. Proposed models, with many filter shapes in the
first layer, have proven to be more effective than the state-of-the-art model
based in a single filter shape in the first layer.

Moreover, we proposed an a priori duration model that describes the
probability of a syllable boundary given the score. The likelihood of a
syllable boundary is shaped with the a priori duration model and incorporated
into Viterbi’s algorithm as state transition probabilities – this being the
core of the proposed score-informed Viterbi algorithm.

We validated the proposed method on a jingju a cappella singing dataset, which
achieved better performance than the state-of-the-art. Although the proposed
ameliorations helped improving our results, the proposed method has not yet
solved the task. To this end, we plan to investigate incorporating more domain
knowledge into the decoding process, and to further improve the ODF with RNNs
– that has proven to be very useful for similar tasks [1, 4].

##  7 Acknowledgments

We are grateful for the GPUs donated by NVidia. This work is partially
supported by the Maria de Maeztu Programme (MDM-2015-0502) and the European
Research Council under the European Union’s Seventh Framework Program, as part
of the CompMusic project (ERC grant agreement 267583).

## References

  * [1] Sebastian Böck, Florian Krebs, and Gerhard Widmer.  Joint beat and downbeat tracking with recurrent neural networks.  In ISMIR, New York City, USA, 2016.
  * [2] Paul Boersma.  Praat, a system for doing phonetics by computer.  Glot International, 5(9/10):341–345, 2001.
  * [3] Djork-Arné Clevert, Thomas Unterthiner, and Sepp Hochreiter.  Fast and accurate deep network learning by exponential linear units (elus).  CoRR, abs/1511.07289, 2015.
  * [4] Florian Eyben, Sebastian Böck, Björn W Schuller, and Alex Graves.  Universal onset detection with bidirectional long short-term memory neural networks.  In ISMIR, Utrecht, Netherlands, 2010.
  * [5] Brett Kessler and Rebecca Treiman.  Syllable structure and the distribution of phonemes in english syllables.  Journal of Memory and language, 37(3):295–311, 1997.
  * [6] D. Kingma and J. Ba.  Adam: A method for stochastic optimization.  arXiv:1412.6980, 2014.
  * [7] A. Klapuri.  Sound onset detection by applying psychoacoustic knowledge.  In ICASSP, Phoenix, USA, 1999.
  * [8] Anna M. Kruspe.  Keyword spotting in a-capella singing.  In ISMIR, Taipei, Taiwan, 2014.
  * [9] Emilio Molina, Ana M. Barbancho, Lorenzo J. Tardón, and Isabel Barbancho.  Evaluation Framework for Automatic Singing Transcription.  In ISMIR, Taipei, Taiwan, 2014.
  * [10] N. Obin, F. Lamare, and A. Roebel.  Syll-O-Matic: An adaptive time-frequency representation for the automatic segmentation of speech into syllables.  In ICASSP, Vancouver, Canada, 2013.
  * [11] Jordi Pons and Xavier Serra.  Designing efficient architectures for modeling temporal features with convolutional neural networks.  In ICASSP, New orleans, USA, 2017.
  * [12] Jordi Pons, Olga Slizovskaia, Rong Gong, Emilia Gómez, and Xavier Serra.  Timbre analysis of music audio signals with convolutional neural networks.  arxiv:1703.06697, 2017.
  * [13] L. Rabiner.  A tutorial on hidden Markov models and selected applications in speech recognition.  Proceedings of the IEEE, 77(2):257–286, February 1989.
  * [14] Rafael Caro Repetto and Xavier Serra.  Creating a Corpus of Jingju (Beijing Opera) Music and Possibilities for Melodic Analysis.  In ISMIR, Taipei, Taiwan, 2014.
  * [15] Odette Scharenborg, Vincent Wan, and Mirjam Ernestus.  Unsupervised speech segmentation: An analysis of the hypothesized phone boundaries.  The Journal of the Acoustical Society of America, 127(2):1084–1095, 2010.
  * [16] J. Schlüter and S. Böck.  Improved musical onset detection with convolutional neural networks.  In ICASSP, Florence, Italy, 2014.
  * [17] Chee-Chuan Toh, Bingjun Zhang, and Ye Wang.  Multiple-feature fusion based onset detection for solo singing voice.  In ISMIR, Philadelphia, USA, 2008.
  * [18] Elizabeth Wichmann.  Listening to Theatre: The Aural Dimension of Beijing Opera.  University of Hawaii Press, 1991.