# Self-organized PT-symmetry of exciton-polariton condensate in a double-well
potential

P.A. Kalozoumis  Materials Science Department, School of Natural Sciences,
University of Patras, GR-26504 Patras, Greece  Hellenic American University,
436 Amherst st, Nashua, NH 0306 USA  Institute of Electronic Structure and
Laser, FORTH, GR-70013 Heraklion, Crete, Greece  D. Petrosyan  Institute of
Electronic Structure and Laser, FORTH, GR-70013 Heraklion, Crete, Greece  A.
Alikhanyan National Science Laboratory (YerPhI), 0036 Yerevan, Armenia

###### Abstract

We investigate the dynamics and stationary states of a semiconductor exciton-
polariton condensate in a double well potential. We find that upon the
population build up of the polaritons by above-threshold laser pumping,
coherence relaxation due to the phase fluctuations of the polaritons drives
the system into a stable fixed point corresponding to a self-organized PT-
symmetric phase.

##  I Introduction

One of the prominent research directions in semiconductor optics is the study
of exciton-polariton condensation in microcavities. Exciton-polaritons are
hybrid quasi-particles of strongly coupled quantum well excitons and cavity
photons [1, 2] which retain the properties of both matter and light. The
excitonic part mediates effective interactions between the polaritons, giving
rise to interesting nonlinear properties, whereas the small effective mass of
the photonic component enables Bose-Einstein condensation even at ambient
temperatures [3, 4, 5], in contrast to their ultracold atomic counterparts
[6]. The short lifetime of the polariton condensate renders it an open system
that requires continuous replenishing from the excitonic reservoir via
external pumping. After their experimental realization [7, 8], the polariton
condensates have been shown to be ideal system for studies of many effects at
the interface of non-equilibrium physics and nonlinear dynamics.

The intrinsic nonlinear dynamics of polariton systems lead to a variety of
effects, such as the appearance of a Mach-Cherenkov cone in a supersonic flow
[9], the formation of quantized vortices [10], and dark solitons [11].
Moreover, the polariton condensates can be engineered with high precision by
the external laser fields [8, 12, 13, 14]. Finally, such systems are promising
candidates for various applications in photonic devices, such as switches,
gates and transistors [3], as well as for quantum simulators of interacting
spin models [15].

The “open” nature of the system, featuring gain and loss, leads to interesting
implications when the dissipative dynamics become pseudo-Hermitian. This is
the case in parity-time (PT) symmetric setups, where dissipation losses are
exactly balanced by the pumping gain. Systems with PT-symmetry has been a
flourishing and broad research field, extending from quantum mechanics [16]
and field theory [17] to optics [18] and acoustics [19].

The interplay between the inherent losses and the laser pumping in such a way
as to preserve the PT symmetry of the system provides an effective framework
where a polariton system can exhibit coherent, Hermitian-like dynamics for
relatively long times. Recent works have shown promising results, such as
permanent Rabi oscillations [20], multistability and condensation below
threshold [21], exceptional points in polaritonic cavities below lasing
threshold [22], and coherent oscillations of a two species polariton mixture
in a double well [23]. The latter has been shown to be able to simulate the
dynamics of a pair of spin-1/2 particles (qubits) in the presence of exchange
interaction. Yet, polariton structures in the framework of PT symmetry have
not been extensively studied yet, and more efforts are required to understand
the rich landscape of phenomena which emerge from this framework.

In this work we study the dynamics of an exciton-polariton condensate in a
double well potential, in the presence of time-varying exciton populations and
phase fluctuations. We consider the coupled-mode equations for the polaritons
supplemented by the rate equations for the laser-pumped exciton reservoirs,
and derive analytically the steady state solutions for the exciton and
polariton populations as well as their coherence. We find that, when the total
pumping rate is above threshold, the system automatically attains the PT
symmetric state, independently of the pumping rates of the individual sites.
Employing numerical simulations for several different pumping rates and
initial conditions, we verify our analytical findings. We also study the
stability and robustness of our results in the presence of phase noise caused
by the unavoidable phase fluctuations of the polaritons.

##  II The exciton-polariton system

Figure 1: Schematic top view (left panel) and side view (right panel) of a
polariton system in a double quantum well. Spatially shaped pumping lasers
populate with rates $P_{L}$ and $P_{R}$ the reservoir excitons $n_{L,R}$,
which decay via recombination with rates $\Gamma$ and energy-relax and scatter
into the polariton condensate with rate $R$. The pumping lasers also create
the confining potentials for the polaritons $\psi_{L}$ and $\psi_{R}$, which
decay with rates $\kappa$, are continuously replenished by reservoir excitons
with rates $Rn_{L,R}$, while interacting with each other via the Josephson
coupling $J$.

The system under consideration is schematically illustrated in Fig. 1. One or
more layers of semiconductor quantum wells are placed inside the semiconductor
microcavity near the antinode of the resonant cavity field mode. Spatially
shaped pumping lasers replenish continuously the exciton reservoirs and
simultaneously create confining potentials for the polariton condensate.
Assuming a tight-binding double-well potential, the exciton-polariton system
can be described by the following set of equations for the polariton
condensate wavefunctions $\psi_{L}$ and $\psi_{R}$ in the left ($L$) and right
($R$) wells [24]:

$\displaystyle i\partial_{t}\psi_{L}$ $\displaystyle=$
$\displaystyle\left[\epsilon_{L}+\eta|\psi_{L}|^{2}\right]\psi_{L}+\frac{i}{2}\left[Rn_{L}-\kappa\right]\psi_{L}-J\;\psi_{R},\qquad$
(1a) $\displaystyle i\partial_{t}\psi_{R}$ $\displaystyle=$
$\displaystyle\left[\epsilon_{R}+\eta|\psi_{R}|^{2}\right]\psi_{R}+\frac{i}{2}\left[Rn_{R}-\kappa\right]\psi_{R}-J^{*}\psi_{L},$
(1b)

where $\epsilon_{L,R}$ are the single-particle energies, $\eta$ is the
nonlinear interaction strength, $\kappa$ is the decay rate of the polaritons
due to the exciton recombination and cavity photon losses (assumed the same
for both wells), and $J$ is the Josephson (tunnel) coupling between the wells.
The polariton equations are supplemented by the equations for the populations
$n_{L,R}$ of the reservoir excitons,

$\displaystyle\partial_{t}n_{L}$ $\displaystyle=$ $\displaystyle P_{L}-\Gamma
n_{L}-Rn_{L}|\psi_{L}|^{2},$ (2a) $\displaystyle\partial_{t}n_{R}$
$\displaystyle=$ $\displaystyle P_{R}-\Gamma n_{R}-Rn_{R}|\psi_{R}|^{2},$ (2b)

which are created by laser pumping with rates $P_{L,R}$, decay with rate
$\Gamma$, and scatter into the polariton condensate with rate $R$.

In Appendix A we briefly outline the PT-symmetry conditions for a condensate
in a double well potential. Neglecting for the moment the non-linearity
$\eta$, the PT-symmetry condition is satisfied when
$\epsilon_{L,R}=\epsilon\,(=0)$ and the gain in one well exactly compensates
the losses in the other,

$\gamma_{L}\equiv\frac{1}{2}[Rn_{L}-\kappa]=-\frac{1}{2}[Rn_{R}-\kappa]\equiv-\gamma_{R},$
(3)

as per Eqs. (1a) and (1b), which leads to $n_{L}+n_{R}=\frac{2\kappa}{R}$. The
threshold pumping at which the polariton condensate starts to form can be
obtained from the condition that the sum of the gain and loss in both wells is
non-negative, $\gamma_{L}+\gamma_{R}\geq 0$. With Eq. (3), this condition is
equivalent to

$n_{L}+n_{R}\geq\frac{2\kappa}{R}.$ (4)

Note that exactly at the threshold, this is the same condition as for the PT-
symmetry. If we consider the stationary regime for the reservoir excitons,
$\partial_{t}n_{L,R}=0$, we find from Eqs. (2) the steady-state values

$n_{L,R}=\frac{P_{L,R}}{\Gamma+R|\psi_{L,R}|^{2}}.$ (5)

Exactly at the threshold for condensate formation, the values of the polariton
populations in both wells, $|\psi_{L,R}|^{2}$, are marginally equal to zero
and we have $n_{L,R}\simeq P_{L,R}/\Gamma$. Substituting these values into Eq.
(4), we find the threshold pumping condition

$P_{L}+P_{R}\geq\frac{2\kappa\Gamma}{R},$ (6)

above which the condensate begins to form, while for
$P_{L}+P_{R}<2\kappa\Gamma/R$ the condensate decays to zero.

Figure 2: Dynamics of the polariton populations $|\psi_{L,R}|^{2}$ (upper
panels), the reservoir excitons $n_{L,R}$ (insets), and the coherence $\Theta$
(lower panels), as obtained from the numerical solution of Eqs. (1-2) for the
parameters $\epsilon=0$, $\kappa=10J$, $\Gamma=2J$, $R=0.02J$ (time is in
units of $J^{-1}$), corresponding to the threshold values of pumping
$P_{L}+P_{R}=2\kappa\Gamma/R=2000J$ and the steady-state exciton populations
$n_{L}+n_{R}=2\kappa/R=1000$ (green dashed lines in the insets). The initial
exciton populations are always taken as $n_{L,R}=P_{L,R}/\Gamma$, while the
initial polariton amplitudes $\psi_{L,R}$ are seeded with small (complex)
values. (a) Linear case ($\eta=0$) with the pumping rates $P_{L}=1000J$ and
$P_{R}=990J$ slightly below threshold, leading to decay of the initial
polariton populations $|\psi_{L}(t)|^{2}$ and $|\psi_{R}(t)|$ and steady-state
exciton population $n_{L}+n_{R}<2\kappa/R$. (b) Same as in (a), but for
stronger pumping rates $P_{L}=1080$ and $P_{R}=1020$ above threshold and
initial conditions $\textrm{Re}\Theta(0)=0$, leading to the initial build up
of the polariton populations and their continuous Rabi-like oscillations,
while $\textrm{Re}\Theta(t)=0\;\forall\;t$, and the exciton population
$n_{L}+n_{R}$ oscillating slightly above the threshold value $2\kappa/R=1000$.
(c) Same as in (b) but for the initial conditions $\textrm{Re}\Theta(0)\neq
0$, leading to a steady-state of the system. (d) Same as in (b) with the
initial conditions $\textrm{Re}\Theta(0)=0$ but in the presence of
nonlinearity $\eta=0.3J$ that couples $\textrm{Re}\Theta$ and
$\textrm{Im}\Theta$, leading to a steady-state of the system.

In the upper panels of Fig. 2 we show the polariton populations
$|\psi_{L,R}|^{2}$ for different pumping rates and initial conditions, with
and without non-linear interaction, as obtained from the numerical solution of
Eqs. (1-2). The insets show the evolution of the exciton populations $n_{L}$
and $n_{R}$ and their sum $n_{L}+n_{R}$. For pumping below threshold, we
observe a decay of the initial (seed) polariton populations with rate
$\gamma_{L}+\gamma_{R}\simeq\frac{R}{2\Gamma}(P_{L}+P_{R})-\kappa<0$,
accompanied by Rabi-like oscillations, while the exciton populations settle to
$n_{L,R}=P_{L,R}/\Gamma$. For pumping above threshold, the polariton
populations grow until reaching certain values $|\psi_{L,R}|^{2}$ at which
$n_{L}+n_{R}\simeq\frac{2\kappa}{R}$, while the Rabi-like oscillations persist
or are eventually damped, depending on the initial conditions or presence of
non-linear interaction, as discussed below. Remarkably, the polariton and
exciton populations increase and decrease, respectively, reaching the same
stationary values which satisfy the PT-symmetry conditions, independently of
the pumping rates, as long as pumping is retained above threshold.

##  III Equivalence of the PT-symmetry and steady state conditions

To understand the dynamics of the system, it is convenient to express Eqs. (1)
in terms of the polariton populations $|\psi_{L,R}|^{2}$ and coherence
$\Theta\equiv\psi_{L}\psi_{R}^{*}$ as

$\displaystyle\partial_{t}|\psi_{L}|^{2}$ $\displaystyle=$ $\displaystyle
2\gamma_{L}|\psi_{L}|^{2}+2J\textrm{Im}\Theta,$ (7a)
$\displaystyle\partial_{t}|\psi_{R}|^{2}$ $\displaystyle=$ $\displaystyle
2\gamma_{R}|\psi_{R}|^{2}-2J\textrm{Im}\Theta,$ (7b)
$\displaystyle\partial_{t}\Theta$ $\displaystyle=$
$\displaystyle-i\left[\epsilon_{L}-\epsilon_{R}+\eta(|\psi_{L}|^{2}-|\psi_{R}|^{2})\right]\Theta$
(7c) $\displaystyle+(\gamma_{L}+\gamma_{R})\Theta-
iJ\left(|\psi_{L}|^{2}-|\psi_{R}|^{2}\right).$

Note that below threshold, $(\gamma_{L}+\gamma_{R})<0$, both the polariton
populations and their coherence decay to zero, as already mentioned above.

Let us assume $\epsilon_{L,R}=0$ and consider first the case of vanishing
nonlinearity $\eta=0$. Equation (7c) indicates that the coherence decays only
if its real part is nonzero. In turn, the solution for the real part of the
coherence is

$\textrm{Re}\Theta(t)=\textrm{Re}\Theta(0)\;e^{\int_{0}^{t}(\gamma_{L}+\gamma_{R})dt^{\prime}}.$
(8)

Hence, if initially $\textrm{Re}\Theta(0)=0$, it will remain so at later
times, $\textrm{Re}\Theta(t)=0\;\forall\;t>0$. Then the dynamics of the
system, if pumped above threshold, will exhibit continuous Rabi-like
oscillations with frequency $J$, while no steady state will be attained, as in
Fig. 2(b).

In practice, however, even if initially we have $\textrm{Re}\Theta(0)=0$
[e.g., either $\psi_{L}(0)=0$ or $\psi_{R}(0)=0$], the unavoidable phase
fluctuations of the polaritons will eventually lead to the appearance of
finite $\textrm{Re}\Theta(t)\neq 0$, which in turn will result in the decay of
coherence and drive the system to the steady state. Equivalently, if we have
initially $\textrm{Re}\Theta(0)\neq 0$, the system can still exhibit initially
Rabi-like oscillations, but then it will eventually attain the steady state,
as in Fig. 2(c). Finally, as seen from Eq. (7c) the nonlinear interaction
couples the real and imaginary parts of the coherence $\Theta$ with the rate
$\eta(|\psi_{L}|^{2}-|\psi_{R}|^{2})$. Hence, in the presence of nonlinearity
$\eta\neq 0$, we expect the eventual decay of the coherence with the system
attaining the steady state, for any initial conditions and independent on the
phase fluctuations, as in Fig. 2(d).

Setting the time derivative in the left-hand side of the Eq. (7c) equal to
zero, we find the steady state is reached when

$R[n_{L}+n_{R}]-2\kappa=0,\quad\mathrm{and}\quad|\psi_{L}|^{2}=|\psi_{R}|^{2}.$
(9)

Remarkably, the first equation corresponds exactly to the PT-symmetry
condition $n_{L}+n_{R}=\frac{2\kappa}{R}$ discussed above. Moreover, this
condition is satisfied even in the presence of nonlinear interaction $\eta\neq
0$, because the equal polariton populations as per the second equation lead to
exactly the same energy shifts $\eta|\psi_{L,R}|^{2}$ of the polaritons in
both wells. In other words, for any initial conditions, and provided the total
pumping is above threshold as per Eq. (6) but otherwise arbitrary $P_{L}$ and
$P_{R}$, the system attains a stable fixed point corresponding to the PT-
symmetric state. Even when no steady state exists or is yet reached, the PT
condition in Eq. (9) is approximately satisfied, as seen in the insets of Fig.
2.

Figure 3: Dynamics of the polariton populations $|\psi_{L,R}|^{2}$ (upper
panels), the coherence $\Theta$ (middle panels), and the left-well polariton
field correlation function $g^{(1)}(t)$ (lower panels), for the same
parameters as in Fig. 2 with the addition of phase fluctuations causing
decoherence with rate $\xi=0.05J$, as obtained from the the ensemble averaged
solution of Eqs. (1-2). (a) Linear case ($\eta=0$) with the below threshold
pumping $P_{L}=1000J$ and $P_{R}=990J$ leading to the exponential decay of the
correlation function $g^{(1)}(t)\propto
e^{-(\xi-|\gamma_{L}+\gamma_{R}|/2)t}$. (b) Same as in (a), but for stronger
pumping above threshold $P_{L}=1080$ and $P_{R}=1020$. Now, independently of
the initial value of $\textrm{Re}\Theta(0)$, the phase fluctuations cause
exponential decay of $g^{(1)}(t)\propto e^{-\xi t}$ while the system
approaches the steady state. (c) Same as in (b) but in the presence of
nonlinearity $\eta=0.3J$, causing accelerated decay of $g^{(1)}(t)\propto
e^{-0.057t}$ and faster approach of the system to the steady state. In the
lower panels, we also show the correlation functions $g^{(1)}(t)$ obtained
from the long-time average of the system dynamics (the oscillating tail of
$g^{(1)}$ is due to the finite length of the time series).

Combining Eqs. (5) and (9), we find that steady state polariton populations
are

$|\psi_{L}|^{2}=|\psi_{R}|^{2}=\frac{P_{L}+P_{R}}{2\kappa}-\frac{\Gamma}{R},$
(10)

while the exciton populations are

$n_{L,R}=\frac{2\kappa P_{L,R}}{R(P_{L}+P_{R})}.$ (11)

Using these stationary values for $n_{L,R}$ and $|\psi_{L,R}|^{2}$ in Eqs.
(7a) and (7b) in the steady state we obtain

$\textrm{Im}(\Theta)=\frac{\kappa\Gamma(P_{L}-P_{R})}{2RJ(P_{L}+P_{R})}-\frac{P_{L}-P_{R}}{4J},$
(12)

and from
$|\psi_{L}\psi_{R}^{*}|^{2}=[\textrm{Re}(\psi_{L}\psi_{R}^{*})]^{2}+[\textrm{Im}(\psi_{L}\psi_{R}^{*})]^{2}$
we obtain

$\textrm{Re}(\Theta)=\mathcal{D}\sqrt{4J^{2}(P_{L}+P_{R})^{2}-\kappa^{2}(P_{L}-P_{R})^{2}}$
(13)

where

$\mathcal{D}=\frac{P_{L}+P_{R}-2\kappa\Gamma/R}{4J\gamma(P_{L}+P_{R})}.$

These results are verified by the numerical simulations illustrated in Fig. 2
and they equally hold for any value the nonlinearity strength $\eta$.

##  IV Phase fluctuations

As mentioned above, the coherence of the polariton condensate will decay due
to the phase fluctuations that are always present in realistic quantum
systems. We therefore incorporate the phase noise in our numerical
calculations and investigate how it modifies the dynamics of the polaritons
and the coherence. We model the phase fluctuations as the standard Wiener
process for stochastic differential equations. Thus, the single-particle
energies $\epsilon_{L,R}$ in Eqs. (1) become Gaussian stochastic variables
with the mean $\braket{\epsilon_{L,R}}=0$ and variance $\sigma^{2}=2\xi/\delta
t$, where $\xi$ is the decoherence rate and $\delta t$ is the time step for
picking a new random energy.

In Fig. 3 we show the results of our numerical simulations as obtained upon
the ensemble average over $N=1000$ independent realizations of the system
dynamics. We compute the first-order correlation functions $g^{(1)}(t)$, which
quantify the coherence for the polaritonic fields, via

$g^{(1)}(t)=\frac{\langle\psi(t_{0})\psi(t)\rangle}{\sqrt{\langle|\psi(t_{0})|^{2}\rangle\langle|\psi(t)|^{2}\rangle}},$
(14)

where $\psi=\psi_{L}$ or $\psi_{R}$, and $\braket{\ldots}$ denotes the
ensemble average.

Below the pumping threshold, the polariton fields decay with rate
$\gamma_{L}+\gamma_{R}<0$, but the phase fluctuations with rate $\xi$ causes
even faster decay of coherence, $g^{(1)}(t)\propto
e^{-(\xi-|\gamma_{L}+\gamma_{R}|/2)t}$, as seen in Fig. 3(a). For pumping
above threshold, the phase noise causes exponential decay of the correlation
function $g^{(1)}(t)\propto e^{-\xi t}$, while the system approaches the
steady state independently of the initial value of coherence
$\textrm{Re}\Theta(0)$, as seen in Fig. 3(b). Including also the nonlinear
interaction $\eta\neq 0$ further accelerates the decay of the correlation
function and the system approaches the steady state even faster.

We finally note that for an ergodic system the ensemble-averaged and time-
averaged correlation functions are equivalent. To verify whether our polariton
system is ergodic, we also compute the field correlation function

$g^{(1)}(t)=\frac{\int_{t_{\textrm{i}}}^{t_{\textrm{f}}}d\tau\psi(\tau)\psi^{*}(\tau+t)}{\sqrt{\int_{t_{\textrm{i}}}^{t_{\textrm{f}}}d\tau|\psi(\tau)|^{2}\int_{t_{\textrm{i}}}^{t_{\textrm{f}}}d\tau|\psi(\tau+t)|^{2}}}$
(15)

resulting from a single, long-time trajectory with
$t_{\textrm{f}}-t_{\textrm{i}}=3000/J$. As seen in Fig. 3 (lower panels), the
computed ensemble-averaged and time-averaged correlation functions coincide to
a very approximation, attesting to the ergodicity of our system.

##  V Conclusions

To summarize, we have studied an exciton-polariton system in a double-well
potential, taking into account the dynamics of the reservoir excitons and the
polaritons. We have found that for pumping of the excitons above the total
threshold value for the formation of the polariton condensate, the exciton
populations attain the values that satisfy the PT-symmetry condition for the
polariton condensate, independent of the pumping rates of the individual
wells. Employing the population-coherence equations, we interpreted the
corresponding dynamics and revealed the stable fixed point, or the steady
state, that the system approaches. To make our analysis experimentally
relevant, we have taken also into account the phase fluctuations present in
any realistic system, and computed the first-order correlation functions for
the polariton fields, which revealed the coherence decay with the
corresponding rate.

We note that our results apply to moderate non-linear interaction strength and
small differences in pumping rates of the two wells. For large difference in
the pumping rates, the strong non-linear energy shift of the polariton
condensate energy may lead to self-trapping and break-up of the PT symmetry
[25]

##  VI Acknowledgments

We thank P.G. Savvidis, H. Ohadi, and A.F. Tzortzakakis for fruitful
discussions. This work was co-financed by Greece (General Secretariat for
Research and Technology), and the European Union (European Regional
Development Fund), in the framework of the bilateral Greek-Russian Science and
Technology collaboration on Quantum Technologies (POLISIMULATOR project.)

##  Appendix A Polariton condensate in a PT-symmetric double well

Consider a polariton condensate in a double well potential described by the
coupled-mode equations

$\displaystyle i\partial_{t}\psi_{L}$ $\displaystyle=$
$\displaystyle(\epsilon_{L}+i\gamma_{L})\psi_{L}+\eta|\psi_{L}|^{2}\psi_{L}-J\;\psi_{R},$
(16a) $\displaystyle i\partial_{t}\psi_{R}$ $\displaystyle=$
$\displaystyle(\epsilon_{R}+i\gamma_{R})\psi_{R}+\eta|\psi_{R}|^{2}\psi_{R}-J^{*}\\!\psi_{L},$
(16b)

where $\epsilon_{L,R}$ are the single-particle energies, $\gamma_{L,R}$ are
the incoherent loss ($\gamma<0$) or gain ($\gamma>0$) rates at each well,
$\eta$ is the nonlinear interaction strength, and $J$ is the Josephson
coupling between the wells.

Figure 4: Imaginary part of the eigenvalues $\lambda_{\pm}$ in Eq. (18). For
$\gamma<J(=1)$ we have a PT symmetric phase with real eigenvalues. Above the
bifurcation point at $\gamma=J$, the system enters the PT broken phase with
imaginary eigenvalues.

If we set $\epsilon_{L,R}=0$, assume negligibly weak nonlinearity,
$\eta|\psi|^{2}\ll J$, and set $\gamma_{L}=-\gamma_{R}=\gamma$ so that the
loss at the right well is exactly compensated by the gain at the left well, we
obtain a PT-symmetric Hamiltonian matrix corresponding to Eqs. (16) [23]:

$\mathcal{H}=\begin{pmatrix}i\gamma&-J\\\ -J^{*}&-i\gamma\end{pmatrix}.$ (17)

Its eigenvalues and the corresponding eigenvectors are given by

$\lambda_{\pm}=\pm\sqrt{|J|^{2}-\gamma^{2}}$ (18)

and

$|\pm\rangle=\left[\left(\sqrt{|J|^{2}-\gamma^{2}}\pm
i\gamma\right)|L\rangle\mp J^{*}|R\rangle\right]/N_{\pm}$ (19)

with $N_{\pm}$ the normalization factors. For $\gamma<|J|$, the eigenvalue
spectrum is real and the dynamics is Hermitian-like. For $\gamma>|J|$, the
eigenvalues become imaginary and the system enters the PT-broken phase. The
case $|J|=\gamma$ corresponds to the exceptional point of the system where the
eigenvalues become degenerate and the eigenstates coalesce. Figure 4
illustrates the dependence of imaginary part of the eigenvalues on the
loss/gain parameter $\gamma$.

## References

  * [1] H. Deng, H. Haug, and Y. Yamamoto, Rev. Mod. Phys. 82, 1489 (2010).
  * [2] I. Carusotto and C. Ciuti, Rev. Mod. Phys. 85, 299 (2013).
  * [3] K. G. Lagoudakis, M. Wouters, M. Richard, A. Baas, I. Carusotto, R. Andre, Le Si Dang, and B. Deveaud-Pledran, Nature Phys. 4, 706 (2008).
  * [4] M. Yamaguchi, K. Kamide, R. Nii, T. Ogawa, and Y. Yamamoto, Phys. Rev. Lett. 111, 026404 (2013).
  * [5] C. Schneider, A. Rahimi-Iman, N. Y. Kim, J. Fischer, I. G. Savenko, M. Amthor, M. Lermer, A. Wolf, L. Worschech, V. D. Kulakovskii, I. A. Shelykh, M. Kamp, S. Reitzenstein, A. Forchel, Y. Yamamoto, and S. Hofling, Nature 497, 348 (2013).
  * [6] V. Goblot, H.S. Nguyen, I. Carusotto, E. Galopin, A. Lematre, I. Sagnes, A. Amo, and J. Bloch, Phys. Rev. Lett. 117, 217401 (2016).
  * [7] J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. Keeling, F. M. Marchetti, M. H. Szymanska, R. André, J. L. Staehli, V. Savona, P. B. Littlewood, B. Deveaud, and L. S. Dang, Nature 443, 409 (2006).
  * [8] R. Balili, V. Hartwell, D. Snoke, L. Pfeiffer, K. West, Science 316, 1007 (2007).
  * [9] A. Amo, J. Lefrére, S. Pigeon, C. Adrados, C. Ciuti, I. Carusotto, R. Houdre, E. Giacobino, and A. Bramati, Nat. Phys. 5, 805 (2009).
  * [10] L. Dominici, R. Carretero-Gonzalez, J. Cuevas-Maraver, A. Gianfrate, A. S. Rodrigues, D. J. Frantzeskakis, P. G. Kevrekidis, G. Lerario, D. Ballarini, M. De Giorgi, G. Gigli, and D. Sanvitto, Nature Comm. 9, 1467 (2018).
  * [11] R. Carretero-Gonzalez, J. Cuevas-Maraver, D. J. Frantzeskakis, T. P. Horikis, P. G. Kevrekidis, and A. S. Rodrigues, Phys. Lett. A 381, 3805 (2017).
  * [12] H. Ohadi, A. J. Ramsay, H. Sigurdsson, Y. del Valle-Inclan Redondo, S. I. Tsintzos, Z. Hatzopoulos, T. C. H. Liew, I. A. Shelykh, Y. G. Rubo, P. G. Savvidis, and J. J. Baumberg, Phys. Rev. Lett. 119, 067401 (2017).
  * [13] H. Ohadi, Y. del Valle-Inclan Redondo, A. J. Ramsay, Z. Hatzopoulos, T. C. H. Liew, P. R. Eastham, P. G. Savvidis, and J. J. Baumberg, Phys. Rev. B 97, 195109 (2018).
  * [14] K. Orfanakis, A.F. Tzortzakakis, D. Petrosyan, P.G. Savvidis, and H. Ohadi, Phys. Rev. B (2021).
  * [15] N. G. Berloff, M. Silva, K. Kalinin, A. Askitopoulos, J. D. Töpfer, P. Cilibrizzi,W. Langbein, and P. G. Lagoudakis, Nature Mat. 16, 1120 (2017).
  * [16] C. M. Bender, D. C. Brody, and H. F. Jones, Phys. Rev. Lett. 89, 270401 (2002).
  * [17] C. M. Bender, D. C. Brody, and H. F. Jones, Phys. Rev. D 70, 025001 (2004).
  * [18] R. El-Ganainy, K. G. Makris, D. N. Christodoulides and Z. H. Musslimani, Opt. Lett. 32, 2632 (2007).
  * [19] R. Fleury, D. Sounas, and A. Alú, Nat. Comm. 6, 5905 (2015).
  * [20] I. Y. Chestnov, S. S. Demirchyan, A. P. Alodjants, Y. G. Rubo, and A. V. Kavokin, Sci. Rep. 6, 19551 (2016).
  * [21] J.-Y. Lien, Y.-N. Chen, N. Ishida, H.-B. Chen, C.-C. Hwang, and F. Nori Phys. Rev. B 91, 024511 (2015).
  * [22] J. B. Khurgin, Optica 7(8), 1015 (2020).
  * [23] P. A. Kalozoumis, G. M. Nikolopoulos, and D. Petrosyan, EPL 129, 37003 (2020).
  * [24] M. Wouters and I. Carusotto, Phys. Rev. Lett. 99, 140402 (2007).
  * [25] A. A. Sukhorukov, Z. Xu, and Y.S. Kivshar, Phys. Rev. A 82, 043818 (2010).


* Hinton et al. (2009) Hinton, J. A., Skilton, J. L., Funk, S., et al. 2009, ApJ, 690, L101, doi: 10.1088/0004-637X/690/2/L101
  * Johnston et al. (1992) Johnston, S., Manchester, R. N., Lyne, A. G., et al. 1992, ApJ, 387, L37, doi: 10.1086/186300
  * Kambe et al. (2013) Kambe, E., Yoshida, M., Izumiura, H., et al. 2013, PASJ, 65, 15, doi: 10.1093/pasj/65.1.15
  * Khangulyan et al. (2008) Khangulyan, D., Aharonian, F., & Bosch-Ramon, V. 2008, MNRAS, 383, 467, doi: 10.1111/j.1365-2966.2007.12572.x
  * Khangulyan et al. (2011) Khangulyan, D., Aharonian, F. A., Bogovalov, S. V., & Ribó, M. 2011, ApJ, 742, 98, doi: 10.1088/0004-637X/742/2/98
  * Khangulyan et al. (2007) Khangulyan, D., Hnatic, S., Aharonian, F., & Bogovalov, S. 2007, MNRAS, 380, 320, doi: 10.1111/j.1365-2966.2007.12075.x
  * Kirk et al. (1999) Kirk, J. G., Ball, L., & Skjæraasen, O. 1999, Astroparticle Physics, 10, 31, doi: 10.1016/S0927-6505(98)00041-3
  * Krause et al. (2017) Krause, M., Pueschel, E., & Maier, G. 2017, Astroparticle Physics, 89, 1
  * Kravtsov et al. (2020) Kravtsov, V., Berdyugin, A. V., Piirola, V., et al. 2020, A&A, 643, A170, doi: 10.1051/0004-6361/202038745
  * Li et al. (2017) Li, J., Torres, D. F., Cheng, K.-S., et al. 2017, ApJ, 846, 169, doi: 10.3847/1538-4357/aa7ff7
  * Li et al. (2011) Li, J., Torres, D. F., Zhang, S., et al. 2011, ApJ, 744, L13, doi: 10.1088/2041-8205/744/1/l13
  * Li & Ma (1983) Li, T.-P., & Ma, Y.-Q. 1983, ApJ, 272, 317, doi: 10.1086/161295
  * Lomb (1976) Lomb, N. R. 1976, Ap&SS, 39, 447, doi: 10.1007/BF00648343
  * Maier & Holder (2017) Maier, G., & Holder, J. 2017, International Cosmic Ray Conference, 35, 747\.  https://arxiv.org/abs/1708.04048
  * Maier & VERITAS Collaboration (2017) Maier, G., & VERITAS Collaboration. 2017, International Cosmic Ray Conference, 35, 729.  https://arxiv.org/abs/1708.04045
  * Malyshev et al. (2019) Malyshev, D., Chernyakova, M., Santangelo, A., & Pühlhofer, G. 2019, Astronomische Nachrichten, 340, 465, doi: 10.1002/asna.201913605
  * Manset & Donati (2003) Manset, N., & Donati, J.-F. 2003, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4843, ESPaDOnS; an exhelle spectro-polarimetric device for the observation of stars, ed. S. Fineschi, 425–436, doi: 10.1117/12.458230
  * Martí-Devesa & Reimer (2020) Martí-Devesa, G., & Reimer, O. 2020, A&A, 637, A23, doi: 10.1051/0004-6361/202037442
  * Massi et al. (2001) Massi, M., Ribó, M., Paredes, J. M., Peracaula, M., & Estalella, R. 2001, A&A, 376, 217, doi: 10.1051/0004-6361:20010953
  * Mochol & Kirk (2014) Mochol, I., & Kirk, J. G. 2014, Astronomische Nachrichten, 335, 256, doi: 10.1002/asna.201312028
  * Moldón et al. (2011a) Moldón, J., Johnston, S., Ribó, M., Paredes, J. M., & Deller, A. T. 2011a, ApJ, 732, L10, doi: 10.1088/2041-8205/732/1/L10
  * Moldón et al. (2011b) Moldón, J., Ribó, M., & Paredes, J. M. 2011b, A&A, 533, L7, doi: 10.1051/0004-6361/201117764
  * Moldon et al. (2011) Moldon, J., Ribo, M., & Paredes, J. M. 2011, Astrophysics and Space Science Proceedings, 21, 1.  https://arxiv.org/abs/1101.4153
  * Morgan et al. (1955) Morgan, W. W., Code, A. D., & Whitford, A. E. 1955, ApJS, 2, 41, doi: 10.1086/190016
  * Moritani et al. (2018) Moritani, Y., Kawano, T., Chimasu, S., et al. 2018, PASJ, 70, 61, doi: 10.1093/pasj/psy053
  * Moritani et al. (2015) Moritani, Y., Okazaki, A. T., Carciofi, A. C., et al. 2015, ApJ, 804, L32, doi: 10.1088/2041-8205/804/2/L32
  * Nievas & VERITAS Collaboration (2021) Nievas, M., & VERITAS Collaboration. 2021, International Cosmic Ray Conference
  * Papitto et al. (2012) Papitto, A., Torres, D. F., & Rea, N. 2012, ApJ, 756, 188, doi: 10.1088/0004-637X/756/2/188
  * Paredes-Fortuny et al. (2015a) Paredes-Fortuny, X., Bosch-Ramon, V., Perucho, M., & Ribó, M. 2015a, A&A, 574, A77, doi: 10.1051/0004-6361/201424672
  * Paredes-Fortuny et al. (2015b) Paredes-Fortuny, X., Ribó, M., Bosch-Ramon, V., et al. 2015b, A&A, 575, L6, doi: 10.1051/0004-6361/201425361
  * Park & VERITAS Collaboration (2015) Park, N., & VERITAS Collaboration. 2015, in International Cosmic Ray Conference, Vol. 34, 34th International Cosmic Ray Conference (ICRC2015), 771\.  https://arxiv.org/abs/1508.07070
  * Parsons & Hinton (2014) Parsons, R. D., & Hinton, J. A. 2014, Astroparticle Physics, 56, 26, doi: 10.1016/j.astropartphys.2014.03.002
  * Rea & Torres (2011) Rea, N., & Torres, D. F. 2011, ApJ, 737, L12, doi: 10.1088/2041-8205/737/1/L12
  * Robertson et al. (2015) Robertson, D. R. S., Gallo, L. C., Zoghbi, A., & Fabian, A. C. 2015, MNRAS, 453, 3455, doi: 10.1093/mnras/stv1575
  * SALT Transient Program (2020) SALT Transient Program. 2020, personal communication (L. Townsend)
  * Scargle (1982) Scargle, J. D. 1982, ApJ, 263, 835, doi: 10.1086/160554
  * Shafter et al. (1986) Shafter, A. W., Szkody, P., & Thorstensen, J. R. 1986, ApJ, 308, 765, doi: 10.1086/164549
  * Skilton et al. (2009) Skilton, J. L., Pandey-Pommier, M., Hinton, J. A., et al. 2009, MNRAS, 399, 317, doi: 10.1111/j.1365-2966.2009.15272.x
  * Stellingwerf (1978) Stellingwerf, R. F. 1978, ApJ, 224, 953, doi: 10.1086/156444
  * Stoyanov et al. (2018) Stoyanov, K. A., Zamanov, R. K., & Iliev, I. K. 2018, The Astronomer’s Telegram, 11233, 1
  * Sushch & van Soelen (2017) Sushch, I., & van Soelen, B. 2017, ApJ, 837, 175, doi: 10.3847/1538-4357/aa62ff
  * Timmer & Koenig (1995) Timmer, J., & Koenig, M. 1995, A&A, 300, 707
  * Tody (1986) Tody, D. 1986, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 627, The IRAF Data Reduction and Analysis System, ed. D. L. Crawford, 733, doi: 10.1117/12.968154
  * Tokayer et al. (2021) Tokayer, Y., An, H., Halpern, J., et al. 2021, Contemporaneous Multi-Wavelength Campaign to Study HESS J0632+057s Distinctive Light Curve, to be published in ApJ
  * Torres et al. (2012) Torres, D. F., Rea, N., Esposito, P., et al. 2012, ApJ, 744, 106, doi: 10.1088/0004-637X/744/2/106
  * Uttley et al. (2003) Uttley, P., Edelson, R., McHardy, I. M., Peterson, B. M., & Markowitz, A. 2003, ApJ, 584, L53, doi: 10.1086/373887
  * Volkov et al. (2021) Volkov, I., Kargaltsev, O., Younes, G., Hare, J., & Pavlov, G. 2021, arXiv e-prints, arXiv:2103.04403.  https://arxiv.org/abs/2103.04403
  * Weekes et al. (2002) Weekes, T. C., Badran, H., Biller, S. D., et al. 2002, Astroparticle Physics, 17, 221, doi: 10.1016/S0927-6505(01)00152-9
  * Weng et al. (2021) Weng, S.-S., Pan, Z., Qian, L., et al. 2021, The Astronomer’s Telegram, 14297, 1
  * Yoneda et al. (2020) Yoneda, H., Makishima, K., Enoto, T., et al. 2020, Phys. Rev. Lett., 125, 111103, doi: 10.1103/PhysRevLett.125.111103
  * Yudin et al. (2017) Yudin, R. V., Potter, S. B., & Townsend, L. J. 2017, MNRAS, 464, 4325, doi: 10.1093/mnras/stw2676
  * Zamanov et al. (2021) Zamanov, R. K., Stoyanov, K. A., Mart, J., Marchev, V. D., & Nikolov, Y. M. 2021, Astronomische Nachrichten, 342, 531, doi: 10.1002/asna.202123856
  * Zamanov et al. (2016) Zamanov, R. K., Stoyanov, K. A., Martí, J., et al. 2016, A&A, 593, A97, doi: 10.1051/0004-6361/201628735
  * Zanin et al. (2013) Zanin, R., Carmona, E., J., S., et al. 2013, in Proc. of the 33th International Cosmic Ray Conference (ICRC), 773, doi: http://inspirehep.net/record/1412925/files/icrc2013-0773.pdf

##  Appendix A Orbital period determination

Literature provides a variety of methods for periodicity analysis of sparse
astronomical data (e.g. Graham et al., 2013). The performance of the different
techniques depend on quality, coverage, and shape of the given light curves.
There is no clear guidance which technique is best for a given data set. For
these reasons, the following methods are evaluated using Monte Carlo-generated
light curves: discrete correlation functions (DCF, Edelson & Krolik, 1988;
Robertson et al., 2015)777The implementation of the DCF method in the python
package pyDCF (https://github.com/astronomerdamo/pydcf) is used., correlation
analysis comparing the light curves with a binned-average light curve (PCC,
Malyshev et al., 2019), Lomb-Scargle periodograms (Lomb, 1976; Scargle,
1982)888The implementation of the Lomb-Scargle method in astropy is used
(https://docs.astropy.org/en/stable/timeseries/lombscargle.html), and phase
dispersion minimisation method (PDM, Stellingwerf, 1978)999The implementation
of the PDM method in the python package PyAstronomy
(https://github.com/sczesla/PyAstronomy) is used..

Figure 11:  Outcome of the analysis of 1000 Monte Carlo-generated light curves
using the following data sets as templates for the toy MC: Top: _Swift_ -XRT;
Bottom: VERITAS and H.E.S.S.. In all cases a true orbital period of 321 days
is assumed. Different techniques for periodicity analysis, as indicated in the
figure panels, have been applied to the same MC-generated light curves.
Results are given both as mean and 68% fiducial interval for all methods.

The Monte Carlo-generated light curves are based on the phase-binned average
profiles of the observed fluxes in X- or gamma rays and are generated assuming
the same distribution in time as the measurements, with flux values randomly
altered according to the corresponding uncertainties. As the phase-binned
averaged profile removes variability which is not due to statistical
uncertainties of the measurement, additional power-law noise following a
$(1/f)^{1.2}$ spectrum is added (Timmer & Koenig, 1995).

Figure 12:  Results of the periodicity analysis of _Swift_ -XRT (a-d) and
gamma-ray (e,f) measurements applying different methods: a) Phase dispersion
(PDM; _Swift_ -XRT); b) Discrete correlation coefficient (DCF; _Swift_ -XRT);
c) Pearson correlation coefficient (PCC; _Swift_ -XRT); d) Search for periods
other than the orbital period in the range 50 to 1000 days; Pearson
correlation coefficient (PCC; _Swift_ -XRT); the red-dashed line indicates the
99% containment level obtained from Monte Carlo sets of the XRT light curve
(to take spectral leakage into account, an orbital periodicity of 317.3 days
is assumed in the Monte Carlo sets); e) Phase dispersion (PDM; gamma-ray
energies); f) Pearson correlation coefficient (PCC; gamma-ray energies). All
coefficients are plotted as a function of orbital period. The green shaded
areas indicate the 68% fiducial interval obtained by the analysis of MC-
generated light curves (see text for details).

Figure 11 shows the results of the different techniques applied to 1000 Monte
Carlo-generated light curves based on the _Swift_ -XRT and gamma-ray data sets
assuming a true orbital period of 321 days. The search region is restricted to
a $\pm 30$ day-interval around the true orbital period. This prior is
necessary, as the PDM method tends towards reconstructing from a significant
fraction ($>20$%) of light curves an orbital period roughly half the true
orbital period. Applying the method of Lomb-Scargle to the generated light
curves leads to inconsistent results, as the largest peak in the periodograms
for most light curves is at half of the true period. This result is also
obtained when applying this method to the _Swift_ -XRT data. Folding the light
curve with such short periods leads to inconsistent solutions, as both the PDM
and Lomb-Scargle methods ignore its non-sinusoidal shape and in particular the
very low fluxes respectively non-detections of the system visible at orbital
phases of $\approx 0.4$ (assuming 317.3 days of orbital period).

All shown methods are able to reconstruct, for the majority of the MC light
curves, the true orbital period with an uncertainty of less than 2 days for
the _Swift_ -XRT and gamma-ray data set. The DCF method does not provide
reliable estimations when applied to the gamma-ray measurements or on the toy
MC based on the gamma-ray data. This is probably due to the sparsity of the
gamma-ray dataset and the larger uncertainties of the flux measurements.
Statistical uncertainties for the given datasets are derived from the 68%
fiducial intervals of the corresponding toy MC analysis. Systematic
uncertainties are derived from the largest difference to the expected values
(0.6 days) and the impact of the choice of the bin width on the calculation of
the averaged light curves. Bin widths from 0.025 to 0.1 in orbital phase are
tested and the largest difference is used to estimate the contribution to the
systematic uncertainties. The total systematic uncertainty for the orbital
period determination is estimated to be $\delta_{sys}^{Swift\
\mathrm{XRT}}=1.5$ days and $\delta_{sys}^{\mathrm{Gamma}}=2.5$ days.

Confidence limits on the correlation coefficients are obtained in a similar
fashion: 1000 toy MC light curves with similar data structure to the one
obtained by observations are used to calculate the 95 and 99% quantiles. Two
null hypotheses are distinguished: a constant flux is assumed for the
determination of the statistical significance of the correlation coefficients
for the orbital period analysis (see Figure 12, c)). For the search for
modulation periods other than the orbital period, toy MC light curves as
described above, including the average observed orbital modulation pattern,
are generated and analysed with the different period determination methods
(see Figure 12, d)).

Figure 12 shows the results of application of all three techniques on the
_Swift_ -XRT and gamma-ray light curve measurements (see Figure 2). A summary
of all obtained orbital periods together with those reported in the literature
is given in Table 3.

##  Appendix B Impact of orbital period uncertainty on phase-folded light
curves

Uncertainties on the orbital period determination might lead to significant
differences in the shape of the phase-folded light curves given the long total
observation time of $\approx$15 years for gamma-ray and $\approx$10 years for
X-ray measurements. To test this, four different orbital periods are assumed:
$P_{+}=319.5$ days, $P_{-}=315.1$, $P_{M}=313$ days, and $P_{H\alpha}=308$
days. $P_{+}$ and $P_{-}$ correspond to a change of the orbital period by 1
$\sigma$ statistical error plus the systematic uncertainty; $P_{M}$ and
$P_{H\alpha}$ correspond to the solutions presented in Moritani et al. (2018).

Figures 13 and 14 show the impact of the variation of the assumed orbital
period: the three most prominent features of a peak around phases 0.3, a
minimum around phases 0.4, and a second maximum region around phases 0.6 are
clearly visible in all cases. This shows that the choice of $P_{orbit}=317.3$
does not influence the discussion of the physical properties presented in this
work, except for the shortest period of 308 days. It should be noted that
assuming such a short orbital period, as derived from optical H$\alpha$
observations (Moritani et al., 2018), would change this picture significantly,
as most of the discussed features disappear (Figures 13 (d) and 14 (d).

Figure 13:  X-ray (0.3–10 keV) light curves as function of orbital phase
assuming orbital periods as indicated below the figures. For further details,
see Figure 3 and text. Figure 14:  Gamma-ray ($>$350 GeV) light curves as
function of orbital phase assuming orbital periods as indicated below the
figures. For further details, see Figure 3 and text.

##  Appendix C Light curves per orbital cycle

Figure 15:  X-ray (0.3–10 keV) light curves as function of orbital phase for
each of the observed orbital cycles. Orbits are numbered following the start
of the gamma-ray observations (no X-ray observations are available for the
first four orbits; empty panels are shown for easier comparisons with Figure
16). The MJD given in each panel indicates the start of each orbit. An orbital
period of 317.3 days is assumed.Vertical lines show statistical uncertainties;
note that these are smaller than the marker size for all instruments but
_Swift_ -XRT. The thin blue line and gray-shaded band in each canvas indicates
the average _Swift_ -XRT light curve and its 68% containment region calculated
from all measurements and smoothed by applying a cubic spline fit.  Figure 16:
Gamma-ray ($>$ 350 GeV) light curves as function of orbital phase for each of
the observed orbital cycles (see Figure 15 for further details). The thin blue
line and gray-shaded band in each canvas indicates the average gamma-ray light
curve and its 68% containment region calculated from all measurements and
smoothed by applying a cubic spline fit.

The gamma-ray observations with H.E.S.S., MAGIC, and VERITAS sum up to a total
observation time of $\approx 450$ h spanning 18 orbits of HESS J0632+057
covering the period of 2004–2019. The X-ray observations by the Swift _XRT_ ,
Chandra, XMM-Newton, Suzaku, and NuSTAR are obtained during 14 orbits of the
binary system. Despite this large amount of data, there is no good coverage
along most of these orbital cycles due to observational constraints and the
long binary period of 317 days. This is illustrated in figures 15 and 16,
which shows the light curve per orbital cycle both at X-ray and gamma-ray
energies assuming an orbital period of 317.3 days.

##  Appendix D Contemporaneous spectral energy distributions

Contemporaneous X-ray and gamma-ray spectral energy distributions for 38
periods of VERITAS and _Swift_ -XRT observations are available through the
Zenodo data repository101010https://doi.org/10.5281/zenodo.5157848.


# Sequence-level self-learning with multiple hypotheses

###### Abstract

In this work, we develop new self-learning techniques with an attention-based
sequence-to-sequence (seq2seq) model for automatic speech recognition (ASR).
For untranscribed speech data, the hypothesis from an ASR system must be used
as a label. However, the imperfect ASR result makes unsupervised learning
difficult to consistently improve recognition performance especially in the
case that multiple powerful teacher models are unavailable. In contrast to
conventional unsupervised learning approaches, we adopt the _multi-task
learning_ (MTL) framework where the $n$-th best ASR hypothesis is used as the
label of each task. The seq2seq network is updated through the MTL framework
so as to find the common representation that can cover multiple hypotheses. By
doing so, the effect of the _hard-decision_ errors can be alleviated. We first
demonstrate the effectiveness of our self-learning methods through ASR
experiments in an accent adaptation task between the US and British English
speech. Our experiment results show that our method can reduce the WER on the
British speech data from 14.55% to 10.36% compared to the baseline model
trained with the US English data only. Moreover, we investigate the effect of
our proposed methods in a federated learning scenario.

Index Terms: unsupervised learning, self-learning, encoder-decoder, multi-task
learning, end-to-end acoustic modeling

##  1 Introduction

Unsupervised learning of an acoustic model (AM) has a long history from the
Bayesian model to the deep neural network (DNN) system in the field of
automatic speech recognition (ASR) [1, 2, 3, 4, 5]. It is typically done by
transferring knowledge from stronger teacher model(s) to the student model [6,
7, 8, 9, 10, 11, 12] or adapting the seed model pre-trained with a sufficient
amount of labeled data [1, 2, 3, 13, 14, 15]. Because of no need of heavy
decoding processes with teacher models, the latter self-training approach is
perhaps preferable in many situations such as on-device personalization and
federated learning scenarios [16, 17, 18, 19]. In this work, we thus address
the self-training task for AM without an additional bigger teacher model and
parallel adaptation data111The definitions of unsupervised learning and self-
learning vary in different applications and fields. For presentation purposes,
we will refer any model optimization methods without manual transcripts as
unsupervised learning even if the initial seed model is pre-trained with
labeled data. We will say self-learning is a special class of unsupervised
learning techniques which optimizes the seed model without another teacher
model..

The optimization criteria of unsupervised learning for DNN-based ASR systems
can fall into two categories: frame-level [8, 10, 11, 12], and sequence-level
criteria [7, 15, 9]. It is well known that direct sequence-level optimization
[15] and sequence training from the model initialized with the frame-level
cross-entropy objective [10, 11] provide better recognition accuracy. Many
successful results have been reported by adapting the hybrid model of the DNN
and hidden Markov model (HMM) in an unsupervised way.

In the hybrid DNN-HMM systems, acoustic, pronunciation, and language models
have been trained separately, each with a different objective. While such an
approach can use each type of training data efficiently, this disjoint
modeling method leads to the sub-optimal solution, especially in terms of
model size. This issue has been addressed by designing the _end-to-end_
network, which directly generates a word or character sequence from a sequence
of speech features. Two popular approaches in this area are perhaps recurrent
NN transducer (RNN-T) [20] and attention-based encoder-decoder (AED) network
[21, 22, 23]. Unlike connectionist temporal classification (CTC) [24], both
methods do not make an unreasonable assumption for ASR that the label outputs
are conditionally independent of each other. It is shown in [25] that the AED
has the potential of outperforming the RNN-T arguably because of a more
flexible assumption on the alignment between input and output, given a large
number of training samples [26, 27].

Based on the reasons above, we focus on the development of unsupervised
learning methods for the AED network. For that, knowledge distillation (KD)
[28, 29, 9] is a promising approach. In particular, sequence-level KD [29, 9]
is more suitable for the AED network to improve modeling capability for
sequential input and output. Kim and Rush indeed showed in [29] that the
sequence-level KD technique can achieve better machine translation performance
with lower complexity. Sequence-level KD is also applied to the ASR tasks:
model compression [9] and domain adaptation with a slight modification of
information injected to the student [30]. Notice that the unsupervised
learning method proposed in [30] requires parallel data for the teacher and
student models, a pair of clean and noisy data. Again, the use of the bigger
teachers or additional parallel adaptation data is not possible in some
application scenarios, such as a federated learning scenario [18].

In contrast to prior work on unsupervised learning for the AED [9, 30], we
develop the new self-learning techniques which neither require gigantic
teacher models nor the parallel data. In such a scenario, we must use the ASR
output of the seed model as the target labels. The use of the inaccurate ASR
result as the ground truth may degrade recognition accuracy after adaptation.
In order to alleviate such a noisy transcript issue, we use the $N$-best
hypotheses from the ASR output for adaptation. More specifically, we optimize
the AED with multiple objective functions associated with the $N$-best
hypotheses. We implement the optimization algorithm with the multi-label
learning (MLL) or multi-task learning (MTL) framework [31, 32]. Our techniques
proposed here are evaluated through two kinds of ASR tasks: an accent
adaptation task and a federated learning situation.

The rest of the paper is organized as follows. Section 2 briefly reviews the
sequence-to-sequence model with the AED and its optimization criterion in
supervised and unsupervised manners. Section 3 describes our unsupervised
learning algorithms for the attention-based encoder-decoder network. Section 4
describes the ASR experiments in two kinds of tasks: accent adaptation and
federated transfer learning tasks. In section 5, we conclude this work.

##  2 Background

###  2.1 Attention-based encoder-decoder (AED)

Let $\mbox{$\mathbf{x}$}=[x_{1},...,x_{K}]$ and
$\mbox{$\mathbf{y}$}=[y_{1},...,y_{M}]$ be the feature input and target word
embeddings with $K$ and $M$ being the source and target lengths, respectively.
ASR involves finding the most probable target sequence given the feature
input:

$\hat{\mbox{$\mathbf{y}$}}=\mathop{\rm
argmax}_{\mbox{$\mathbf{y}$}\in\mathcal{Y}}\log
p(\mbox{$\mathbf{y}$}|\mbox{$\mathbf{x}$})\vspace{-0.1cm}$ (1)

where $\mathcal{Y}$ is the set of all possible sequences. As shown in Figure
1, $p(\mbox{$\mathbf{y}$}|\mbox{$\mathbf{x}$})$ is modeled with an encoder, a
decoder and an attention component. The encoder network reads the feature
sequence and the decoder network produces a distribution over the target
words, one word embedding at a time given the source [23, 26]. We employ the
attention architecture from [21] by replacing the scoring function with the
ReLU unit.

Figure 1: Attention-based encoder-decoder (AED) network

###  2.2 Supervised learning

We briefly review supervised learning for the sequence-to-sequence model
before formulating the unsupervised learning problem. First, consider the
sequence-level distribution specified with the model parameters $\theta$ over
all possible sequences $\mbox{$\mathbf{y}$}\in\mathcal{Y}$ in the log domain,

$\log p(\mbox{$\mathbf{y}$}|\mbox{$\mathbf{x}$})=\sum_{i=1}^{M}\log
p(y_{i}|\mbox{$\mathbf{x}$},\mbox{$\mathbf{y}$}_{<i};\theta)\vspace{-0.1cm}$
(2)

where $\mbox{$\mathbf{y}$}_{<i}$ is either the ground truth of the previous
characters or sampled from the model with a certain sampling rate in order to
make the model robust against prediction errors during inference [23]. The
sampling rate can be constant [23] or scheduled [33]. For the experiment, we
use the scheduled sampling method. Sequence-level training then involves
estimating $\theta$ that maximizes the log-likelihood (2) over all the
sequences.

###  2.3 Unsupervised learning

The self-learning framework of the attention-based sequence-to-sequence model
is very similar with sequence-level knowledge distillation (KD) [29]. The only
difference is that there is no better teacher model which transfers knowledge
to the seed model in the self-learning scenario [34]. In such a scenario, the
seed model must adapt itself to new data without knowledge distillation from
the more powerful teacher model; the sequence distribution of target words can
be very noisy. The erroneous labels can degrade accuracy of the seed model
severely [5, 34]. In this section, we describe how we approximate the
objective function in the case of the unsupervised learning scenario.

First, let us consider the case of the semi-supervised setting with sequence-
level KD. As shown in [29], the distribution from the data is replaced with a
probability distribution derived from another teacher model
$q(\mbox{$\mathbf{y}$}|\mbox{$\mathbf{x}$})$. The loss function of sequence-
level KD can be written as

$\mathcal{L}_{\text{SEQ-
KD}}=-\sum_{\mbox{$\mathbf{y}$}\in\mathcal{Y}}q(\mbox{$\mathbf{y}$}|\mbox{$\mathbf{x}$})\log
p(\mbox{$\mathbf{y}$}|\mbox{$\mathbf{x}$})\vspace{-0.1cm}$ (3)

Calculating (3) is computationally intractable. Kim and Rush showed in [29]
that the loss (3) could be well approximated with the best hypothesis from the
beam search output with the teacher model in a machine translation task. With
the first best hypothesis $\hat{\mbox{$\mathbf{y}$}}_{\text{T,1}}$, they
approximated (3) as

$\mathcal{L}_{\text{SEQ-KD}}\approx-\log
p(\mbox{$\mathbf{y}$}=\hat{\mbox{$\mathbf{y}$}}_{\text{T,1}}|\mbox{$\mathbf{x}$})\vspace{-0.1cm}$
(4)

It is worth noting that the same algorithm was also applied to an ASR model
compression task in [9] although their improvement was not as significant as
that reported in [29].

Now consider the case of self-learning. Again, we only have the seed model in
this case. Thus, we use the best hypothesis of the seed model
$\hat{\mbox{$\mathbf{y}$}}_{\text{S,1}}$ instead of that of the teacher model
$\hat{\mbox{$\mathbf{y}$}}_{\text{T,1}}$. It is straightforward to express the
loss approximation function:

$\mathcal{L}_{\text{SEQ-KD}}\approx-\log
p(\mbox{$\mathbf{y}$}=\hat{\mbox{$\mathbf{y}$}}_{\text{S,1}}|\mbox{$\mathbf{x}$})\vspace{-0.1cm}$
(5)

The model can still be improved with self-learning since the seed model learns
features from unseen data. However, the seed model needs to provide reasonably
good accuracy on the unseen data. Otherwise, the model may degrade due to the
erroneous label estimates.

##  3 Proposed self-learning method

###  3.1 Loss function approximation with multiple hypotheses

Figure 2: MT network with encoder-decoder shared.

Figure 3: MT network with encoder shared.

Figure 4: Flow chart of unsupervised training with the $N$-best hypotheses

In many cases, the best ASR output contains errors. Thus, approximating (3)
with the only best hypothesis may degrade recognition performance after
unsupervised training. Instead, we use the $N$-best hypotheses weighted based
on the ASR confidence score for approximation. First, let us denote the weight
normalized with the Softmax function for the $n$-th best candidate
$\mbox{$\mathbf{y}$}_{S,n}$ as

$q_{n}(\mbox{$\mathbf{y}$}_{\text{S,n}})=\frac{\exp(s_{n}/T)}{\sum_{i=1}^{N}\exp(s_{i}/T)}$
(6)

where $s_{n}$ is the $n$-th best confidence score and $T$ is a temperature
which controls a probability distribution over the classes [28]. We set $T=1$
for the experiments. We then rewrite the loss function (3) as

$\mathcal{L}_{\text{SEQ-
KD}}\approx-\sum_{n=1}^{N}q_{n}(\hat{\mbox{$\mathbf{y}$}}_{\text{S,n}})\log
p(\mbox{$\mathbf{y}$}=\hat{\mbox{$\mathbf{y}$}}_{\text{S,n}}|\mbox{$\mathbf{x}$})\vspace{-0.1cm}$
(7)

Minimizing this loss function can be done by creating multiple labels
associated with weights for one utterance. We refer to this method as multi-
label learning (MLL) in this paper. The loss function can be also optimized
through multi-task learning (MTL) [31, 32] as it will be explained in the next
section.

###  3.2 Optimization with multi-task learning

Optimization of the loss function (7) can be done with the multi-task (MT)
network where each separate Softmax layer is associated with the loss of the
$n$-th hypothesis. We consider two network architectures in this work. Figure
4 shows the first architecture proposed here. In the architecture shown in
Figure 4, the parameters of the encoder, attention and decoder network are
shared during the optimization process. Only the output layer is separated for
each task associated with the term in (7). The weights of the new branch for
the $n$-th best hypothesis are initialized with those of the main task.

Figure 4 shows the second architecture investigated in this work. In contrast
to the architecture in Figure 4, the second network architecture shares the
encoder and attention module only; the set of the decoder and the output layer
is separated for each task. After self-training, the encoder is expected to
generate a better acoustic embedding that covers the multiple hypotheses. This
can be viewed as acoustic feature-space adaptation with the AED network.

Figure 4 shows a flow chart of our unsupervised training scheme with the
$N$-best hypotheses. As shown in Figure 4, we first decode the sequence of the
speech features with the beam search in order to obtain the $N$-best
hypotheses, $\hat{\mbox{$\mathbf{y}$}}_{\text{S,1}}$ $\cdots$
$\hat{\mbox{$\mathbf{y}$}}_{\text{S,N}}$. Then, the seed model is adapted with
the $N$ hypotheses. This can be done by either multi-label learning (MLL) or
MTL framework, as described before. The normalized probability computed with
(6) is used for weighting the loss of each utterance in MLL and the task in
MTL. After adopting the model with the $N$-th hypotheses, we train the model
with the best hypothesis only. The only single task branch is used in the case
of MTL. These steps will be repeated until there is no improvement in a
character error rate on the validation data. During decoding, we can remove
the sub-task branches and keep the main task branch only.

##  4 ASR experiment

In this section, we describe the ASR experiments in order to compare the
unsupervised learning techniques with the AED.

###  4.1 Accent Adaptation Task

First, we addressed the English accent adaptation task from American English
to British English. The data used for the experiment were collected through
the live traffic and stored with the sampling rate of 16 kHz. For the
experiments, we train the seed model with approximately 75,000 hours of
American English speech data. About 1,000 hours of the British English speech
data were used for adaptation. The adaptation data are treated as
untranscribed in the case of unsupervised learning experiments. The British
Cortana corpus was used as a test set. The test data were also collected
through the live traffic and consists of 14,000 utterances spoken by hundreds
of users.

The feature extraction front-end computes 40 log mel-filterbank energy
features at a frame rate of 10 milliseconds (msec), concatenates three
adjacent feature frames, and then downsamples it into a 30 msec frame rate.
Each feature was rescaled to have zero mean and unit variance over the
training set. The global mean and variance calculated with the training data
were applied to adaptation and evaluation. In order to augment adaptation
data, we randomly scale the amplitude of a signal, change the speed in the
same manner implemented in Sox software [35] and perform the spec-augment on
the feature domain [36]. The network used here has 6 layers of 1024
bidirectional GRU nodes in the encoder and 2 layers of 1024 unidirectional GRU
nodes in the decoder. All dropout types were applied with $0.1$ dropout
probability. We used label smoothing and Adam optimizer in the experiments.

Table 1 shows the word error rates (WERs) obtained with different models, the
seed model trained with American English data, British AMs adapted in the
supervised and unsupervised learning manners. In the self-learning scenario,
we performed 4 types of sequence-level self-knowledge distillation methods,
conventional training with the 1-best hypothesis (5), multi-label learning
(MLL) with the 4-best hypotheses (7), multi-task learning (MTL) with the
attention, encoder and decoder (AED) shared and that with the attention and
encoder (AE) shared. For unsupervised MTL, we use the 4-best hypotheses. It is
clear from Table 1 that the WER of the American English model can be reduced
from 14.55% to 7.8% by supervised learning with the manual transcripts. It is
also clear from Table 1 that the WER of the seed model can be reduced by the
self-learning methods, although the improvement is not as significant as that
obtained with supervised learning. Table 1 also shows that the use of the
4-best hypotheses provides better recognition accuracy than the case of the
1-best hypothesis only. The WERs of the 1-best and N-best results on the
adaptation data obtained with the seed model were 22.6% and 13.7%,
respectively. Thus, we consider using multiple hypotheses to provide better
recognition accuracy. We did not observe a significant difference among
different self-learning methods with the 4-best hypotheses after decoding the
adaptation data twice, but it is worth mentioning that MTL with the shared AE
network provided the lowest WER with the first decoding pass.

Adaptation Type | Adaptation Method | WER(%)
---|---|---
Seed model | US baseline model | 14.55
Supervised learning | Training with transcription | 7.8
Self-learning | Training with the 1-best hypo | 12.09
MLL with 4-best hypos | 10.66
MTL with shared AED | 10.83
MTL with shared AE | 10.36
Table 1: Word error rates (WER) on the British Cortana data for each training
method

###  4.2 Librispeech Task

We also conducted ASR experiments with the Librispeech (LS) data [37] for this
study. For the LS task, we used the same front-end as described in Section 4.1
but slightly different model configuration. The AED network used here consists
of 6 layers of 1024 bidirectional LSTM nodes and 2 layers of 1024
unidirectional LSTM nodes. For the LS task, we used byte-pair encoding (BPE)
[38] to create 16,000 subword units. The optimizer settings are the same as
described in [18].

The LS corpus contains approximately 1000 hours of read speech for training.
We split the whole LS training set into two subsets. A seed model was first
trained with the first half of the training set. For adaptation, we further
divide the second half dataset into the per-speaker dataset. The adaptation
dataset contains approximately 1100 speakers. Since the seed model does not
see the LS vocabulary and word sequences of the adaptation dataset, the ASR
result on the adaptation data will not be perfect; self-learning has to be
done with noisy labels. The seed model was then adapted with the speaker-
partitioned dataset through our federated learning simulator [18, 19].
Briefly, the federated learning algorithm iteratively performs the following
processes on a server and clients: (1) broadcasting the global (seed) model
from the server to the clients, (2) decoding speech with the global model at
each client, (3) adapting the global model with each client’s data and the ASR
hypotheses, (4) encrypting each client’s model so as to ensure privacy222We
did not use encryption here for the sake of simplicity although it is
straightforward to apply multi-party computation (MPC) secure algorithms or
differential privacy techniques such as [17]., (5) aggregating the linearly-
weighted models from a pool of clients, and (6) optimizing the global model
through the aggregated gradients. Those processes are repeated until the loss
on the validation data is converged. We only decode speech every 512
aggregation steps for unsupervised training. We found that this could avoid
training divergence due to inaccurate ASR results with an intermediate model
trained halfway. Notice that the adapted model is only uploaded to the server
while privacy data such as speech and transcripts are secured in the client.
The adapted models will be sent to the server after processing all the
utterances per speaker on each client. This saves a tremendous amount of
network bandwidths.

Table 2 shows the WERs on the LS test clean set obtained with different models
in the same way as Table 1. We performed supervised and unsupervised learning
algorithms on the seed model with the second half of the training data. The
last row in Table 2 shows the upper bound WER in the case that the model is
trained from scratch with the full amount of the LS training data. It is clear
from Table 2 that the WER of the seed model can be significantly reduced by
supervised-training with the rest 50% of the training data from 5.66%to 4.40%.
It is also clear from Table 2 that recognition accuracy can still be improved
by self-learning with MTL even when the seed model does not know the
vocabulary and word sequences of the adaptation data. We observed that the
validation loss diverged in MLL, which very likely led to recognition accuracy
degradation. The results also show that the convergence performance is
improved by separating more layers relevant to the language model, the decoder
and output layers, for each hypothesis. This suggests that the $N$ best
hypotheses contain little language information shareable over the decoders and
output layers. This may be because language knowledge of the seed model on the
adaptation dataset is poor. In order to check how much the N-best results can
potentially help recognition accuracy, we decoded the second half of the LS
training data with our LS seed model. The WERs of the 1-best and N-best
results on the LS adaptation data were 10.4% and 7.62%. The difference between
the 1-best and N-best results was not as significant as that in the case of
the accent adaptation task. We assume that many word errors in this LS task
cannot be recovered by simply increasing the $N$-best hypotheses because the
seed model does not observe the language of the adaptation dataset. On the
other hand, the seed model used in the British Cortana task should cover the
vocabulary and grammar for the adaptation data very well because it is trained
with a large amount of the American English data containing the Cortana task.
Of course, there are some different word usages between the American and
British English. We are, however, led to conclude that the short phrase ASR
task such as the Cortana task was probably not affected significantly.

Adaptation Type | Adaptation Method | WER(%)
---|---|---
Seed model | Trained from scratch with 50% of data | 5.66
Supervised learning | Training with the rest 50% of data | 4.40
Self-learning | Training with the 1-best hypo | 5.58
MLL with 4-best hypos | 5.72
MTL with shared AED | 5.55
MTL with shared AE | 5.22
Reference | The whole LS training data | 4.00
in centralized training |
Table 2: WERs on the LS test clean dataset for each training method in the
privacy preserving scenario (where the seed model has no language information
on the adaptation data)

##  5 Conclusions

In this work, we have investigated the ASR accuracy of self-learning methods
with the attention-based sequence-to-sequence model. We proposed the new self-
learning methods with the $N$ best hypotheses; the seed model is adapted with
MLL or MTL in an end-to-end manner. We conducted two kinds of ASR experiments:
the English accent adaptation task and the federated-learning scenario. In the
accent adaptation task, our experiment results show that our best self-
learning method can reduce the WER on the British speech data from 14.55% to
10.36% compared to the baseline model trained with the US English data only.
In the federated-learning scenario, the improvement on the LS dataset was not
as significant as that on the Cortana dataset. In contrast to the Cortana
task, the seed model used in the LS task does not contain language information
on the adaptation data well. Thus, misrecognition results may not be entirely
recoverable by increasing the N-best candidates. We plan to combine our AM
self-learning and the language data augmentation method a text-to-speech
module [39]

##  6 Acknowledgements

The authors would like to thank Masaki Itagaki, Heiko Rahmel, Naoyuki Kanda,
Lei He, Ziad Al Bawab, Jian Wu, and Xuedong Huang for their project support
and technical discussions.

## References

  * [1] P. C. Woodland, D. Pye, and M. J. F. Gales, “Iterative unsupervised adaptation using maximum likelihood linear regression,” in _The 4th International Conference on Spoken Language Processing, Philadelphia, PA, USA, October 3-6, 1996_. ISCA, 1996.
  * [2] G. Zavaliagkos and T. Colthurst, “Utilizing untranscribed training data to improve performance,” in _DARPA Broadcast News Transcription and Understanding Workshop, Landsdowne_ , 1998, pp. 301–305.
  * [3] L. Lamel, J. luc Gauvain, and G. Adda, “Lightly supervised and unsupervised acoustic model training,” _Computer Speech and Language_ , vol. 16, pp. 115–129, 2002.
  * [4] F. Wessel and H. Ney, “Unsupervised training of acoustic models for large vocabulary continuous speech recognition,” _IEEE Transactions on Speech and Audio Processing_ , vol. 13, no. 1, pp. 23–31, 2005.
  * [5] Y. Huang, Y. Wang, and Y. Gong, “Semi-supervised training in deep learning acoustic model,” in _Proc. Interspeech_ , September 2016.
  * [6] J. Li, R. Zhao, J.-T. Huang, and Y. Gong, “Learning small-size dnn with output-distribution-based criteria,” in _Proc. Interspeech_ , September 2014\.
  * [7] J. H. M. Wong and M. J. F. Gales, “Sequence student-teacher training of deep neural networks,” in _Proc. Interspeech_ , September 2016.
  * [8] J. Li, M. L. Seltzer, X. Wang, R. Zhao, and Y. Gong, “Large-scale domain adaptation via teacher-student learning,” in _Proc. Interspeech_ , 2017, pp. 2386–2390.
  * [9] R. M. Mun’im, N. Inoue, and K. Shinoda, “Sequence-level knowledge distillation for model compression of attention-based sequence-to-sequence speech recognition,” _CoRR_ , vol. abs/1811.04531, 2018.
  * [10] S. H. K. Parthasarathi, N. Sivakrishnan, P. Ladkat, and N. Strom, “Realizing petabyte scale acoustic modeling,” _IEEE Journal on Emerging and Selected Topics in Circuits and Systems_ , vol. 9, no. 2, pp. 422–432, 2019.
  * [11] L. Mosner, M. Wu, A. Raju, S. H. K. Parthasarathi, K. Kumatani, S. Sundaram, R. Maas, and B. Hoffmeister, “Improving noise robustness of automatic speech recognition via parallel data and teacher-student learning,” in _Proc. ICASSP_ , 2019, pp. 6475–6479.
  * [12] Z. Meng, J. Li, Y. Zhao, and Y. Gong, “Conditional teacher-student learning,” _CoRR_ , vol. abs/1904.12399, 2019.
  * [13] J. Huang, E. Marcheret, K. Visweswariah, V. Libal, and G. Potamianos, _The IBM Rich Transcription 2007 Speech-to-Text Systems for Lecture Meetings_. Berlin, Heidelberg: Springer, 2007, vol. 4625.
  * [14] K. Kumatani, T. Arakawa, K. Yamamoto, J. W. McDonough, B. Raj, R. Singh, and I. Tashev, “Microphone array processing for distant speech recognition: Towards real-world deployment,” in _Proc. APSIPA ASC_ , 2012.
  * [15] V. Manohar, P. Ghahremani, D. Povey, and S. Khudanpur, “A teacher-student learning approach for unsupervised domain adaptation of sequence-trained ASR models,” in _2018 IEEE Spoken Language Technology Workshop, SLT 2018, Athens, Greece, December 18-21, 2018_. IEEE, 2018, pp. 250–257.
  * [16] J. Konecny, B. McMahan, and D. Ramage, “Federated Optimization: Distributed Optimization Beyond the Datacenter,” _arXiv preprint arXiv:1511.03575v1_ , 2015.
  * [17] B. Jayaraman, L. Wang, D. Evans, and Q. Gu, “Distributed learning without distress: Privacy-preserving empirical risk minimization,” in _Advances in Neural Information Processing Systems 31_ , S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, Eds. Curran Associates, Inc., 2018, pp. 6343–6354.
  * [18] D. Dimitriadis, K. Kumatani, R. Gmyr, Y. Gaur, and S. E. Eskimez, “A federated approach in training acoustic models,” in _Proc. Interspeech_ , 2020.
  * [19] ——, “Federated transfer learning with dynamic gradient aggregation,” vol. ArXiv: 2008.02452, 2020.
  * [20] A. Graves, “Sequence transduction with recurrent neural networks,” _CoRR_ , vol. abs/1211.3711, 2012.
  * [21] J. K. Chorowski, D. Bahdanau, D. Serdyuk, K. Cho, and Y. Bengio, “Attention-based models for speech recognition,” in _Advances in Neural Information Processing Systems 28_ , C. Cortes, N. D. Lawrence, D. D. Lee, M. Sugiyama, and R. Garnett, Eds. Curran Associates, Inc., 2015, pp. 577–585.
  * [22] D. Bahdanau, J. Chorowski, D. Serdyuk, P. Brakel, and Y. Bengio, “End-to-end attention-based large vocabulary speech recognition,” in _ICASSP_ , 2016, pp. 4945–4949.
  * [23] W. Chan, N. Jaitly, Q. V. Le, and O. Vinyals, “Listen, attend and spell,” _CoRR_ , vol. abs/1508.01211, 2015.
  * [24] A. Graves, S. Fernández, and F. Gomez, “Connectionist temporal classification: Labelling unsegmented sequence data with recurrent neural networks,” in _In Proceedings of the International Conference on Machine Learning, ICML 2006_ , 2006, pp. 369–376.
  * [25] R. Prabhavalkar, K. Rao, T. N. Sainath, B. Li, L. Johnson, and N. Jaitly, “A comparison of sequence-to-sequence models for speech recognition,” in _Proc. Interspeech_ , 2017, pp. 939–943.
  * [26] E. Battenberg, J. Chen, R. Child, A. Coates, Y. Gaur, Y. Li, H. Liu, S. Satheesh, A. Sriram, and Z. Zhu, “Exploring neural transducers for end-to-end speech recognition,” in _2017 IEEE Automatic Speech Recognition and Understanding Workshop, ASRU 2017, Okinawa, Japan, December 16-20, 2017_. IEEE, 2017, pp. 206–213.
  * [27] C.-C. Chiu, W. Han, Y. Zhang, R. Pang, S. Kishchenko, P. Nguyen, A. Narayanan, H. Liao, S. Zhang, A. Kannan, R. Prabhavalkar, Z. Chen, T. N. Sainath, and Y. Wu, “A comparison of end-to-end models for long-form speech recognition,” 2019, pp. 889–896.
  * [28] G. E. Hinton, O. Vinyals, and J. Dean, “Distilling the knowledge in a neural network,” _ArXiv_ , vol. abs/1503.02531, 2015.
  * [29] Y. Kim and A. M. Rush, “Sequence-level knowledge distillation,” _CoRR_ , vol. abs/1606.07947, 2016.
  * [30] Z. Meng, J. Li, Y. Gaur, and Y. Gong, “Domain adaptation via teacher-student learning for end-to-end speech recognition,” in _ASRU_. IEEE, December 2019.
  * [31] R. Caruana, “Multitask learning,” Ph.D. dissertation, Carnegie Mellon University, Pittsburgh, PA 15213, 9 1997.
  * [32] S. Ruder, “An overview of multi-task learning in deep neural networks,” _ArXiv_ , vol. abs/1706.05098, 2017.
  * [33] S. Bengio, O. Vinyals, N. Jaitly, and N. Shazeer, “Scheduled sampling for sequence prediction with recurrent neural networks,” _CoRR_ , vol. abs/1506.03099, 2015.
  * [34] Q. Xie, E. H. Hovy, M.-T. Luong, and Q. V. Le, “Self-training with noisy student improves imagenet classification,” _ArXiv_ , vol. abs/1911.04252, 2019.
  * [35] “Sox - sound exchange, audio manipulation tool,” Available at http://sox.sourceforge.net/ (2020/04/28).
  * [36] D. S. Park, W. Chan, Y. Zhang, C.-C. Chiu, B. Zoph, E. D. Cubuk, and Q. V. Le, “SpecAugment: A Simple Data Augmentation Method for Automatic Speech Recognition,” in _Proc. Interspeech_ , 2019, pp. 2613–2617.
  * [37] V. Panayotov, G. Chen, D. Povey, and S. Khudanpur, “Librispeech: An asr corpus based on public domain audio books,” _Proc. ICASSP_ , pp. 5206–5210, 2015\.
  * [38] R. Sennrich, B. Haddow, and A. Birch, “Neural machine translation of rare words with subword units,” in _Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics_. Berlin, Germany: Association for Computational Linguistics, Aug. 2016, pp. 1715–1725.
  * [39] N. Rossenbach, A. Zeyer, R. Schlüter, and H. Ney, “Generating synthetic audio data for attention-based speech recognition systems,” in _Proc. ICASSP_ , 2020, pp. 7069–7073.


# Bayesian tensor factorization for predicting clinical outcomes using
integrated human genetics evidence

Onuralp Soylemez

###### Abstract

The approval success rate of drug candidates is very low with the majority of
failure due to safety and efficacy. Increasingly available high dimensional
information on targets, drug molecules and indications provides an opportunity
for ML methods to integrate multiple data modalities and better predict
clinically promising drug targets. Notably, drug targets with human genetics
evidence are shown to have better odds to succeed. However, a recent tensor
factorization-based approach found that additional information on targets and
indications might not necessarily improve the predictive accuracy. Here we
revisit this approach by integrating different types of human genetics
evidence collated from publicly available sources to support each target-
indication pair. We use Bayesian tensor factorization to show that models
incorporating all available human genetics evidence (rare disease, gene
burden, common disease) modestly improves the clinical outcome prediction over
models using single line of genetics evidence. We provide additional insight
into the relative predictive power of different types of human genetics
evidence for predicting the success of clinical outcomes.

Machine Learning, Drug Discovery, Human Genetics, Bayesian Tensor
Factorization, CompBio, ICML

##  1 Motivation

The approval success rate of drug candidates is less than 10% with most of
clinical trial failures attributed to safety concerns and a lack of clinical
efficacy (Cook et al., 2014). Increasingly available high dimensional
information on targets, drug molecules and indications provides an opportunity
for ML methods to better predict clinically promising drug targets. Notably,
drug targets with human genetics evidence are twice more likely to be approved
[(Nelson et al., 2015), (King et al., 2019)] and the most recent drug
approvals from FDA corroborate this strong trend (Ochoa et al., 2022).
However, a recent tensor factorization-based approach from (Yao et al., 2019)
found that additional information on targets and indications might not
necessarily improve the predictive accuracy underscoring the importance of
feature selection and evidence data quality. Here we revisit this approach by
integrating different lines of human genetics evidence collated from publicly
available sources and assess the relative predictive performance of models
incorporating different types of human genetics evidence.

##  2 Data curation

Open Targets Platform curates and maintains target-disease evidence by
harmonizing information on genetic diseases, genetic variation, clinical trial
outcomes, gene expression data and biomedical literature (Ochoa et al., 2021).
We used three lines of human genetics evidence based on disease variant
frequency to support the statistical and biological association between human
genetic variation in a drug target and their impact on medical outcomes (see
Table 1).

###  2.1 Rare genetic diseases

We used a list of curated genes with reasonably well-established causal link
with a disease. This class of genetics evidence is enriched for genes
associated with rare Mendelian diseases whereby very rare variants with large
effect and high penetrance are the causative genetic alteration underlying the
disease or clinical manifestation. We leveraged expert manual curation of
diagnostic-grade gene-disease relationships from ClinGen (Strande et al.,
2017) and Genomics England (Martin et al., 2019) to annotate target-indication
(gene-disease) pairs.

###  2.2 Gene-level burden association

Gene-level collapsing methods combine information from genetic variants found
at appreciable frequencies in the general population and assess the
statistical strength of association between the aggregate variation and health
outcomes. We collated a list of significant gene burden associations across
thousands of target-disease pairs from UK Biobank. Public-private partner
institutes analyzed whole-exome sequencing data from more than 400,000 UK
Biobank participants to identify genes with coding variants that are
collectively enriched in individuals with a selected set of medical outcomes
[(Backman et al., 2021), (Wang et al., 2021), (Karczewski et al., 2022)]. Any
target-indication pair that did not reach the empirical significance threshold
in the respective study was labeled as negative association rather than
missing.

###  2.3 Evidence from GWAS

Genome-wide association studies (GWAS) test the statistical association
between common genetic variants and diseases, and have identified many
disease-variant associations that recapitulate known disease biology as well
as nominate novel therapeutic hypotheses. While GWAS have become an
indispensable tool in drug discovery for novel drug target identification and
validation, prioritization of causal genes among dozens of candidate genes
with equally compelling biological explanation remains a significant
challenge.

We leveraged a recently developed prioritization model, ’locus-to-gene’ (L2G)
model, that integrates human genetics from GWAS Catalog and UK Biobank, known
target and disease biology and multi-omic datasets (Ghoussaini et al., 2020),
(Mountjoy et al., 2021). We used the scoring predictions from the L2G model
from Open Targets Platform to annotate drug targets that are predicted to be
the causal genetics factor for a disease or phenotype.

###  2.4 Clinical trial outcomes

For clinical trial outcome data, we followed a similar label annotation
procedure as employed in (Yao et al., 2019). Specifically, we labeled every
drug target - indication pair as ”approved” (positive label) if there is at
least one drug molecule that has been approved for the corresponding
indication. For the remaining target-indication pairs, a ”failure” status
(negative label) was assigned if there is at least one clinical trial for the
pair that was either terminated or suspended. Additionally, we leveraged the
data from Open Target Platform’s NLP-based classification of clinical trials
to annotate the reason for trial failure. Clinical trials that were inferred
to be unfavorable due to safety concerns or efficacy were also assigned
negative labels.

###  2.5 Disease ontology

To facilitate the integration of indication data from multiple sources, we
mapped the EFO (experimental factor ontology) IDs for each indication to MeSH
IDs using EBI-EMBL ontology cross-reference database (OxO) (Malone et al.,
2010).

Table 1: Description of the three lines of human genetics evidence used in
this analysis.
Evidence type |  Description
---|---
Rare disease |  List of curated genes with established causal link between gene and disease.
Gene burden |  Gene-based rare variant associations in UK Biobank using whole exome sequencing data.
GWAS |  Prioritization of causal genes at GWAS locus based on genetic and functional genomics features using locus-to-gene (L2G) model.
Combined evidence |  Integrating human genetics evidence from all three types of evidence.

##  3 Model description

Given a binary matrix of drug targets (genes) and clinical outcomes (success,
failure, unknown), our goal is to impute the unknown cells or missing entries
using inter-relationships among targets and indications. We considered four
models incorporating human genetics evidence either individually or
altogether. Specifically, we created rank-3 tensors with each mode referring
to drug targets, indications and human genetics evidence, respectively, and
used Bayesian probabilistic matrix factorization using MCMC (Salakhutdinov &
Mnih, 2008) to factorize the binary matrices as implemented in SMURFF, a
highly optimized framework for Bayesian tensor factorization (Vander Aa et
al., 2019).

For each tensor factorization, we built a model with 32 latent dimensions and
used burn-in of 500 samples for the Gibbs sampler. We collected 3500 samples
from the model, and kept every 350th sample and averaged the predictions from
these samples for the final prediction.

##  4 Results

We evaluated the predictive performance of each model using AUROC, and the
model with combined evidence across three lines of human genetics evidence
performed slightly better than the other models (see Table 2). NLP-based
classification of clinical trial stop reasons yielded a small conservative set
of negative outcomes resulting in significant class imbalance between clinical
success and failure. To address the class imbalance, we also computed F1
scores for each model. In particular, the discrepancy between AUROC and F1
scores for the gene burden model highlights the dramatic class imbalance for
this model. It is very likely that a non-trivial fraction of target-indication
pairs may reach statistical significance when larger sample sizes and more
refined definition of indications are considered. Alternatively, more nuanced
set of rare and common variants with overlapping burden signal can be
considered (Weiner et al., 2022).

Table 2: Classification accuracies for the models considered in this study. F1 score was calculated using a threshold of 0.5. Class imbalance shows the proportion of positive labels out of total labels for the respective model. Model/Evidence | AUROC | F1 score | Imbalance
---|---|---|---
Rare disease |  93.2 $\pm$ 0.3 |  96.6 $\pm$ 0.2 | 87.2%
Gene burden |  92.6$\pm$ 0.3 |  81 $\pm$ 0.6 | 2.5%
GWAS |  93.3$\pm$ 0.2 |  95.4 $\pm$ 0.2 | 39.4%
Combined |  94.5$\pm$ 0.2 |  98.1$\pm$ 0.1 | 29.3%

We corroborate the previous finding that target-indication pairs from Phase 3
are enriched for validated or de-risked drug targets and therefore have higher
probability of success. clinical trials at later stages are more likely to
succeed (Yao et al., 2019) (see Figure 1).

Figure 1: Bayesian tensor factorization model prediction scores from the best
performing model (’combined model’). Each target-indication pair was grouped
by the maximum clinical phase reached. Preclinical phase refers to research
compounds that have not made to Phase I clinical trials yet. P-values were
calculated using two-sided Mann-Whitney-Wilcoxon test with Bonferroni
correction.

Interestingly, we also find that registered trials at preclinical stage
(research compounds) appear to have better odds of success than trials that
have progressed (i.e., Phase 1 and Phase 2) suggesting that sponsors consider
validated drug-indication pairs increasingly more often in their early drug
discovery research programs.

##  5 Discussion

Here we used publicly available data on approved drug indications and human
genetics evidence for the corresponding drug targets to predict the outcome of
ongoing clinical trials. Notably, we used Bayesian tensor factorization with
additional information from different lines of human genetics evidence to
support these targets. Our results show that the model with combined evidence
(rare disease, gene burden, common disease) modestly improves the accuracy of
predicting clinically promising drug targets when compared to alternative
models with single line of evidence. While this finding is encouraging for the
increasing appreciation for the role of human genetics in drug target
discovery and validation efforts, substantial class imbalance due to necessity
of expert manual curation poses a significant challenge for model comparison
and establishing benchmarks for further method development.

While the approved drug indications may be considered as true positive labels,
there are numerous reasons why the outcome of a clinical trial was not
favorable. Even when a clinical trial meets its primary objectives, trial
sponsors may choose not to move forward with the trial due to business
reasons, rapidly changing standard of care or anticipation of difficulty to
enroll eligible patients. Here we relied on NLP classification for labeling
the clinical trial outcome data, however, it is very likely that text-based
classifications may not completely capture the nature of a particular trial
failure. There is significant need for better documentation and structure of
clinical trial data to improve the effectiveness of text-based classification
and semantic analysis.

Choosing the most appropriate strategy to integrate different lines of human
genetics evidence remains to be an active area of research in drug discovery
(e.g., (Stacey et al., 2019), (Dornbos et al., 2022)). Historically, Mendelian
genetics evidence has proven to be a convenient source of strong causal
evidence between drug target and indication where the evidence implicates a
single mutation or gene as the molecular cause underlying the disease.
However, as DNA sequencing has become increasingly inexpensive, genomics
studies in much larger populations and more complex medical conditions have
begun to yield human genetics evidence for gene-disease associations in the
form of hundreds of mutations with individually marginal effects on health
outcomes and complex traits. While there are significant advances in
statistical and computational approaches to combining these effects for
patient stratification (see (Khera et al., 2018)), integrating evidence from
genetics associations from rare and common diseases has proven to be
difficult.

Understanding the relevant importance of different sources of human genetics
evidence for predicting clinically promising drug targets will significantly
help develop safe and effective therapies. In the experimental setup presented
in this study, we built multiple models to consider each source of human
genetics evidence separately and combined to compare the relative predictive
performance of each model. Notably, the burden model performed the poorest. It
is conceivable that the poor predictive performance is largely due to high
class imbalance in this model as well as relatively few available labels.
Further research is necessary to probe whether this class of genes with burden
evidence biologically represent difficult drug targets for therapeutic
modulation (e.g., highly selective targeting) or empirical significance
thresholds for these genes are too conservative.

Integration of matrix (or tensor) factorization approaches with neural
networks has proven to be very useful for predicting gene expression from
highly structured data modalities such as genomic and epigenomic data
(Schreiber et al., 2020). In case that there are substantial non-linear
relationship between different sources of human genetics data, this approach
can be useful for predicting trial outcome data using informative latent
representation of genetics evidence.

## Data Availability

All the data used in this analysis are publicly available on Open Targets
Platform (Ochoa et al., 2021): https://platform.opentargets.org/downloads.
Data on human genetics evidence and clinical trial outcomes were downloaded
from the latest release of the platform (v22.06). Detailed information on each
data source is available at https://github.com/cx0/icml-human-genetics.

## Acknowledgements

We are grateful to the Open Targets team and public/private partner
institutions for their commitment to open data sharing.

## References

  * Backman et al. (2021) Backman, J. D., Li, A. H., Marcketta, A., Sun, D., and Center, R. G.  Exome sequencing and analysis of 454,787 uk biobank participants.  _Nature_ , 599:628–634, 2021.
  * Cook et al. (2014) Cook, D., Brown, D., Alexander, R., March, R., Morgan, P., Satterthwaite, G., and Pangalos, M. N.  Lessons learned from the fate of astrazeneca’s drug pipeline: a five-dimensional framework.  _Nature Reviews Drug Discovery_ , 13:419–431, 2014.
  * Dornbos et al. (2022) Dornbos, P., Singh, P., Jang, D.-K., Mahajan, A., Biddinger, S. B., Rotter, J. I., McCarthy, M. I., and Flannick, J.  Evaluating human genetic support for hypothesized metabolic disease genes.  _Cell Metabolism_ , 34(5):661–666, 2022.
  * Ghoussaini et al. (2020) Ghoussaini, M., Mountjoy, E., Carmona, M., Peat, G., Schmidt, E., Hercules, A., Fumis, L., Miranda, A., Carvalho-Silva, D., Buniello, A., Burdett, T., Hayhurst, J., Baker, J., Ferrer, J., Gonzalez-Uriarte, A., Jupp, S., Karim, M., Koscielny, G., Machlitt-Northen, S., Malangone, C., Pendlington, Z. M., Roncaglia, P., Suveges, D., Wright, D., Vrousgou, O., Papa, E., Parkinson, H., MacArthur, J. A. L., Todd, J., Barrett, J. C., Schwartzentruber, J., Hulcoop, D., Ochoa, D., McDonagh, E. M., and Dunham, I.  Open Targets Genetics: systematic identification of trait-associated genes using large-scale genetics and functional genomics.  _Nucleic Acids Research_ , 49:1311–1320, 2020.
  * Karczewski et al. (2022) Karczewski, K. J., Solomonson, M., Chao, K. R., Goodrich, J. K., Tiao, G., Lu, W., Riley-Gillis, B. M., Tsai, E. A., Kim, H. I., Zheng, X., Rahimov, F., Esmaeeli, S., Grundstad, A. J., Reppell, M., Waring, J., Jacob, H., Sexton, D., Bronson, P. G., Chen, X., Hu, X., Goldstein, J. I., King, D., Vittal, C., Poterba, T., Palmer, D. S., Churchhouse, C., Howrigan, D. P., Zhou, W., Watts, N. A., Nguyen, K., Nguyen, H., Mason, C., Farnham, C., Tolonen, C., Gauthier, L. D., Gupta, N., MacArthur, D. G., Rehm, H. L., Seed, C., Philippakis, A. A., Daly, M. J., Davis, J. W., Runz, H., Miller, M. R., and Neale, B. M.  Systematic single-variant and gene-based association testing of thousands of phenotypes in 426,370 uk biobank exomes.  _medRxiv_ , 2022.
  * Khera et al. (2018) Khera, A. V., Chaffin, M., Aragam, K. G., Haas, M. E., Roselli, C., Choi, S. H., Natarajan, P., Lander, E. S., Lubitz, S. A., Ellinor, P. T., and Kathiresan, S.  Genome-wide polygenic scores for common diseases identify individuals with risk equivalent to monogenic mutations.  _Nature Genetics_ , 50:1219–1224, 2018.
  * King et al. (2019) King, E., Davis, W., and Degner, J.  Are drug targets with genetic support twice as likely to be approved? revised estimates of the impact of genetic support for drug mechanisms on the probability of drug approval.  _PLOS Genetics_ , 15(12), 2019.
  * Malone et al. (2010) Malone, J., Holloway, E., Adamusiak, T., Kapushesky, M., Zheng, J., Kolesnikov, N., Zhukova, A., Brazma, A., and Parkinson, H.  Modeling sample variables with an experimental factor ontology.  _Bioinformatics_ , 26(8):1112–1118, 2010.
  * Martin et al. (2019) Martin, A. R., Williams, E., Foulger, R. E., Leigh, S., Daugherty, L. C., Niblock, O., Leong, I. U., Smith, K. R., Gerasimenko, O., Haraldsdottir, E., et al.  Panelapp crowdsources expert knowledge to establish consensus diagnostic gene panels.  _Nature genetics_ , 51(11):1560–1565, 2019.
  * Mountjoy et al. (2021) Mountjoy, E., Schmidt, E. M., Carmona, M., Schwartzentruber, J., Peat, G., Miranda, A., Fumis, L., Hayhurst, J., Buniello, A., Karim, M. A., et al.  An open approach to systematically prioritize causal variants and genes at all published human gwas trait-associated loci.  _Nature Genetics_ , 53(11):1527–1533, 2021.
  * Nelson et al. (2015) Nelson, M. R., Tipney, H., Painter, J. L., Shen, J., Nicoletti, P., Shen, Y., Floratos, A., Sham, P. C., Li, M. J., Wang, J., Cardon, L. R., Whittaker, J. C., and Sanseau, P.  The support of human genetic evidence for approved drug indications.  _Nature Genetics_ , 47:856–860, 2015.
  * Ochoa et al. (2021) Ochoa, D., Hercules, A., Carmona, M., Suveges, D., Gonzalez-Uriarte, A., Malangone, C., Miranda, A., Fumis, L., Carvalho-Silva, D., Spitzer, M., Baker, J., Ferrer, J., Raies, A., Razuvayevskaya, O., Faulconbridge, A., Petsalaki, E., Mutowo, P., Machlitt-Northen, S., Peat, G., McAuley, E., Ong, C. K., Mountjoy, E., Ghoussaini, M., Pierleoni, A., Papa, E., Pignatelli, M., Koscielny, G., Karim, M., Schwartzentruber, J., Hulcoop, D. G., Dunham, I., and McDonagh, E. M.  Open targets platform: supporting systematic drug–target identification and prioritisation.  _Nucleic Acid Research_ , 49:1302–1310, 2021.
  * Ochoa et al. (2022) Ochoa, D., Karim, M., Ghoussaini, M., Hulcoop, D. G., McDonagh, E. M., and Dunham, I.  Human genetics evidence supports two-thirds of the 2021 fda-approved drugs.  _Nature reviews. Drug discovery_ , 2022.
  * Salakhutdinov & Mnih (2008) Salakhutdinov, R. and Mnih, A.  Bayesian probabilistic matrix factorization using markov chain monte carlo.  In _Proceedings of the 25th international conference on Machine learning_ , ICML ’08, pp. 880–887, 2008.  URL https://doi.org/10.1145/1390156.1390267.
  * Schreiber et al. (2020) Schreiber, J., Durham, T., Bilmes, J., and Noble, W. S.  Avocado: a multi-scale deep tensor factorization method learns a latent representation of the human epigenome.  _Genome biology_ , 21(1):1–18, 2020.
  * Stacey et al. (2019) Stacey, D., Fauman, E., Ziemek, D., Sun, B., Harshfield, E., and Wood, A.  Progem: a framework for the prioritization of candidate causal genes at molecular quantitative trait loci.  _Nucleic Acid Research_ , 47, 2019.
  * Strande et al. (2017) Strande, N. T., Riggs, E. R., Buchanan, A. H., Ceyhan-Birsoy, O., DiStefano, M., Dwight, S. S., Goldstein, J., Ghosh, R., Seifert, B. A., Sneddon, T. P., et al.  Evaluating the clinical validity of gene-disease associations: an evidence-based framework developed by the clinical genome resource.  _The American Journal of Human Genetics_ , 100(6):895–906, 2017.
  * Vander Aa et al. (2019) Vander Aa, Tom and Chakroun, I., Ashby, T. J., Simm, J., Arany, A., Moreau, Y., Le Van, T., Golib Dzib, J. F., Wegner, J., Chupakhin, V., Ceulemans, H., Wuyts, R., and Verachtert, W.  Smurff: a high-performance framework for matrix factorization, 2019.  URL https://arxiv.org/abs/1904.02514.
  * Wang et al. (2021) Wang, Q., Dhindsa, R., Carss, K., Harper, A., and Initiative, A. G.  Rare variant contribution to human disease in 281,104 uk biobank exomes.  _Nature_ , 597:527–532, 2021.
  * Weiner et al. (2022) Weiner, D. J., Nadig, A., Jagadeesh, K. A., Dey, K. K., Neale, B. M., Robinson, E. B., Karczewski, K. J., and O’Connor, L. J.  Polygenic architecture of rare coding variation across 400,000 exomes.  _medRxiv_ , 2022.  URL https://www.medrxiv.org/content/early/2022/07/07/2022.07.06.22277335.
  * Yao et al. (2019) Yao, J., Hurle, M. R., Nelson, M. R., and Agarwal, P.  Predicting clinically promising therapeutic hypotheses using tensor factorization.  _BMC Bioinformatics_ , 20(69), 2019.