# Coherent $J/\psi$ electroproduction on 4He and 3He at the Electron-Ion
Collider:
probing nuclear shadowing one nucleon at a time

Vadim Guzey  National Research Center “Kurchatov Institute”, Petersburg
Nuclear Physics Institute (PNPI),Gatchina, 188300, Russia  Matteo Rinaldi
Dipartimento di Fisica e Geologia, Università degli Studi di Perugia and
Istituto Nazionale di Fisica Nucleare, Sezione di Perugia, via A. Pascoli, I -
06123 Perugia, Italy  Sergio Scopetta  Dipartimento di Fisica e Geologia,
Università degli Studi di Perugia and Istituto Nazionale di Fisica Nucleare,
Sezione di Perugia, via A. Pascoli, I - 06123 Perugia, Italy  Mark Strikman
Pennsylvania State University, University Park, PA, 16802, USA  Michele
Viviani  INFN-Pisa, 56127 Pisa, Italy

###### Abstract

While the phenomenon of gluon nuclear shadowing at small $x$ has been getting
confirmation in QCD analyses of various LHC measurements involving heavy
nuclei, it has not been possible so far to establish experimentally the number
of target nucleons responsible for nuclear shadowing in a given process. To
address this issue, we study coherent $J/\psi$ electroproduction on 4He and
3He in the kinematics of a future electron-ion collider (EIC). Using an
accurate treatment of the wave functions of 4He and 3He, we show that this
process has the power to disentangle the contributions of the interaction with
a specific number of nucleons, in particular, the contribution of two nucleons
at the momentum transfer $t\neq 0$. This challenging, first ever calculation
of coherent vector meson production on light nuclei, using realistic nuclear
wave functions, provides an additional test of models of nuclear shadowing and
a novel probe of the gluon distribution in impact parameter representation,
allowing for a 3D imaging of gluons in light nuclei. Besides, thanks to this
analysis, unique information on the real part of the corresponding scattering
amplitude could be accessed.

QCD in nuclei, gluon distribution, nuclear shadowing, light nuclei, electron-
ion collider

###### pacs:

13.60.Hb,14.20.Dh,27.10.+h

Introduction. Studies of the phenomenon of nuclear shadowing have a long
history Frankfurt:1988nt ; Arneodo:1992wf ; Piller:1999wx ; Armesto:2006ph ;
Frankfurt:2011cs . In quantum mechanics and in the eikonal limit, it is
manifested in the total hadron-nucleus cross section being smaller than the
sum of individual hadron-nucleon cross sections. In essence, this is due to
the simultaneous interactions of the projectile with $k\geq 2$ nucleons of the
nuclear target, leading to a reduction (shadowing) of the total nuclear cross
section. In this framework the interaction of the projectile with a nucleus is
described by a sum of diagrams corresponding to the potential interaction with
individual nucleons, giving rise to the Glauber model Glauber:1955qq ;
Glauber:1970jm . However, it was demonstrated by Mandelstam Mandelstam:1963cw
and Gribov Gribov:1968jf  that the contribution of eikonal diagrams in quantum
field theory models tends to zero at high energies because, qualitatively,
there is not enough time between interactions with two nucleons for the
projectile to transform back into itself. As a result, shadowing in the high
energy limit is determined by the totality of diffractive interactions of the
projectile in different configurations Gribov:1968jf .

More recently the issue of nuclear shadowing in hard processes with nuclei was
discussed in the context of modifications of nuclear parton distribution
functions (PDFs) and an eventual onset of the regime of high gluon densities
(saturation), which are relevant for the physics programs of the Large Hadron
Collider (LHC) Salgado:2011wc ; Citron:2018lsq  and of a future Electron-Ion
Collider (EIC) Accardi:2012qut . Specifically, combining the Gribov-Glauber
approach to nuclear shadowing with the QCD factorization theorems for
diffractive and inclusive leading-twist processes in deep inelastic scattering
(DIS), a large leading-twist (LT) gluon nuclear shadowing at small $x$ was
predicted Frankfurt:1998ym  ($x$ is the nucleus momentum fraction carried by
the gluons). It was later confirmed by analyses of coherent photoproduction of
charmonia in ultraperipheral collisions (UPCs) of heavy ions at the LHC, which
showed Guzey:2013xba ; Guzey:2013qza ; Guzey:2020ntc  that $R_{g}(x=6\times
10^{-4}-10^{-3},\mu^{2}\approx 3\ {\rm
GeV}^{2})=g_{A}(x,\mu^{2})/[Ag_{N}(x,\mu^{2})]\approx 0.6$ for lead nuclei
($g_{A}(x,\mu^{2})$ and $g_{N}(x,\mu^{2})$ refer to the gluon density in the
nucleus and the nucleon, respectively). Alternative calculations performed in
the eikonal dipole models, where nuclear shadowing is a higher twist effect,
lead to a somewhat weaker shadowing, see e.g. Lappi:2013am .

The large magnitude of nuclear shadowing means that the gluons in heavy nuclei
probed in this process likely belong to more than one nucleon of the target,
due to their overlap in the transverse plane, see e.g. Refs. Frankfurt:2011cs
; Guzey:2016qwo . However, it is rather difficult to discriminate between
different mechanisms of shadowing using scattering off heavy nuclei since the
shadowing correction in this case is a result of summing a sign-alternating
and slowly converging series. In addition, measurements of coherent scattering
with the momentum transfer squared $t\neq 0$ represent a challenging task
because of a steep $t$ dependence. As a consequence, it is impossible to
establish the exact number of target nucleons involved in the process, which
hinders an access to information on the gluon dynamics in nuclei.

In this paper, starting from the successful description of the effect of
nuclear shadowing for heavy nuclei, we propose an alternative, complementary
strategy of studying coherent production of $J/\psi$ in DIS off 4He and 3He
light nuclei, which should be feasible at the EIC. Using specific features of
their response functions, namely, the presence of a zero in the one-body form
factor at moderate $\sqrt{|t|}=0.7$ GeV/c for 4He and $\sqrt{|t|}=0.8$ GeV/c
for 3He, we argue that it is possible to separate the contributions to nuclear
shadowing coming from the interaction with two and three nucleons of the
nuclear target. Besides, the ions under investigation have no excited states,
so that it is easy to select coherent events. This is based on an old idea
proposed initially in Ref. Levin:1975pc . Indeed, since the differential
$p+^{4}{\rm He}\to p+^{4}{\rm He}$ cross section does not present a minimum at
$-t\simeq 0.6$ GeV2, where the 4He charge form factor has a minimum, it has to
be dominated by effects beyond the impulse approximation (IA), namely, by the
interaction with several nucleons leading to nuclear shadowing. Supplementing
this with accurate calculations of one-, two-, and three-body form factors
based on exact solutions of the Schrödinger equation with realistic
potentials, which have greatly improved in the last decade, we make
predictions for the $d\sigma_{\gamma^{\ast}+^{4}{\rm He}\to J/\psi+^{4}{\rm
He}}/dt$ and $d\sigma_{\gamma^{\ast}+^{3}{\rm He}\to J/\psi+^{3}{\rm He}}/dt$
differential cross sections in a broad range of $t$, including the region of
finite $t$, where the IA is largely suppressed and the cross sections are
unambiguously sensitive to the contributions of the interaction with exactly
two nucleons of the target (the contribution corresponding to the interaction
with three and four nucleons near the minimum and for lower $|t|$ are
numerically small).

Note that runs with polarized 3He beams are planned at an EIC in order to
study the neutron spin structure Accardi:2012qut .

In the case of electron-deuteron scattering, the IA induced by the quadrupole
form factor dominates up to large $|t|$. Hence, we do not discuss this
reaction here. However, in a long run, if polarized deuteron beams for an EIC
become available, experiments using such beams would provide an independent
measurement of the double scattering amplitude (the interaction with $k=2$
nucleons). At the same time, the strategy discussed in our paper is probably
the only one allowing one to measure nuclear shadowing in light nuclei at
colliders, since its effect for the total electron-nucleus cross sections is
at the level of few %, see, e.g. Ref. Frankfurt:2003jf . Its smallness can be
readily seen by examining our predictions in Figs. 3, 4, and 5, where the
difference between the IA and full results at $t=0$ is twice the shadowing
effect for the total cross section.

Multiple scattering formalism for coherent electroproduction of $J/\psi$ on
light nuclei. As explained in the Introduction, at high energies projectiles
interact coherently with all nucleons of the nuclear target. The contributions
to the nuclear scattering amplitude corresponding to the interaction with
$k=1,2,3,\dots$ nucleons of the target is shown in Fig. 1; they interfere
destructively leading to the suppression of the nuclear cross section (nuclear
shadowing) Glauber:1955qq . In the Gribov-Glauber approach to nuclear
shadowing Gribov:1968jf , the contribution of the interaction with $k=2$
nucleons is unambiguously given by the diffractive (elastic) cross section on
the nucleon.

Figure 1: Multiple scattering series for $\gamma^{\ast}+A\to J/\psi+A$
scattering amplitude. Graphs $a$, $b$, $c$ correspond to the interaction with
$k=1$, 2, and 3 target nucleons, respectively. The zigzag lines labeled
$I\\!\\!P$ denote diffractive exchanges; the solid blob in graph $c$ stands
for the interaction with cross section $\sigma_{3}(x)$.

At the same time, the contributions corresponding to $k\geq 3$ nucleons cannot
be in general expressed in terms of diffraction on the nucleon and, hence,
need to be modeled. A convenient way to do it is offered by the Good-Walker
formalism of eigenstates of the scattering operator Good:1960ba ;
Blaettel:1993ah , which allows one to characterize the interaction with $k$
nucleons by the $k$th moment $\langle\sigma^{k}\rangle=\int d\sigma
P_{h}(\sigma)\sigma^{k}$. Here, the distribution $P_{h}(\sigma)$ parametrizes
the hadronic structure of the virtual photon and gives the probability for the
photon to interact with a nucleon with the cross section $\sigma$.

In this approach, the differential cross section of electroproduction
(photoproduction) of $J/\psi$ mesons on a nucleus $A$ can be written in terms
of the $d\sigma_{\gamma^{\ast}N\rightarrow J/\psi N}/dt$ cross section on the
proton at $t=0$ and the slope of its $t$ dependence $B_{0}(x)$, in the
following form

$\displaystyle\frac{d\sigma_{\gamma^{\ast}A\rightarrow J/\psi A}}{dt}(t)$
$\displaystyle=$ $\displaystyle\frac{d\sigma_{\gamma^{\ast}N\rightarrow J/\psi
N}}{dt}(t=0)$ (1) $\displaystyle\times$
$\displaystyle\left|F_{1}(t)e^{(B_{0}(x)/2)t}+\sum_{k=2}^{A}F_{k}(t)\right|^{2}\,,$

where the first term is proportional to the impulse approximation (Born
contribution), $F_{1}(t)=A{{\Phi_{1}(q)}}$, and the second term gives the
contribution of the interaction with $2\leq k\leq A$ nucleons, see Ref.
Levin:1975pc ,

$\displaystyle F_{k}(t)$ $\displaystyle=$
$\displaystyle\left(-\frac{1}{8\pi^{2}}\right)^{k-1}\binom{A}{k}\frac{\langle\sigma^{k}\rangle}{\langle\sigma\rangle}\frac{(1-i\eta)^{k}}{1-i\eta_{0}}\
\int\prod_{l=1}^{k}d^{2}\vec{q}_{l}$ (2) $\displaystyle\times$
$\displaystyle\,{{f(q_{l})}}\,{{\Phi_{k}(\vec{q},\vec{q}_{l})}}\,\delta^{2}\left(\sum_{l=1}^{k}\vec{q}_{l}-\vec{q}\right)\,,$

where all the transferred momenta $\vec{q}_{i}$ have been taken to be purely
transverse; $t=-|\vec{q}|^{2}$. The nuclear structure is taken into account
via the $k$-body form factors $\Phi_{k}$

$\displaystyle{{\Phi_{k}(\vec{q}_{1},...\vec{q}_{k})}}$ $\displaystyle=$
$\displaystyle\int\prod_{i=1}^{A}\left\\{\frac{d\vec{p}_{i}}{(2\pi)^{3}}\right\\}$
(3) $\displaystyle\times$
$\displaystyle\psi^{\ast}_{P}(\vec{p}_{1}+\vec{q}_{1},...\vec{p}_{k}+\vec{q}_{k},...,\vec{p}_{A})\,$
$\displaystyle\times$
$\displaystyle\psi_{P}(\vec{p}_{1},...,\vec{p}_{k},...\vec{p}_{A})\,\delta\left(\sum_{i=1}^{A}\vec{p}_{i}\right)\,,$

which represent the probability amplitude for $k$ nucleons in the nucleus $A$
to interact with the probe absorbing the momenta
$\vec{q}_{1},\cdots,\vec{q}_{k}$ and then going back to the nucleus, which
recoils elastically with a total momentum transfer
$\vec{q}=\sum_{k}\vec{q}_{k}$. Besides, $\psi_{P}$ is the intrinsic wave
function of the nucleus $A$ with total momentum $P$. A pictorial
representation of the two-body form factor $\Phi_{2}$ is shown in Fig. 2. One
should notice that Eq. (3) represents an integral over $3^{A-1}$ variables of
complicated wave functions for each ${q_{i},...q_{k}}$ set of values. Thus, to
evaluate accurately Eq. (2), the function $\Phi_{k}$ is needed on an enough
dense grid of these $k$ variables. As a consequence, the computation time
grows dramatically with $A$ and $k$.

Figure 2: A sketch of the two-body form factor $\Phi_{2}$, see Eq. (3).

In Eq. (2), we used that the $t$ dependence and the energy dependence of the
$\gamma^{\ast}N\rightarrow J/\psi N$ amplitude entering the IA are in general
different from those of the soft hadronic fluctuation-nucleon amplitudes
responsible for the interaction with $k\geq 2$ nucleons. The latter are
characterized by the factor parametrizing its dependence on the momentum
transfer, $f(q_{l})$, and the ratio of the real to imaginary parts,
$\eta\neq\eta_{0}$, where $\eta_{0}$ is the ratio of the real to imaginary
parts of the $\gamma^{\ast}p\to J/\psi\,p$ amplitude.

Note that the multiple scattering formalism presented in this section provides
a good description of the data on proton-4He elastic scattering Bujak:1978rla
.

Results for 4He and 3He. Below, we specify input for our calculations and show
our predictions for the differential cross sections of $J/\psi$
electroproduction on 4He and 3He in generic kinematics of an EIC. We start
with 4He, for which preliminary results have been presented in the EIC Yellow
Report AbdulKhalek:2021gbh .

The $k$-body form factors $\Phi_{1}$, $\Phi_{2}$, and $\Phi_{3}$ have been
calculated using a realistic nuclear wave function evaluated along the lines
of Marcucci:2019hml  and using the N4LO chiral potential of Ref. Entem:2017gor
with a cutoff of 500 MeV and three-body forces. Within this approach, the
binding energy of 4He turns out to be $B(^{4}{\rm He})=-28.15$ MeV, which is
very close to the experimental value. A consistency check of the numerics is
provided by the relation
$\Phi_{2}(\vec{q}_{2}=0,\vec{q}_{1})=\Phi_{1}(\vec{q}_{1})$, which has been
numerically successfully tested.

A cumbersome realistic calculation of the $\Phi_{4}$ form factor has not been
performed since a very small contribution is predicted. This expectation is
supported by an estimate carried out within a harmonic oscillator shell model
and treating properly the center-of-mass motion, as suggested in Ref.
Levin:1975pc . We verified that this procedure is able to yield reasonably
well the relative sizes of the $\Phi_{2}$ and $\Phi_{3}$ contributions
obtained within the realistic analysis.

In the leading logarithmic approximation of perturbative QCD, coherent
electro- and photoproduction of $J/\psi$ on nuclei probes the gluon density of
the target Guzey:2013xba ; Guzey:2013qza , see Fig. 1. Therefore, it is
natural to assume that hadronic fluctuations of the photon described by the
distribution $P_{h}(\sigma)$ are similar to those for the nuclear gluon
distribution. In this case, one can express the ratio of moments
$\langle\sigma^{k}\rangle/\langle\sigma\rangle$ entering Eq. (2) in terms of
two effective cross sections $\sigma_{2}(x)$ and $\sigma_{3}(x)$
Frankfurt:2011cs ,

$\displaystyle\frac{\langle\sigma^{2}\rangle}{\langle\sigma\rangle}$
$\displaystyle=$
$\displaystyle\sigma_{2}(x)\,,\quad\frac{\langle\sigma^{3}\rangle}{\langle\sigma\rangle}=\sigma_{3}(x)\sigma_{2}(x)\,,$
$\displaystyle\frac{\langle\sigma^{4}\rangle}{\langle\sigma\rangle}$
$\displaystyle=$ $\displaystyle(\sigma_{3}(x))^{2}\sigma_{2}(x)\,,$ (4)

where we have explicitly indicated the dependence on Bjorken $x$ (the
invariant photon-nucleon energy). The $\sigma_{2}(x)$ cross section is
unambiguously determined by the probability of diffraction in DIS on the
proton in the gluon channel H1:2006zyl  and the slope of the $t$ dependence of
the $\gamma^{\ast}p\to Xp$ cross section ($X$ denotes the diffractively
produced final state), $B=6$ GeV-2 H1:2006uea . At $x=10^{-3}$, one finds
$\sigma_{2}(x)=25$ mb.

On the other hand, the $\sigma_{3}(x)$ cross section needs to be modeled
through $P_{h}(\sigma)$. Using two plausible models Frankfurt:2011cs , one
finds $\sigma_{3}(x)=30-50$ mb. However, one of the key advantages of light
nuclei compared to heavy nuclei is that the sensitivity to the value of
$\sigma_{3}$ is negligible since the contribution of the $k=3$ and $k=4$ terms
in Eq. (1) is very small in the studied range of $t$. In our numerical
analysis and in the plots, we used $\sigma_{3}(x)=35$ mb. Thus, the effect of
nuclear shadowing is determined by the interaction with two nucleons, whose
strength is controlled by $\sigma_{2}(x)$.

At the same time, away from the minimum of the cross section and for large
$|t|$, the three-body contribution becomes significant and reduces the cross
section by approximately a factor of $2$. Hence, accurate measurements at
large $|t|$ will allow one to extract $\sigma_{3}(x)$ as well.

Finally, for the slope of the $t$ dependence the $\gamma^{\ast}p\to J/\psi p$
cross section, we use $B_{0}(x)=4.5$ GeV-2 measured by H1 and ZEUS
collaborations at HERA (see Guzey:2013qza  for references). This value
corresponds to $x\simeq 10^{-3}$, which is typical for the EIC kinematics. In
addition, using the Gribov-Migdal relation, we estimate the $\eta_{0}$ and
$\eta$ ratios of the real and imaginary parts of the corresponding amplitudes
exploiting their measured energy dependence: $\eta_{0}=(\pi/2)\times 0.1\simeq
0.16$ and $\eta=(\pi/2)\times 0.2\simeq 0.3$. In our analysis, we neglected
the $t$ dependence of $\eta$ and $\eta_{0}$ since the slopes of the
corresponding scattering amplitudes weakly depend on energy (i.e., the slopes
of the Regge trajectories $\alpha^{\prime}(0)$ are small).

Figure 3 shows our predictions for the ratio of the differential cross section
for $J/\psi$ coherent production on 4He to that for the nucleon target at
$t=0$ as a function of $-t$ at $x=10^{-3}$.

Figure 3:  Ratio of the differential cross section for $J/\psi$ coherent
production on 4He to the same quantity for the nucleon target at $t=0$ as a
function of $-t$ at $x=10^{-3}$.

One can see from the figure that the cross section is dominated by the 1-body
(IA) and the 2-body rescattering dynamics. The first minimum is clearly
shifted from $-t=$ 0.45 GeV2 to $-t=$ 0.27 GeV2, essentially due to the 2-body
contribution. Since 1-body dynamics is under remarkable theoretical control,
it allows one to disentangle 2-body dynamics and unambiguously relate it to
leading-twist gluon nuclear shadowing. Note also that the clear minimum of the
$t$ dependence in the IA case is filled because $\eta_{0}\neq\eta\neq 0$ in
the full calculation. This represents a unique opportunity to measure the
ratios of the real to imaginary parts of the corresponding scattering
amplitudes.

The quality of the IA result can be tested at $x=0.05$, where it is expected
to be dominating in a broad range of $t$, due to a vanishingly small
contribution of the shadowing correction. This is illustrated in Fig. 4, which
presents the $x$ evolution of the gluon shadowing correction in 4He. It shows
the ratio of the differential cross section for $J/\psi$ coherent production
on 4He to the same quantity at $t=0$ as a function of $-t$. At $x=10^{-3}$,
the full result is shown. At $x=0.05$, the IA result is presented. In the
latter case, the parameters of the model of $J/\psi$ production have been
properly changed (in particular, we used $B_{0}(x)=3$ GeV-2 Frankfurt:2010ea
).

Figure 4:  Ratio of the differential cross section for $J/\psi$ coherent
production on 4He to the same quantity at $t=0$ as a function of $-t$: the IA
result at $x=0.05$ is compared to the full one at $x=10^{-3}$.

Note that this $x$ evolution of the $t$ dependence agrees with that predicted
in Ref. Frankfurt:2010ea , which was obtained considering HERA data; a check
of this model will be possible at the EIC. Since the 1-body contribution
dominates the cross section at $x=0.05$, where no shadowing is expected in a
wide range of $t$, the emergence of LT gluon shadowing at lower $x$ points to
a significant broadening in the impact parameter space of the nuclear gluon
distribution, as discussed in Ref. Guzey:2016qwo  for heavy nuclei. If
confirmed, this observation would be a relevant step towards a 3D imaging of
gluons in nuclei.

We have also repeated our analysis for the 3He system, which will be
systematically used at an EIC. In this case, the nuclear form factors
$\Phi_{1},\Phi_{2}$, and $\Phi_{3}$ have been calculated using a realistic
wave function developed along the lines of Ref. Marcucci:2019hml  and using
the AV18 nucleon-nucleon potential Wiringa:1994wb  and including UIX three-
body forces Pudliner:1995wk . Again, as a consistency check of the numerics,
the relation $\Phi_{2}(\vec{q}_{2}=0,\vec{q}_{1})=\Phi_{1}(\vec{q}_{1})$ has
been successfully tested.

Figure 5:  Ratio of the differential cross section for $J/\psi$ coherent
production on 3He to the same quantity for the nucleon target at $t=0$ as a
function of $-t$ at $x=10^{-3}$.

Figure 5 shows the ratio of the differential cross section for $J/\psi$
coherent production on 3He to the same quantity for the nucleon target at
$t=0$ as a function of $-t$ at $x=10^{-3}$. One can see from the figure that
the pattern of the $t$ dependence is similar to that found for 4He. Again, the
1-body and 2-body scatterings dominate the cross sections. A relevant shift in
the minimum momentum transfer due to 2-body dynamics is predicted and, since
the IA is under theoretical control, there are very good chances to
disentangle 2-body dynamics, i.e., LT gluon shadowing, from the 1-body
contribution. In addition, accurate measurements at large $|t|$ would also
allow one to extract the three-nucleon contribution.

One should also mention that within our approach, the average number of
participating nucleons (wounded nucleons) is
$\nu=A\sigma_{\gamma^{\ast}N}/\sigma_{\gamma^{\ast}A}$ Bertocchi:1976bq ,
which leads to $\nu\sim 1.7$ for $x=10^{-3}$. Thus, using the set of 3He and
4He data, one would be able to test the consistency of our description.

Conclusions. Measurements of coherent $J/\psi$ electroproduction at finite
values of the momentum transfer $t$ with light ion beams at an EIC can nicely
complement investigations performed at the LHC with ultraperipheral collisions
of heavy nuclei. In particular, it will be possible to establish how many
nucleons contribute to the impressive gluon shadowing effects seen at the LHC,
which constitutes important information hardly accessible in the LHC data
collected with heavy nuclei probing predominantly the $t\simeq 0$ values. We
demonstrated this by performing a realistic calculation for the 3He and 4He
systems at $t\neq 0$ and considering contributions coming from different
numbers of nucleons involved in the process. We have clearly shown that the
first diffraction minimum is shifted with respect to that predicted by the
impulse approximation calculation. Since the latter contribution is under good
theoretical control, very good opportunities to disentangle multinucleon
dynamics, in particular two-nucleon dynamics contributing to the gluon
shadowing phenomenon, are expected. It will also be possible to obtain unique
information on the real part of the corresponding scattering amplitude.
Analyzing the $x$ evolution of the $t$ dependence predicted in our
calculation, the emergence of LT gluon shadowing at low $x$ points to a
significant broadening of the gluon distribution in impact parameter space.
This is just an example of many possibilities offered by the process under
scrutiny towards a novel 3D imaging of gluons in nuclei.

The use of light ion beams would greatly expand the EIC potential for probing
the small-$x$ dynamics. We will perform further investigations considering
additional light ions (for example, deuteron beams), other vector mesons in
the final state allowing for a sizable longitudinal momentum transfer, deeply
virtual Compton scattering (DVCS) with light nuclei, and the $Q^{2}$
dependence of cross sections of these processes.

Acknowledgements. Supported in part by the STRONG-2020 project of the European
Union’s Horizon 2020 research and innovation programme: Fund no 824093. The
research of M.S was supported by the US Department of Energy Office of
Science, Office of Nuclear Physics under Award No. DE-FG02-93ER40771.

## References

  * (1) L. L. Frankfurt and M. I. Strikman, Phys. Rept. 160, 235-427 (1988).
  * (2) M. Arneodo, Phys. Rept. 240, 301-393 (1994).
  * (3) G. Piller and W. Weise, Phys. Rept. 330, 1-94 (2000) [arXiv:hep-ph/9908230 [hep-ph]].
  * (4) N. Armesto, J. Phys. G 32, R367-R394 (2006) [arXiv:hep-ph/0604108 [hep-ph]].
  * (5) L. Frankfurt, V. Guzey and M. Strikman, Phys. Rept. 512, 255-393 (2012) [arXiv:1106.2091 [hep-ph]].
  * (6) R. J. Glauber, Phys. Rev. 100, 242-248 (1955).
  * (7) R. J. Glauber and G. Matthiae, Nucl. Phys. B 21, 135-157 (1970).
  * (8) S. Mandelstam, Nuovo Cim. 30, 1148-1162 (1963).
  * (9) V. N. Gribov, Sov. Phys. JETP 29 (1969), 483-487 (Zh. Eksp. Teor.Fiz. 56 (1969) 892-901).
  * (10) C. A. Salgado, J. Alvarez-Muniz, F. Arleo, N. Armesto, M. Botje, M. Cacciari, J. Campbell, C. Carli, B. Cole and D. D’Enterria, et al. J. Phys. G 39, 015010 (2012) [arXiv:1105.3919 [hep-ph]].
  * (11) Z. Citron, A. Dainese, J. F. Grosse-Oetringhaus, J. M. Jowett, Y. J. Lee, U. A. Wiedemann, M. Winn, A. Andronic, F. Bellini and E. Bruna, et al. CERN Yellow Rep. Monogr. 7, 1159-1410 (2019) [arXiv:1812.06772 [hep-ph]].
  * (12) A. Accardi, J. L. Albacete, M. Anselmino, N. Armesto, E. C. Aschenauer, A. Bacchetta, D. Boer, W. K. Brooks, T. Burton and N. B. Chang, et al. Eur. Phys. J. A 52, no.9, 268 (2016) [arXiv:1212.1701 [nucl-ex]].
  * (13) L. Frankfurt, V. Guzey and M. Strikman, Phys. Rev. Lett. 91, 202001 (2003) [arXiv:hep-ph/0304149 [hep-ph]].
  * (14) L. Frankfurt and M. Strikman, Eur. Phys. J. A 5, 293-306 (1999) [arXiv:hep-ph/9812322 [hep-ph]].
  * (15) V. Guzey, E. Kryshen, M. Strikman and M. Zhalov, Phys. Lett. B 726, 290 (2013) [arXiv:1305.1724 [hep-ph]].
  * (16) V. Guzey and M. Zhalov, JHEP 1310, 207 (2013) [arXiv:1307.4526 [hep-ph]].
  * (17) V. Guzey, E. Kryshen, M. Strikman and M. Zhalov, Phys. Lett. B 816 (2021), 136202 [arXiv:2008.10891 [hep-ph]].
  * (18) T. Lappi and H. Mantysaari, Phys. Rev. C 87, no.3, 032201 (2013) [arXiv:1301.4095 [hep-ph]].
  * (19) V. Guzey, M. Strikman and M. Zhalov, Phys. Rev. C 95 (2017) no.2, 025204. [arXiv:1611.05471 [hep-ph]].
  * (20) E. M. Levin and M. I. Strikman, Sov. J. Nucl. Phys. 23 (1976), 216 LENINGRAD-75-203.
  * (21) M. L. Good and W. D. Walker, Phys. Rev. 120, 1857-1860 (1960).
  * (22) B. Blaettel, G. Baym, L. L. Frankfurt, H. Heiselberg and M. Strikman, Phys. Rev. D 47, 2761-2772 (1993).
  * (23) A. Bujak, A. Kuznetsov, B. Morozov, V. A. Nikitin, P. Nomokonov, Y. Pilipenko, V. Smirnov, E. Jenkins, E. Malamud and M. Miyajima, et al. Phys. Rev. D 23 (1981), 1895.
  * (24) R. Abdul Khalek, A. Accardi, J. Adam, D. Adamiak, W. Akers, M. Albaladejo, A. Al-bataineh, M. G. Alexeev, F. Ameli and P. Antonioli, et al. [arXiv:2103.05419 [physics.ins-det]].
  * (25) L. E. Marcucci, J. Dohet-Eraly, L. Girlanda, A. Gnech, A. Kievsky and M. Viviani, Front. in Phys. 8 (2020), 69 [arXiv:1912.09751 [nucl-th]].
  * (26) D. R. Entem, R. Machleidt and Y. Nosyk, Phys. Rev. C 96, no.2, 024004 (2017) [arXiv:1703.05454 [nucl-th]].
  * (27) A. Aktas et al. [H1], Eur. Phys. J. C 48, 715-748 (2006) [arXiv:hep-ex/0606004 [hep-ex]].
  * (28) A. Aktas et al. [H1], Eur. Phys. J. C 48, 749-766 (2006) [arXiv:hep-ex/0606003 [hep-ex]].
  * (29) L. Frankfurt, M. Strikman and C. Weiss, Phys. Rev. D 83, 054012 (2011) [arXiv:1009.2559 [hep-ph]].
  * (30) R. B. Wiringa, V. G. J. Stoks and R. Schiavilla, Phys. Rev. C 51, 38-51 (1995) [arXiv:nucl-th/9408016 [nucl-th]].
  * (31) B. S. Pudliner, V. R. Pandharipande, J. Carlson and R. B. Wiringa, Phys. Rev. Lett. 74, 4396-4399 (1995) [arXiv:nucl-th/9502031 [nucl-th]].
  * (32) L. Bertocchi and D. Treleani, J. Phys. G 3, 147 (1977)


11institutetext:  Department of Vascular Medicine, Academic Medical Center,
University of Amsterdam, 1105 AZ Amsterdam, the Netherlands22institutetext:
HORAIZON BV, Delft, The Netherlands 33institutetext: Department of Clinical
Epidemiology, Biostatistics and Bioinformatics,
University Medical Center Amsterdam, The Netherlands

# Interpretable Models via Pairwise permutations algorithm

Troy Maasland 2⋆2⋆ João Pereira Equal contribution to this work1122 Diogo
Bastos 1122 Marcus de Goffau 11 Max Nieuwdorp 11 Aeilko H. Zwinderman 33
Evgeni Levin 1122

###### Abstract

One of the most common pitfalls often found in high dimensional biological
data sets are correlations between the features. This may lead to statistical
and machine learning methodologies overvaluing or undervaluing these
correlated predictors, while the truly relevant ones are ignored. In this
paper, we will define a new method called pairwise permutation algorithm (PPA)
with the aim of mitigating the correlation bias in feature importance values.
Firstly, we provide a theoretical foundation, which builds upon previous work
on permutation importance. PPA is then applied to a toy data set, where we
demonstrate its ability to correct the correlation effect. We further test PPA
on a microbiome shotgun dataset, to show that the PPA is already able to
obtain biological relevant biomarkers.

###### Keywords:

permutation importance correlation PPA Diabetes.

##  1 Introduction

Measuring feature importance has often been plagued by high feature
correlations. One important drawback is the lack of a theoretical definition
for variable importance, in case variables are correlated [12] [14], even in
linear models [12]. From a clinical perspective, correlated biomarkers are of
high interest because they both may play a role in a shared biological pathway
identified by the model and yet exhibit different behaviour in other
circumstances. The method proposed in this paper, which will be referred to as
pairwise permutation algorithm (PPA), allows us to calculate the importance of
features without having to rely on the previously mentioned selection
approaches. Highly correlated features, which have a similar relation with the
output value, should have close importance ranks since they explain the same
variability in the data. The pairwise permutation algorithm aims to provide
feature importance values while avoiding the use of aggressive pre-selection
techniques, since these techniques might remove relevant information from the
data. It also manages to retain model interpretability by generating an
importance value per feature, even when applied to black box models. Moreover,
when working with highly dimensional biological data sets, it is simply not
feasible to try and address each of the correlations in the data individually.

##  2 Related Work

In this work, we focus on model-agnostic procedures which can be divided into
local and global methods. Local-based methods such as LIME and its variants
[6, 41] attempt to explain predictions on single data points by perturbing it
and building a simple, yet interpretable model on the perturbed predictions.
Similarly, SHAP [5], offers a local explanation based on the additional
prediction value each feature has when adding it to all possible feature
subsets. Unlike local-based methods, global methods are concerned with
determining the overall model behaviour and what features it values for its
prediction. For example, in clinical research, the goal is to determine
biomarkers that can identify a condition in the general population, or
potential targets for novel drug development. Therefore, in this setting, we
are mainly concerned with a more holistic view of feature importance i.e.
global. A notable example is that of permutation importance which was first
introduced by Breiman [2] in random forests as a way to understand the
interaction of variables that is providing the predictive accuracy. Suppose
that for a certain feature $i$ in data-set $\mathbb{X}$, we randomly permute
the instances’ values, and denote the resultant data-set by
$\mathbb{X}_{i}^{\pi}$. Permutation importance is defined as the difference in
the expected model loss on the original dataset and the original one:

$PI_{\\{i\\}}(f):=E\left[f(\mathbb{X}_{i}^{\pi})\right]-E\left[f(\mathbb{X})\right]$
(1)

For random forests, there is already available work that analyzes the
behaviour of this permutation importance, including the cases when high
correlations are present. Gregorutti et al [14] provided a theoretical
description of the effect of correlations on the permutation importance, a
phenomenon already observed by Toloși and Lengauer [14] [13]. Furthermore, a
feature selection procedure was introduced, which was more efficient in
selecting important, highly correlated variables[14]. Strobl et al showed that
the larger feature importance values for correlated predictors in random
forests were due to the preference for such predictors in the early splits of
the trees. A new conditional permutation-based feature importance calculation
was suggested, in order to circumvent this inflation, as well as the
depreciation for its correlated predictor [4]. Furthermore, Hooker and Mentch
proposed the ’permute and relearn’ approach [11]. Based on this approach we
define the relearned permutation importance as

$PI_{j}^{\pi L}=E\left[f^{\pi
j}(\mathbb{X_{t}})\right]-E\left[f(\mathbb{X_{t}})\right]$ (2)

In which $f^{\pi j}$ is the model trained on the train dataset $\mathbb{X^{\pi
j}}$, in which feature j is permuted, $f$ the model trained on the original
train dataset $\mathbb{X}$ and $\mathbb{X_{t}}$ the test dataset. One drawback
of this approach was also mentioned in the context of correlated features, as
this resulted in the compensation effect, in which the importance of the
correlated features was reduced [11]. Local based methods, such as the ones
introduced earlier, are focused on the contribution of each feature towards
individual predictions, whereas permutation importance gives us a more broad
estimation, since it is based on the overall accuracy of the model. While the
former approach provides a higher degree of interpretability, the latter is
usually more appropriate in a research environment, in which the aim would be
to discover new leads which could help researchers to investigate the
underlying biological mechanisms.

##  3 Pairwise Permutations Algorithm

###  3.1 Notation

We will refer to a single instance of the data-set as instance or point
interchangeably throughout the paper. We denote matrices, 1-dimensional arrays
and scalars with capital bold and regular text, respectively (e.g.
$\mathbf{X},\;\mathbf{x},\;\alpha$). Matrices’ columns and rows will be
denoted by $\mathbf{X}[:,i]$ and $\mathbf{X}[j,:]$, respectively. The expected
loss of a function given by:
$\frac{1}{N}\sum_{i=1}^{N}l\left[y,f(\mathbb{x}_{i})\right]$ will be denoted
by $E\left[f(\mathbb{X})\right]$.

###  3.2 Intuition

Features that are equally important for the output value should have similar
feature importance ranks, and these should not be affected by feature
correlation. In an attempt to prevent the compensation effect for correlated
features mentioned by Hooker and Mentch, we have chosen to permute all the
feature pairs and calculate the corresponding permutation importance of the
pair. A key assumption in our method is that the higher the correlations, the
larger should be the correction to that feature individual importance.

###  3.3 Definition

In this section, we define the pairwise permutation importance (PPI) as the
weighted average of the permutation importance values, computed using the
’permute and relearn’ approach defined in Equation 2. The correlations between
the feature pairs will act as the weights. Let $\mathbb{R}$ be the
correlations matrix between all the features and $\mathbb{R}_{i,j}$ the
correlation value between features i and j. Let $PI_{i,j}$ define the relearn
permutation importance (see Equation 2) when both the feature $i$ and $j$ have
been permuted together, and $PI_{i,i}$ the relearn permutation importance,
when only feature $i$ is permuted.

$PPI_{i}=\frac{1}{\underbrace{\sum\limits_{j=1}^{M}|\mathbb{R}_{i,j}|}_{q}}\underbrace{\left(PI_{i,i}+\sum\limits_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{M}|\mathbb{R}_{i,j}|\cdot PI_{i,j}\right)}_{p}$ (3)

Note that when a feature has no correlations in the data, according to the
previous equation, the PPI will actually follow the relearn permutation
importance, which in our terminology is the single permutation importance
(SPI). Since for complex data sets with thousands of features the
computational time can become infeasible ($\mathcal{O}(N^{2})$), one possible
simplification is to set a threshold and consider only the permutation pairs
with a correlation above it. We define this procedure in algorithm 1.

###  3.4 Expected Difference

It might be tempting to compute the expected loss of the model, perform the
permutation analysis and then compute the difference of the expected losses.
This is actually how Fisher et al. [3] defined the permutation importance.
However, we note that this procedure is sub-par as we show in the following
theorem:

###### Theorem 3.1

For a given function $f:\mathbb{R}^{M}\rightarrow\mathbb{R}$, let $\mathbb{X}$
and $\mathbb{x}$ be a sample and an instance from the domain of $f$,
respectively, $\mathbb{X}^{\epsilon}_{i}$ be $\mathbb{X}$ with permuted values
for the r.v. $X_{i}$ and $\tilde{\mathbb{x}}$ an instance from
$\mathbb{X}^{\epsilon}_{i}$. Then, for any loss function
$l\left[y,f(\mathbb{x})\right]$ and norm function
$||\cdot||:\mathbb{R}^{M}\rightarrow\mathbb{R}$ it holds that:
$\mathbb{E}\left[\lVert
l\left[y,f(\mathbb{x})\right]-l\left[y,f(\tilde{\mathbb{x}})\right]\rVert\right]\geq\lVert\mathbb{E}\left[l\left[y,f(\mathbb{x})\right]\right]-\mathbb{E}\left[l\left[y,f(\tilde{\mathbb{x}})\right]\right]\rVert$

###### Proof

Consider the following convex function $\varphi(x)=\lVert x\rVert$ for
$x=l\left[y,f(\mathbb{x})\right]-l\left[y,f(\tilde{\mathbb{x}})\right]$. Then,
by Jensen’s inequality:

$\displaystyle\mathbb{E}\left[\varphi(x)\right]\geq\varphi\left(\mathbb{E}\left[x\right]\right)\Leftrightarrow$
$\displaystyle\mathbb{E}\left[\lVert
l\left[y,f(\mathbb{x})\right]-l\left[y,f(\tilde{\mathbb{x}})\right]\rVert\right]\geq\lVert\mathbb{E}\left[l\left[y,f(\mathbb{x})\right]-l\left[y,f(\tilde{\mathbb{x}})\right]\right]\rVert$
$\displaystyle=\lVert\mathbb{E}\left[l\left[y,f(\mathbb{x})\right]\right]-\mathbb{E}\left[l\left[y,f(\tilde{\mathbb{x}})\right]\right]\rVert$

∎

This means that computing the expected value of the normed difference of
individual loss values is more robust to non-linear relationships between the
input variables then computing the difference of the normed expected loss
values.

Input: $\mathbb{X}$, $\mathbb{X_{t}}$, $\mathbb{y_{test}}$,
$E\left[f(\mathbb{X_{t}})\right]$, $\mathbb{R}$, $\alpha$

Output: $\mathbb{v}$

for  _feature $i$ in $\mathbb{X}$_ do

$p\leftarrow 0$, $q\leftarrow 0$ (equation 3);

for  _feature $j$ in $\mathbb{X}$_ do

if  _$|\mathbb{R}_{i,j}| >\alpha$_ then

Permute the feature pair ($i$, $j$) together in $\mathbb{X}$;

Retrain the model with the permuted input data $\mathbb{X_{i,j}^{\pi}}$;

Calculate the model’s error, $E\left[f^{\pi,i,j}(\mathbb{X_{t}})\right]$, on
the test data;

Calculate $PI_{i,j}$ through the relearn formula
$E\left[f^{\pi,i,j}(\mathbb{X_{t}})\right]-E\left[f(\mathbb{X_{t}})\right]$;

if  _i=j_ then

$p\leftarrow p+PI_{i,i}$;

else

$p\leftarrow p+|\mathbb{R}_{i,j}|\cdot PI_{i,j}$;

end if

$q\leftarrow q+|\mathbb{R}_{i,j}|$;

end if

end for

$PPI_{i}\leftarrow p/q$;

$\mathbb{v}$.append($PPI_{i}$);

end for

Algorithm 1 Pairwise permutations algorithm

##  4 Simulations with toy dataset

To see how our new PPA would behave for correlated features, we generated a
toy dataset, based on the one used by Hooker and Mentch [11]. The data was
created by assuming a linear regression model:

$y_{i}=x_{i1}+x_{i2}+x_{i3}+x_{i4}+x_{i5}+0x_{i6}+0.5x_{i7}+0.8x_{i8}+1.2x_{i9}+1.5x_{i10}.$
(4)

This was then turned into a classification model, by generating the binary
outcome y with the classification rule in Equation 5 with with
$\varepsilon\sim N(0,0.1)$.

$y_{i}=\begin{cases}1,&\text{for }y_{i}+\varepsilon_{i}\geq\overline{y}\\\
0,&\text{otherwise}\end{cases},$ (5)

All features were generated from a multivariate normal distribution
$N(0,\Sigma)$ with $\Sigma$ equal to the identity matrix, except that
$\Sigma_{12}=\Sigma_{21}=\rho=0.9$. All features were then transformed into a
uniform distribution, mimicking how Hooker and Mentch [11] generated their
data. In total, 1000 samples were generated.
In case the features are in the same scale, the coefficients in the linear
model can be seen as the conditional importance of the feature on all other
variables [4] [11]. Therefore, based on the magnitude of the coefficients, we
can rank the features on their importances, where features with the same
coefficients should be equally important, while a feature with a higher
coefficient should get higher importance than a feature with a lower
coefficient. The order of the features should not be affected by any
correlations between the features.
Using XGBoost with the logistic loss function as the classification algorithm
[17, 18], we performed $50$ stratified shuffle splits ($70\%$train/$30\%$test)
and measured the ROC AUC after adding a noise feature to the dataset and
standard scaling it. We found the optimal hyperparameters using a 5-fold cross
validation grid search. To compute the PPIs, a correlation threshold of 0.3
was used. Also, the SPIs were obtained.

###  4.1 Results

(a) $x_{1}$ and $x_{2}$ with $\rho=0.9$

(b) $x_{1}$ and $x_{2}$ with $\rho=0.9$

(c) $x_{1}$ and $x_{6}$ with $\rho=0.9$

Figure 1: Average rank $\pm$ standard error for each feature based on the
pairwise permutation importance algorithm for (a) and (c) and the single
permutation importance algorithm for (b).

The classification model obtained an average AUC of $0.97\pm 0.01$. As shown
by the average feature ranks in Figure 1(a), our new PPA is able to retrieve
the right order of feature importances, in which $x_{10}$ is clearly the most
important one, followed by $x_{9}$. As expected, $x_{6}$ and the random
variable are identified as the least important features. The results for the
SPI are shown in Figure 1(b). It is clearly shown that the PPA outperforms
this approach, as the SPI decreased the importance of the correlated features
$x_{1}$ and $x_{2}$ and was not able to retrieve the right order of feature
importances. This was also observed for the random forest algorithm by Hooker
and Mentch [11]. The toy dataset showed that in case two features have the
same coefficient in the linear model and are correlated, the PPA is able to
retrieve the right order for the feature importances. We also analysed the
effect of a correlation of $\rho=0.9$ between $x_{1}$ and $x_{6}$, by changing
the covariance matrix $\Sigma$ to $\Sigma_{16}=\Sigma_{61}=0.9$ and setting
the correlation between $x_{1}$ and $x_{2}$ to 0. This represents a case in
which an important feature is correlated to an irrelevant feature. However, we
saw in this case that the importance of $x_{1}$ was decreased by $x_{6}$,
while the importance of $x_{6}$ was increased by $x_{1}$, as shown in Figure
1(c). This could be expected as the grouped importance is shared equally
between both features, while in the case of features with different
importances, this might not be the right assumption. In this case, the PPA may
not be the appropriate choice

##  5 Microbial biomarkers for Type 2 Diabetes Mellitus

In this section we test the PPA on a real-world dataset, specifically
microbiome data. The goal is to obtain biologically relevant markers.
Therefore, we downloaded the Qin 2012 microbiome dataset from MLRepo [20],
[21]. This curated classification dataset contained shotgun data for 124
samples, representing Chinese healthy controls (n = 59) and Type 2 Diabetes
Mellitus (T2D) patients (n = 65). For full details of the preprocessing of the
raw sequence reads for datasets in MLRepo, see [20]. We used the same
procedure as in the previous section with some additional preprocessing.
First, the read counts were rarefied to 28 358 reads per sample, which was the
lowest observed number of reads in a sample. After that, features with less
than 6 reads per sample on average, representing a relative abundance of
0.02%, were removed. The final dataset consisted then of 124 samples with 377
microbial OTUs.

###  5.1 Results

The classification model was able to achieve an average ROC AUC score of 0.92
$\pm$ 0.05, as depicted in figure 2(a). Figure 2(b) represents the top 15 most
predictive OTUs in the classification model. Analyzing these OTUs (and several
more beyond the top 15) primarily highlights 2 main patterns. The strongest
pattern observed in the data, most likely represents an effect that T2D has on
the dietary behavior of these Chinese T2D patients. Lactobacillus acidophilus,
Acidaminococcus intestini and Anaerostipes caccae are strongly associated with
T2D and with each other in this dataset. A regular dose of L. acidophilus is
commonly recommended in Chinese Medicine [22]. Fermented soybean products are
popular in China (i.a.) and various of these products commonly contain L.
acidophilus [23], [24], [25]. Indeed, there is evidence that supports the
beneficial claims regarding these fermented products and T2D [26], [27].
Trans-aconitic acid in the urine is a biomarker for the consumption of soy
products [28] and A. intestini is known to be able to oxidise trans-aconitate
[29] converting it to acetate. A. caccae, is an acetate and lactate consuming
butyrate producer. Cross-feeding interactions between L. acidophilus and A.
caccae have been analyzed in detail in vitro [30]. Other butyrate producing
species, like Roseburia intestinalis, can have similar cross-feeding
interactions [31] but were not part of this specific pattern, but with the 2nd
main pattern (see below), suggesting that A. caccae was part of same fermented
soybean product popular with, or given to, these Chinese T2D patients that
likely also contained L. acidophilus and A. intestini.

(a) The ROC-AUC plot.

(b) Top 15 feature importance ranks.

Figure 2: (a) Individual ROC AUC curves for each shuffle and average ROC AUC
plot for all shuffles. (b) Average rank $\pm$ standard error for the top 15
ranked features based on the pairwise permutation importance algorithm.

The second pattern involves several butyrate producers (the Roseburia,
Faecalibacterium, Coprococcus genera, several Eubacterium species and
Anaerostipes hadrus) in a cross-feeding relationship with various acetate
producing dietary fibre degrading species (Blautia and Ruminococcus
representatives). This cluster of species is generally found to be negatively
associated with T2D, not just in this study throughout the diabetes microbiome
field [21, 32, 33, 34, 35]. Insufficient butyrate production has been
associated with both T1D and T2D development both in rats, mice and in humans
[36], [37], [38], [39]. Besides being used by colonocytes as a primary energy
source [40] butyrate is a powerful inhibitor of histone deacetylase, which has
emerged as a target in the control of insulin resistance [42], [43], [44].
Animal and in vitro studies have generally found a beneficial effect of
butyrate and acetate on glucose homeostasis and insulin sensitivity [45].

##  6 Conclusions

In this paper, we have set a first step in correcting the compensation effect,
observed for ’permute and relearn’ permutation importances in case correlated
features are present. Our new PPA is able to obtain the right ranking of
features, when two features are highly correlated and have the same
importance, stated by the magnitude of their coefficient, in linear models.
Furthermore, while not yet optimal for correlations between more than 2
features or correlated features with unequal importance related to the output
variable, our PPA is already able to obtain relevant biological insights in a
Chinese Type 2 Diabetes microbiome dataset.

##  7 Acknowledgments

We would like to thank Manon Balvers for helping with computational
experiments.

## References

  * [1] Kiers H, Smilde A. Stat. Meth. & Appl. (2007) 16: 193.
  * [2] Breiman L. Random forests. Machine Learning, 45(1):5–32, 2001.
  * [3] Fisher, Aaron & Rudin, Cynthia & Dominici, Francesca. (2018). Model Class Reliance: Variable Importance Measures for any Machine Learning Model Class, from the ”Rashomon” Perspective.
  * [4] Strobl C, Boulesteix AL, Kneib T, Augustin T, Zeileis A: Conditional variable importance for random forests. BMC Bioinformatics 2008, 9: 307. 10.1186/1471-2105-9-307
  * [5] Lundberg, Scott M, and Su-In Lee: A unified approach to interpreting model predictions. Advances in Neural Information Processing Systems. 2017.
  * [6] Ribeiro M, Singh S and Guestrin C: Why Should I Trust You?: Explaining the Predictions of Any Classifier. eprint arXiv:1602.04938. 2016.
  * [7] Pereira J, Groen AK, Stroes ESG, Levin E: Graph Space Embedding. Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence Main track. Pages 3253-3259. https://doi.org/10.24963/ijcai.2019/451. 2019
  * [8] Shrikumar A, Greenside P, and Kundaje A. Learning Important Features Through Propagating Activation Differences. arXiv preprint arXiv:1704.02685 (2017)
  * [9] Kootte R, Levin E, Salojärvi J, Smits L, Hartstra A, Udayappan S, Hermes G, Bouter K, Koopen A, Holst J, et al: Improvement of insulin sensitivity after lean donor feces in metabolic syndrome is driven by baseline intestinal microbiota composition. Cell Metab., 26 (2017), pp. 611-619
  * [10] Ojala M, Garriga GC: Permutation Tests for Studying Classifier Performance. Journal of Machine Learning Research 11 (2010) 1833-1863
  * [11] Hooker G, Mentch L, Please stop permuting features An explanation and alternatives, (2019), arXiv preprint arXiv:1905.03151v1
  * [12] Grömping U,(2009) Variable Importance Assessment in Regression: Linear Regression versus Random Forest, The American Statistician, 63:4, 308-319, DOI: 10.1198/tast.2009.08199
  * [13] Toloși L, Lengauer T, Classification with correlated features: Unreliability of feature ranking and solutions, Bioinformatics, (2011) 27(14): 1986-1994
  * [14] Gregorutti B, Michel B, Saint-Pierre P, Correlation and variable importance in random forests, Stat Comput, (2017) 27, 659-678
  * [15] Imangaliyev S, Keijser B, Crielaard W, Tsivtsivadze E, Personalized microbial network inference via co-regularized spectral clustering, Methods, (2015), 83:28-35
  * [16] Ishwaran H, et al, Variable importance in binary regression trees and forests, Electronic Journal of Statistics 1, (2007), 519-537
  * [17] Caruana R, Niculescu-Mizil A, Crew G, et al, Ensemble selection from libraries of models, Twenty-first international conference on Machine learning - ICML ’04. ACM, 2004:18
  * [18] Chen T, Guestrin C, XGBoost: A Scalable Tree Boosting System, Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. ACM, 2016:785–794
  * [19] Meinshausen N, Bühlmann P, Stability selection, J R Stat Soc Ser B (Statistical Methodol) (2010);72:417–473
  * [20] Vangay P, Hillmann B.M, Knights D, Microbiome Learning Repo (ML Repo): A public repository of microbiome regression and classification tasks, GigaScience (2019), 8:1-12
  * [21] Qin J, Li Y, Cai Z, et al., A metagenome-wide association study of gut microbiota in type 2 diabetes, Nature, (2012), 490:55-60
  * [22] Cohen M. R, The New Chinese Medicine Handbook: An Innovative Guide to Integrating Eastern Wisdom with Western Practice for Modern Healing, Fair Winds Press (2015)
  * [23] Chang S. Y, Kim D-H, Han M.J, Physicochemical and sensory characteristics of soy yogurt fermented with Bifidobacterium breve K-110, Streptococcus thermophilus 3781, or Lactobacillus acidophilus Q509011, Food Sci Biotechnol, (2010), 19:107-113
  * [24] Bedani R, Rossi EA, Isay Saad SM. Impact of inulin and okara on Lactobacillus acidophilus La-5 and Bifidobacterium animalis Bb-12 viability in a fermented soy product and probiotic survival under in vitro simulated gastrointestinal conditions. Food Microbiol. 2013;34(2):382‐389
  * [25] Kanda H, Wang H.L., Hesseltine C.W., et al., Yoghurt production by Lactobacillus fermentation of soybean milk, Process biochemistry, (1976), 11(4):23
  * [26] Kwon D.Y, Daily III J.W, Kim H.J, Antidiabetic effects of fermented soybean products on type 2 diabetes, Nutrition Research (2010), 30(1):1-13
  * [27] Mueller N.T, Odegaard A.O, Gross M. D, et al., Soy intake and risk of type 2 diabetes mellitus in Chinese Singaporeans, European Journal of Nutrition (2012), 51: 1033-1040
  * [28] Münger L.H, Trimigno A, Picone G, et al., Identification of urinary food intake biomarkers for milk, cheese, and soy-based drink by untargeted GC-MS and NMR in healthy humans, J Proetome Res. (2017), 16(9):3321-3335
  * [29] Cook G.M, Wells J.E, Russell J.B, Ability of Acidaminococcus Fermentans to oxidize trans-aconitate and decrease the accumulation of tricarballylate, a toxic end product of ruminal fermentation, Appl Environ Microbiol. (1994), 60(7):2533-7
  * [30] Moens F, Verce M, De Vuyst L, Lactate- and acetate-based cross-feeding interactions betweeen selected strains of Lactobacilli, Bifidobacteria and colon bacteria in the presence of inulin-type fructans, Int J Food Microbiol. (2017), 241:225-236
  * [31] Saulnier D.M.A, Spinler J.K, Gibson G.R, et al., Mechanisms of Probiosis and Prebiosis: Considerations for enhanced functional foods, Curr Opin Biotechnol. (2009), 20(2):135-41
  * [32] Goffau de M.C, Luopajärvi K, Knip M, et al., Fecal microbiota composition differs between children with beta-cell autoimmunity and those without, Diabetes (2013), 62(4):1238-44
  * [33] Hur K.Y, Lee M-S, Gut microbiota and metabolic disorders, Diabetes Metab J. (2015), 39(3): 198-203
  * [34] Hartstra A.V, Bouter K.E.C, Bäckhed F, et al., Insights into the role of the microbiome in obesity and type 2 diabetes, Diabetes Care (2015), 38(1):159-165
  * [35] Murri M, Leiva I, Gomez-Zumaquero J.M, et al., Gut microbiota in children with type 1 diabetes differs from that in healthy children: A case-control study, BMC Medicine (2013), 11(46)
  * [36] Noureldein M. H, Bitar S, Youssef N, Butyrate modulates Diabetes-linked gut dysbiosis: Epigenetic and mechanistic modifications, J Mol Endocrinol (2020), 64(1): 29-42
  * [37] Endesfelder D, Engel M, Davis-Richardson A.G, et al., Towards a functional hypothesis relating anti-islet cell autoimmunity to the dietary impact on microbial communities and butyrate production, Microbiome (2016) 4:17
  * [38] Jia L, Li D, Feng N, et al., Anti-diabetic effects of clostridium butyricum CGMCC0313.1 through promoting the growth of gut butyrate-producing bacteria in Type 2 Diabetic Mice, Sci Rep. (2017), 7(1):7046
  * [39] Khan S, Jena G, Sodium butyrate reduces insulin-resistance, fat accumulation and dyslipidemia in Type-2 Diabetic rat: A comparative study with metformin, Chem Biol Interact. (2016), 254:124-34
  * [40] Donohoe D.R, Garge N, Zhang X, et al., The microbiome and butyrate regulate energy metabolism and autophagy in the mammalian colon, Cell Metab. (2012), 13(5): 517-526
  * [41] Pereira, J., Groen, A. K., Stroes, E. S. G., Levin, E., Graph Space Embedding, Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence, IJCAI-19
  * [42] Sharma S, Taliyan R, Histone deacetylase inhibitors: Future therapeutics for insulin resistance and Type 2 Diabetes, Pharmacol Res. (2016),113(Pt A):320‐326
  * [43] Dirice E, Ng R.W.S, Martinez R, et al., Isoform-selective inhibitor of histone deacetylase 3 (HDAC3) limits pancreatic islet infiltration and protects female nonobese diabetic mice from diabetes, J Biol Chem. (2017), 292(43):17598-17608
  * [44] Khan S, Jena G, The role of butyrate, a histone deacetylase inhibitor in Diabetes Mellitus: Experimental evidence for therapeutic intervention, Epigenomics (2015), 7(4):669-80
  * [45] Canfora E.E, Jocken J.W, Blaak E.E, Short-chain fatty acids in control of body weight and insulin sensitivity, Nat Rev Endocrinol. (2015), 11(10):577-91


# Towards End-to-end Unsupervised Speech Recognition

Alexander H. Liu Wei-Ning Hsu $\ddagger$ Michael Auli $\ddagger$ Alexei
Baevski $\ddagger$

$\dagger$ MIT CSAIL $\ddagger$ Meta AI

<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
Work done during an internship at Meta AI.

###### Abstract

Unsupervised speech recognition has shown great potential to make Automatic
Speech Recognition (ASR) systems accessible to every language. However,
existing methods still heavily rely on hand-crafted pre-processing. Similar to
the trend of making supervised speech recognition end-to-end, we introduce
wav2vec-U 2.0 which does away with all audio-side pre-processing and improves
accuracy through better architecture. In addition, we introduce an auxiliary
self-supervised objective that ties model predictions back to the input.
Experiments show that wav2vec-U 2.0 improves unsupervised recognition results
across different languages while being conceptually simpler. The code is
available at
https://github.com/pytorch/fairseq/tree/master/examples/wav2vec/unsupervised

##  1 Introduction

State-of-the-art ASR models [1, 2] have greatly improved performance over the
recent past. These models have been relying heavily on human supervision -
requiring thousands of hours of transcribed speech. Consequently, the main
limiting factor in building ASR systems for many more languages is the need
for large amounts of labeled data, making it challenging or even impossible to
build such supervised ASR systems for the vast majority of languages of the
world [3].

Fortunately, recent advances [4, 5, 6] have shown that it is possible to train
ASR systems in an unsupervised manner where speech annotation is no longer
required. On a high level, these methods can be viewed as a two-step pipeline:
a pre-processing step followed by the unsupervised training step, with the
goal of transcribing speech into phoneme sequences.

The pre-processing step aims to transform input speech representations into
features that are better aligned to the underlying phonetic content. This is
done through speech segmentation and feature engineering. The goal of speech
segmentation is to discover boundaries in the input sequence to mark the
change from one phone to another. In practice, the boundaries can be
determined through clustering representations [6] or other unsupervised
methods [7, 8]. Besides speech segmentation, prior works also considered
feature engineering a necessary step. From the input speech, they reduced the
frequency by merging intra-segment features and the dimensionality through
k-means clustering [4] or principal component analysis (PCA) [6].

While the motivation of the pre-processing step is clear, the disadvantages
are also obvious. Prior work [9] discovered that the accuracy of segmentation
can dramatically impact subsequent steps, yet unsupervised speech segmentation
is usually far from perfect [10]. This results in cascading errors from
segmentation and other feature engineering steps. To address these issues, we
propose an end-to-end system, which, through a better architecture, is able to
replace all audio-side pre-processing steps with a single speech-to-phoneme
neural network.

In addition to relaxing the aforementioned limitation, we also improve the
unsupervised ASR training step. The core idea of this step is to train a model
which maps features extracted from a model pre-trained on unlabeled speech
data into phoneme sequences by restricting the output vocabulary distribution
and n-gram structures to be similar to real phoneme sequences. Existing
solutions for imposing such restrictions include adversarial training [4] and
empirical output distribution matching [5]. However, these methods do not
ensure that the output is closely related to the spoken content. To make up
for this deficiency, we propose a self-supervised auxiliary objective function
that provides content-based guidance to the model by utilizing pseudo labels
derived from the raw audio signal.

To summarize our work, we propose wav2vec-U 2.0 - an unsupervised ASR
framework with an end-to-end style training and an improved objective
function. In our experiments, we validate the effectiveness of the pipeline
simplification and the proposed modification of the objective function. We
show that wav2vec-U 2.0 performs better than the existing unsupervised
frameworks across different languages from high-resource to low-resource. We
also combine wav2vec-U 2.0 with self-training strategies [11, 6] to compare
against prior works and supervised systems. The implementation is based on
fairseq [12] and the code will be released soon.

##  2 Background

###  2.1 Unsupervised ASR with Adversarial Training

Liu et al. [4] first achieved unsupervised phoneme recognition on TIMIT [13],
a small and clean English dataset for proof-of-concept experiments. Their
system takes MFCC representations of speech as the input and performs the pre-
processing step combining speech segmentation with (a separately trained)
segmental audio word2vec [7] and feature discretization with k-means
clustering. Using the segmented discrete features, they train a phone-based
ASR system through adversarial training inspired by Generative Adversarial
Network (GAN) [14]. In the adversarial training framework, the ASR model
serves as the generator which transcribes the discrete index sequence into
phone sequence with a learnable embedding table followed by softmax
activation. A 2-layer convolution-based discriminator is trained to
distinguish the output of the generator and real phone sequences. The goal of
the generator is to output phone sequences that cannot be distinguished by the
discriminator. By iteratively training the two modules, the generator can
learn to map the pre-processed feature sequences into phone sequences in an
unsupervised manner.

###  2.2 Wav2vec-U

Different from the prior work, Baevski et al. [6] proposed wav2vec-U (Fig. 2)
that first achieved unsupervised speech recognition on benchmark speech
recognition datasets in many different languages. Wav2vec-U takes speech
representations extracted from the wav2vec 2.0 model [15] as the input instead
of human-defined features. Such representations can be obtained by self-
supervised training on audio and have been proven to be much more effective
than human-defined features for ASR.

For the pre-processing step, wav2vec-U applies a sequence of feature
engineering on the speech representations: k-means clustering input frames,
reducing the representation dimension with PCA, merging adjacent PCA features
that share the same k-means cluster, and finally, mean pooling features from
adjacent pairs of timesteps to better match the length of the phone sequences.
To train the ASR model taking the pre-processed features as input, wav2vec-U
followed the adversarial training approach with a single layer CNN generator
and a two-layer discriminator.

Since the proposed wav2vec-U 2.0 model is closely related to wav2vec-U, we
detail the prior work in Section 3.1 and highlight our contributions as
follows: 1) we remove the need for the audio pre-processing step, and 2) we
introduce an auxiliary objective that improves the segment transcriber.

###  2.3 Self-training of Unsupervised ASR

Though training unsupervised ASR through adversarial learning is feasible,
results have shown that the performance of such models is usually far from
supervised models. To close the gap, Chen et al. [11] discovered that GAN
training can be harmonized with Hidden Markov Models (HMMs), iteratively
training GAN and HMM with pseudo labels provided by the previous step and this
can improve the system performance by over 30% after 3 iterations. Following
this method, wav2vec-U proposed a 3-stage pipeline, training GAN followed by
HMM, followed by fine-tuning pre-trained wav2vec2.0 models. For the HMM
training, the unlabeled speech used for adversarial training is pseudo-labeled
by the generator along with a phoneme-to-phoneme Weighted Finite-State
Transducers (WFST; Mohri 16, Mohri et al. 17) to incorporate a phone-based
language model for decoding. Next, the unlabeled speech is re-labeled with the
HMM using a phoneme-to-word WFST to incorporate a word-based language model
for decoding. Finally, the HMM-labeled data is used to fine-tune a pre-trained
wav2vec2.0 model with the target being letters, decoded with a letter-to-word
WFST. The self-training strategy significantly boosted the performance of the
unsupervised system, making it comparable to the state-of-the-art fully-
supervised system from 2019 [18]. In this work, we followed the same self-
training approach to compare wav2vec-U 2.0 against existing ASR systems.

##  3 Methodology

Figure 1: Wav2vec-U [6]. The input wav2vec2.0 feature is pre-processed before
feeding into the generator as described in Section 2.2. The generator is
optimized through adversarial training against the discriminator as described
in Section 3.1.

Figure 2: Proposed wav2vec-U 2.0. The generator takes raw wav2vec 2.0 feature
as input without pre-processing step as described in Section 3.2. In addition
to adversarial training, an auxiliary self-supervised objective is introduced
with pseudo label derived from the raw waveform as described in Section 3.3.

###  3.1 Adversarial Training with wav2vec-U

To introduce the proposed wav2vec-U 2.0, we begin with the adversarial
training defined in wav2vec-U (Fig. 2; Baevski et al. 6). Given $X$ the
sequence of pre-processed features (described in Section 2.2 and Fig. 2) and
$Y_{\text{u}}$ the phonemized unpaired text sequence, the goal of adversarial
training is to train a Generator $\mathcal{G}$ which transcribes $X$ into a
phone sequence that cannot be distinguished from $Y_{\text{u}}$ by a
discriminator $\mathcal{C}$. The full objective function is defined as

$\displaystyle\min_{\mathcal{G}}~{}\max_{\mathcal{C}}{}$
$\displaystyle\mathop{\mathbb{E}}_{Y_{\text{u}}}\left[\log\mathcal{C}(Y_{\text{u}})\right]+\mathop{\mathbb{E}}_{X}\left[\log\left(1-\mathcal{C}(\mathcal{G}(X))\right)\right]$
(1)
$\displaystyle-\lambda\mathcal{L}_{gp}+\gamma\mathcal{L}_{sp}+\eta\mathcal{L}_{pd}$

where the first two terms are the standard GAN objective [14]; the third term

$\displaystyle\mathcal{L}_{gp}$ $\displaystyle=$ (2)
$\displaystyle\mathop{\mathbb{E}}_{X,Y_{\text{u}},\alpha\sim
U(0,1)}\left[\left(\|\nabla\mathcal{C}(\alpha~{}\mathcal{G}(X)+(1-\alpha)~{}Y_{\text{u}})\|-1\right)^{2}\right]$

is the gradient penalty on discriminator calculated with respect to a sample
mixing real and fake input for better convergence [19] where $\alpha$ is the
mixing weight sampled from uniform distribution; the fourth term

$\mathcal{L}_{sp}=\sum_{(p_{t},p_{t+1})\in\mathcal{G}(X)}\|p_{t}-p_{t+1}\|^{2}$
(3)

is the smoothness penalty to penalize locally inconsistent phone prediction
from the generator where $p_{t}$ is the unnormalized generator output at time
$t$; the fifth term

$\mathcal{L}_{pd}=\frac{1}{|B|}\sum_{S\in B}-H_{\mathcal{G}}(\mathcal{G}(S))$
(4)

is the phoneme diversity loss encouraging a diversed output over the phone set
by maximizing entropy of the generator output distribution averaged over every
frame in a mini-batch $B$; and $\lambda,\gamma,\eta$ are the weights for each
term.

###  3.2 Removing Audio Pre-processing

As mentioned in Section 1, existing unsupervised ASR systems including
wav2vec-U rely on hand-crafted pre-processing. Interestingly, we discovered
that the pre-processing step in wav2vec-U is not necessary.

Without the need of speech segmentation, a generator taking raw speech
representation sequence as the input itself can be sufficient for unsupervised
ASR with proper output frequency (i.e. number of phones predicted per second).
While there are various ways to control the output frequency of the generator,
we found that simply changing the stride of the convolution neural network
works well. With the input wav2vec 2.0 features sampled at 50 Hz, we found
that outputting a phone candidate at around 16 Hz by using a stride of 3
worked consistently well across all 9 languages tested in this work. To
further downsample the 16 Hz generator outputs to results closer to phone
sequence, consecutive outputs sharing the identical most possible phone can be
treated as a single output by random sampling as done in wav2vec-U. In
addition, we found that PCA-based dimensionality reduction can be replaced by
a simple batch normalization [20] along the time axis over the wav2vec 2.0
features.

As illustrated in Fig. 2, wav2vec-U 2.0 simply takes the raw representations
extracted from wav2vec 2.0 model as input $X$ instead of pre-processed
features as is done in wav2vec-U (Fig. 2). These modifications allow wav2vec-U
2.0 to be trained without the need of a complicated pre-processing pipeline.
More importantly, learning the mapping from speech representations to phone
sequences in an end-to-end style avoids the defects brought by hand-crafted
pre-processing and performs better as demonstrated in Section 4.2.

###  3.3 Auxiliary Self-supervised Objective

Since speech transcriptions are not available in the unsupervised setting, the
generator relies only on adversarial training to learn the mapping between
speech representations and phone units. The generator can possibly learn to
satisfy the adversarial criterion without learning the correct mapping - for
example, by consistently producing the most common n-grams regardless of the
input speech.

To address this issue, we introduce an auxiliary objective function that aids
self-supervised learning by reconstructing a pseudo label sequence
$Z=(z_{1},z_{2},...,z_{T})$ derived from the input audio. We use pseudo-labels
generated by k-means clustering on Mel-frequency cepstral coefficient (MFCC)
features with 64 clusters. We found that this choice of pseudo labels works
consistently well across different languages, and we show the effect of
different choices of pseudo labels in Section 4.3. The auxiliary self-
supervised objective for the generator is defined as

$\mathcal{L}_{ss}=-\sum_{t}\log P_{\mathcal{G}}(z_{t}|X)$ (5)

where $P_{\mathcal{G}}(z_{t}|X)$, the probability of predicting the pseudo
label $z_{t}$ at time $t$, is obtained by an additional linear transformation
with the softmax function from the intermediate layer of the generator as
illustrated in Fig. 2. The complete objective function of wav2vec-U 2.0 can be
obtained by adding the auxiliary objective to Eq. 1 with a weight $\delta$:

$\displaystyle\min_{\mathcal{G}}~{}\max_{\mathcal{C}}{}$
$\displaystyle\mathop{\mathbb{E}}_{Y_{\text{u}}}\left[\log\mathcal{C}(Y_{\text{u}})\right]-\mathop{\mathbb{E}}_{X}\left[\log\left(1-\mathcal{C}(\mathcal{G}(X))\right)\right]$
(6)
$\displaystyle-\lambda\mathcal{L}_{gp}+\gamma\mathcal{L}_{sp}+\eta\mathcal{L}_{pd}+\delta\mathcal{L}_{ss}$

The benefit of introducing the auxiliary loss is two-fold: a) it provides a
content-based regularization to ensure the output of the generator is closely
related to the input speech and b) it guides the generator with a more
explicit signal by approximating the underlying phone sequence with pseudo
labels derived from speech.

##  4 Experiments

###  4.1 Setup

We follow Baevski et al. [6] and perform experiments on 9 different languages
from LibriSpeech [21], Multi-Lingual LibriSpeech [22], and CommonVoice [23].
Unpaired text data is phonemized with an off-the-shelf phonemizer for English
[24] and other languages†††https://github.com/bootphon/phonemizer with silence
token padded on both sides of each sentence and inserted at word boundaries
with a probability of 0.5. LJSpeech [25] is also considered for preliminary
experiments on end-to-end models outputting characters. We summarize each
dataset below and refer readers to prior work for more details on the dataset
pre-processing.

LibriSpeech & Libri-Light. 960 hours of English speech from LibriSpeech,
discarding the transcriptions, is used for adversarial training. In addition,
53.2k hours of speech from Libri-Light is used for self-training. For the
unpaired text used to train wav2vec-U 2.0 and language model, we used the
official text corpus provided by LibriSpeech excluding the transcriptions of
the speech data [26]. Evaluation is done on the standard dev/test sets.

Multi-Lingual LibriSpeech. For each of the 6 languages (Dutch, French, German,
Italian, Portuguese, and Spanish) selected from Multilingual LibriSpeech, 100
hours of speech is used for adversarial training and the official text corpus
is used as unpaired text. Evaluation is done on the standard test set.

CommonVoice. We follow Baevski et al. [6] to select 2 (Kyrgyz and Tatar) out
of 38 languages from CommonVoice using 1.8 and 4.6 hours of unpaired speech
respectively. For the unpaired text of both languages, we use a mixed text
corpus following Baevski et al. [6]. Evaluation is done on the dev/test set
following Rivière et al. [27].

LJSpeech. This is a single-speaker corpus containing 24 hours of English
speech. Following prior works [28, 29] on weakly-supervised speech
recognition, we randomly select 300 utterances for validation and testing each
and used the rest as unpaired speech for training. For the unpaired text, we
used the same data as the LibriSpeech setup.

Implementation Details For English speech from LibriSpeech, the input features
are extracted from the 15th layer of a pre-trained wav2vec2.0 Large model [15]
encoded from the audio of which silences have been removed using an
unsupervised voice activity detection tool [30]. For non-English speech,
features are extracted from the XLSR-53 [31] model, a multi-lingual version of
wav2vec 2.0.

As illustrated in Fig. 2, the generator is a 2-layer CNN with batch
normalization over the input features. We found that initializing the scaling
factor of batch normalization to 30 for wav2vec 2.0 Large and 35 for XLSR-53
features is critical for convergence. During adversarial training, when
obtaining the output probability distribution from the generator, consecutive
frames that share the same most probable phone are merged into one by random
selection.

The discriminator is a 2-layer CNN that uses the output of the generator or a
one-hot representation of the real phone sequence as input. A single scalar
indicating whether the probability of the input sequence is real or not is
obtained by taking the average over output sequence followed by the sigmoid
activation.

The generator and discriminator are trained for 100k steps in total
interleaving the updates of each (50k updates for each) with a fixed learning
rate of 5e-5 and 3e-4 respectively. Each batch contains 160 audio samples and
160 randomly selected text samples. With the choice of weights in Eq. 6 being
1.0 / 1.5, 1.5 / 2.5, 0 / 3, 0.3 / 0.5 for $\lambda,\gamma,\eta,\delta$
respectively, we run each model with 4 different random seeds and we observe
that over 80% of the models will converge during training. Models are selected
with the Unsupervised Cross-Validation Metric proposed by Baevski et al. [6].

Pykaldi [32] is used to decode the phone-based ASR outputs into words. While
the generator must be trained with a stride of 3, we found decoding using a
language model with a stride of 2 (higher output frequency) leads to a more
accurate output.

###  4.2 Comparison to wav2vec-U

Table 1: Interpolation from wav2vec-U (Fig. 2) to wav2vec-U 2.0 (Fig. 2).
Phone Error Rate (PER) computed with greedy decoding on LibriSpeech dev-other
set averaged over 8 runs. Freq. refers to the frequency of sequence, i.e.
number of tokens per second.

| Pre-processing | Generator configuration | Result
---|---|---|---
| Adjacent | Cluster | PCA | Batch | Linear | Auxiliary | Stride | Freq. | Average
| pooling | pooling | reduction | norm. | proj. | loss | (Hz) | PER
wav2vec-U | ✓ | ✓ | ✓ | - | - | - | 1 | 14 | 18.8 $\pm$ 0.9
$\left.{\begin{gathered}\hfill\\\ \hfill\\\ \hfill\\\ \hfill\\\ \hfill\\\ \hfill\\\ \end{gathered}}\right\downarrow{}$ step (i) | - | ✓ | ✓ | - | - | - | 1 | 28 | $>100$
step (ii) | - | ✓ | ✓ | - | - | - | 2 | 14 | 18.5 $\pm$ 0.6
step (iii) | - | - | ✓ | - | - | - | 2 | 25 | $>100$
step (iv) | - | - | ✓ | - | - | - | 3 | 16 | 19.0 $\pm$ 0.9
step (v) | - | - | - | - | - | - | 3 | 16 | $>100$
step (vi) | - | - | - | ✓ | - | - | 3 | 16 | 16.4 $\pm$ 0.7
step (vii) | - | - | - | ✓ | ✓ | - | 3 | 16 | 15.9 $\pm$ 1.1
wav2vec-U 2.0 | - | - | - | ✓ | ✓ | ✓ | 3 | 16 | 13.6 $\pm$ 0.9
input wav2vec 2.0 feature | 50 | -
ground truth phone sequence |  $\sim$10 | -

We first compare wav2vec-U 2.0 with the closely related prior work wav2vec-U
[6] to show the benefit of removing pre-processing steps and introducing the
auxiliary objective.

Table 1 shows our proposed changes to wav2vec-U (Fig. 2 vs. Fig.2) where
changes are applied gradually from step (i) to step (vii).

Note, steps (i)-(iv) use pooling mechanisms to control the output frequency of
the model. Removing these steps results in higher frequency and the model
fails to converge to a reasonable accuracy. However, accuracy can be recovered
by adjusting the output frequency by changing the stride of the generator
(step (ii) & (iv)). The fact that models with higher output frequency (step
(i) & (iii)) failed to converge shows that having output frequency close to
ground truth (around 10 Hz) is important.

We also observe a similar performance to the original framework after
completely removing all pre-processing (step (iv)) and conclude that pre-
processing speech representations in order to segment them into phone-like
units are not necessary to obtain a working unsupervised ASR system, which can
benefit from learning in an end-to-end fashion.

Comparing steps (iv) and (v), we see that some feature engineering in form of
PCA dimensionality reduction is indeed important for adversarial training to
work, but it can be replaced with a simple batch normalization layer as in
step (vi).

Finally, looking at steps (vii) and the final step (viii), we see that
introducing the auxiliary self-supervised objective improves PER
significantly.

###  4.3 Choice of Pseudo Labels for Self-supervision

Table 2: Different choice of pseudo label for wav2vec-U 2.0. Phone Error Rate (PER) computed with greedy decoding on LibriSpeech dev-other set averaged over 8 runs.  Pseudo Label | # of class | Avg. PER
---|---|---
None | - | 15.9 $\pm$ 1.1
wav2vec2.0 VQ indices2 | 320$\times$2 | 16.6$\pm$2.2
k-means clustering wav2vec2.0 features | 32 | 16.4$\pm$1.4
64 | 15.5$\pm$1.8
128 | 15.9$\pm$0.9
k-means clustering MFCC audio features | 50 | 15.2$\pm$0.9
64 | 13.6$\pm$0.9
100 | 14.8$\pm$1.3
128 | 16.8$\pm$1.7

Next, we investigate different options for pseudo labels that serve as an
approximation of the underlying phone sequence for the auxiliary loss. We
experiment on LibriSpeech but found a similar trend on Multi-Lingual
LibriSpeech and CommonVoice.

To acquire the pseudo label, we need to categorize input features according to
some approximation of the underlying phone sequence from speech. We first
consider the wav2vec 2.0 vector quantization (VQ) indices‡‡‡Indices from both
codebooks in wav2vec2.0 are used as pseudo label with multi-task prediction
which come together with the input feature extraction. We also consider
k-means clustering on wav2vec 2.0 features which are used to detect phone
boundary in wav2vec-U. Finally, we consider performing k-mean clustering on
the Mel-frequency cepstral coefficients (MFCC) audio features, which have been
shown effective for self-supervised speech representation learning [33].

Note that the choice of the pseudo labels is important for the proposed
objective. Table 2 shows the results. VQ indices perform less well than the
baseline method where no pseudo labels are used. MFCC clustering yields the
best results, leading to a reduction of up to 2.3% PER, which is a 14%
relative improvement over the baseline. The number of pseudo label classes can
significantly affect results: 64 k-means clusters perform best on both
wav2vec2.0 PCA features and the MFCC audio features. Since the pseudo labels
are extracted from audio, we believe a smaller label set might not be diverse
enough to cover the phoneme inventory of the language and a larger label set
will encode redundant information in the raw audio.

###  4.4 Results on LibriSpeech

Table 3:  Word Error Rate (WER) on LibriSpeech with different language models
(LM) on the standard LibriSpeech dev/test sets.

Model | Unlabeled | LM | dev | test
---|---|---|---|---
speech (hours) | clean | other | clean | other
Supervised learning w/ 960 hours of speech
DeepSpeech 2 [34] | - | 5-gram | - | - | 5.33 | 13.25
Fully Conv [35] | - | ConvLM | 3.08 | 9.94 | 3.26 | 10.47
TDNN+Kaldi [36] | - | 4-gram | 2.71 | 7.37 | 3.12 | 7.63
SpecAugment [18] | - | RNN | - | - | 2.5 | 5.8
ContextNet [2] | - | LSTM | 1.9 | 3.9 | 1.9 | 4.1
Conformer [1] | - | LSTM | 2.1 | 4.3 | 1.9 | 3.9
Semi-supervised learning w/ 960 hours of speech
Transf. + PL [26] | 54k | CLM+Transf. | 2.00 | 3.65 | 2.09 | 4.11
IPL [37] | 54k | 4-gram+Transf. | 1.85 | 3.26 | 2.10 | 4.01
NST [38] | 54k | LSTM | 1.6 | 3.4 | 1.7 | 3.4
wav2vec 2.0 [15] | 54k | Transf. | 1.6 | 3.0 | 1.8 | 3.3
wav2vec 2.0 + NST [39] | 54k | LSTM | 1.3 | 2.6 | 1.4 | 2.6
Unsupervised learning
wav2vec-U | 54k | 4-gram | 13.3 | 15.1 | 13.8 | 18.0
wav2vec-U 2.0 | 54k | 4-gram | 9.8 | 13.1 | 9.9 | 13.9
Unsupervised learning + Self-Training
wav2vec-U | 54k | 4-gram | 3.4 | 6.0 | 3.8 | 6.5
wav2vec-U 2.0 | 54k | 4-gram | 3.5 | 6.0 | 3.7 | 6.3

Table 4: Word Error Rate (WER) on the Multilingual Librispeech (MLS) for
German (de), Dutch (nl), French (fr), Spanish (es), Italian (it) and
Portuguese (pt).

Model | Labeled | LM | de | nl | fr | es | it | pt | Avg.
---|---|---|---|---|---|---|---|---|---
data used
Labeled training hours (full) | 2k | 1.6k | 1.1k | 918 | 247 | 161 |
Supervised learning |  |
Pratap et al. [22] | full | 5-gram | 6.49 | 12.02 | 5.58 | 6.07 | 10.54 | 19.49 | 10.0
Unsupervised learning |  |
wav2vec-U | 0h | 4-gram | 32.5 | 40.2 | 39.8 | 33.3 | 58.1 | 59.8 | 43.9
wav2vec-U 2.0 | 0h | 4-gram | 23.5 | 35.1 | 35.7 | 25.8 | 46.9 | 48.5 | 35.9
Unsupervised learning + self-training |  |  |
wav2vec-U | 0h | 4-gram | 11.8 | 21.4 | 14.7 | 11.3 | 26.3 | 26.3 | 18.6
wav2vec-U 2.0 | 0h | 4-gram | 11.5 | 17.6 | 12.8 | 10.9 | 18.6 | 20.6 | 15.3

In Table 3, we compare wav2vec-U 2.0 with supervised, semi-supervised, and
unsupervised prior work on the well-studied English benchmark dataset
LibriSpeech. Without self-training, wav2vec-U 2.0 performs reducing word error
rate (WER) by 2% to 4% compared to wav2vec-U. This shows that removing the
pre-processing steps and introducing the self-supervised auxiliary objective
can lead to improved performance on speech recognition. Next, we combine our
method with the 3-stage self-training pipeline described in Section 2.3.
Results show that WER can be significantly reduced with self-training,
reaching below 3.7% WER on the clean data and 6.3% on the noisy data.

On the other hand, we observe the gap between wav2vec-U and wav2vec-U 2.0 is
much smaller after applying the self-training strategy. This suggested that
unsupervised ASR on LibriSpeech might be bottle-necked by the self-training
pipeline. Improving the self-training approach can likely further improve the
performance of unsupervised ASR.

###  4.5 Results on Multi-lingual LibriSpeech

We have also tested our approach on a number of European languages from Multi-
lingual LibriSpeech and show the results in Table 4.

We observe that wav2vec-U 2.0 also outperforms wav2vec-U in this setup,
especially on the languages with higher error rate. Overall, the proposed
framework is able to improve the average WER on MLS by 8% absolute without
self-training. This further demonstrates the benefit of having an end-to-end
learning framework with self-supervision.

The results show that wav2vec-U 2.0 consistently improves over wav2vec-U with
self-training. Average WER improves to 15.3% compared to 18.6% for wav2vec-U.
This is a significant improvement since the supervised topline is 10% WER. On
Portuguese, the WER is only 1.2% worse than the supervised topline.

###  4.6 Results on CommonVoice

Table 5: Phone Error Rate (PER) on CommonVoice.

Model | tt | ky
---|---|---
Supervised learning
Fer et al. [40] | 42.5 | 38.7
m-CPC [27] | 42.0 | 41.2
XLSR-53 [31] | 5.1 | 6.1
Unsupervised learning
wav2vec-U | 25.7 | 24.1
wav2vec-U 2.0 | 22.5 | 24.3
Unsupervised learning + HMM training
wav2vec-U | 13.7 | 14.9
wav2vec-U 2.0 | 10.1 | 15.2

Table 6: Recognition Error Rate on LJSpeech. Unit error rate refers to the
greedy decoding phone error rate for phone-based models or character error
rate for letter-based model. Word error rate is obtained through phone-to-word
or letter-to-word decoding.

Model | Labeled | Phone-based | Unit | Word
---|---|---|---|---
data used | (lexicon required) | error rate | error rate
Semi-supervised learning
Ren et al. [28] | 20min | ✓ | 11.7 | -
10min | ✓ | 64.2 | -
Liu et al. [29] | 10min | ✓ | 35.2 | -
Unsupervised learning
wav2vec-U 2.0 | 0min | ✓ | 16.0 | 14.4
wav2vec-U 2.0 letter-based | 0min | ✗ | 34.6 | 64.0

Other than English and European languages, we experiment on two low-resource
Turkic languages in CommonVoice and the results are shown in Table 6. Since we
only measure the phone error rate, we only use the HMM stage of the self-
training pipeline.

Compared to the supervised topline using 17 hours of labeled speech, wav2vec-U
2.0 is able to achieve 10.1% / 15.2% PER on Tatar / Kyrgyz using just 4.6 /
1.8 hours of speech. This once again shows the value of developing
unsupervised ASR for low-resource languages.

###  4.7 End-to-end Attempt on LJSpeech

While wav2vec-U 2.0 removed the need for human-engineered pre-processing steps
and pipeline design, it is still not fully end-to-end. Like conventional
supervised ASRs and existing unsupervised ASRs, wav2vec-U 2.0 relied on a
lexicon which is an expert-defined word-to-phone mapping. Since the phoneme
inventory is defined by a set of distinctive sounds, learning a speech-to-
phoneme transformation is naturally easier for the model than a speech-to-
grapheme mapping or a speech-to-word mapping. However, collecting a lexicon
for every language can be difficult in practice as it requires expertise in
both linguistics and phonetics.

We, therefore, evaluate wav2vec-U 2.0 without relying on lexicon by replacing
all phones (for both adversarial training and language modeling) with letters
that can be directly inferred from the text data. Results on LJSpeech with the
end-to-end version of wav2vec-U 2.0 are presented in Table 6. We first show
that wav2vec-U 2.0, using no labeled data at all, is already outperforming
previous works on semi-supervised ASR [28, 29] which use small amounts of
human-labeled data.

When switching to an entirely letter-based system without a lexicon, the unit
error rate increases substantially. Such degradation is reasonable as the
former needs to learn to spell alongside the mapping between acoustic and
writing units. Unfortunately, we also discover that the word error rate of the
letter-based model is significantly higher than the phone-based model. These
results suggest that end-to-end unsupervised ASR could be possible but further
research and development are necessary.

##  5 Conclusion

An end-to-end approach for unsupervised ASR is key to increasing applicability
to low-resource languages. In this work, we move towards this goal by removing
the need for human-engineered pre-processing and by improving accuracy. While
wav2vec-U 2.0 is simpler and more accurate, additional work is required to
enable a fully end-to-end approach. Future work includes simplifying the self-
training pipeline and removing the need for a phonemizer. Possible solutions
for the latter are incorporating unsupervised lexicon discovery [41] or
training unsupervised ASR with graphemic text units such as letters directly
as for supervised systems.

## References

  * Gulati et al. [2020] Anmol Gulati, James Qin, Chung-Cheng Chiu, Niki Parmar, Yu Zhang, Jiahui Yu, Wei Han, Shibo Wang, Zhengdong Zhang, Yonghui Wu, and Ruoming Pang.  Conformer: Convolution-augmented transformer for speech recognition.  _Proc. of Interspeech_ , 2020.
  * Han et al. [2020] Wei Han, Zhengdong Zhang, Yu Zhang, Jiahui Yu, Chung-Cheng Chiu, James Qin, Anmol Gulati, Ruoming Pang, and Yonghui Wu.  Contextnet: Improving convolutional neural networks for automatic speech recognition with global context.  _Proc. of Interspeech_ , 2020.
  * Lewis et al. [2016] M. Paul Lewis, Gary F. Simon, and Charles D. Fennig.  Ethnologue: Languages of the world, nineteenth edition.  Online version: http://www.ethnologue.com, 2016.
  * Liu et al. [2018] Da-Rong Liu, Kuan-Yu Chen, Hung-Yi Lee, and Lin shan Lee.  Completely unsupervised phoneme recognition by adversarially learning mapping relationships from audio embeddings.  _Proc. of Interspeech_ , 2018.
  * Yeh et al. [2019] Chih-Kuan Yeh, Jianshu Chen, Chengzhu Yu, and Dong Yu.  Unsupervised speech recognition via segmental empirical output distribution matching.  In _Proc. of ICLR_ , 2019.
  * Baevski et al. [2021] Alexei Baevski, Wei-Ning Hsu, Alexis Conneau, and Michael Auli.  Unsupervised speech recognition.  In A. Beygelzimer, Y. Dauphin, P. Liang, and J. Wortman Vaughan, editors, _Advances in Neural Information Processing Systems_ , 2021.  URL https://openreview.net/forum?id=QmxFsofRvW9.
  * Wang et al. [2018] Yu-Hsuan Wang, Hung-yi Lee, and Lin-shan Lee.  Segmental audio word2vec: Representing utterances as sequences of vectors with applications in spoken term detection.  In _2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 6269–6273. IEEE, 2018.
  * Chen et al. [2019a] Yi-Chen Chen, Sung-Feng Huang, Hung-yi Lee, Yu-Hsuan Wang, and Chia-Hao Shen.  Audio word2vec: Sequence-to-sequence autoencoding for unsupervised learning of audio segmentation and representation.  _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , 27(9):1481–1493, 2019a.
  * Chung et al. [2018] Yu-An Chung, Wei-Hung Weng, Schrasing Tong, and James R. Glass.  Unsupervised cross-modal alignment of speech and text embedding spaces.  _Proc. of NIPS_ , 2018.
  * Kreuk et al. [2020] Felix Kreuk, Joseph Keshet, and Yossi Adi.  Self-supervised contrastive learning for unsupervised phoneme segmentation.  _Proc. of Interspeech_ , 2020.
  * Chen et al. [2019b] Kuan-Yu Chen, Che-Ping Tsai, Da-Rong Liu, Hung-Yi Lee, and Lin shan Lee.  Completely unsupervised speech recognition by a generative adversarial network harmonized with iteratively refined hidden markov models.  In _Proc. of Interspeech_ , 2019b.
  * Ott et al. [2019] Myle Ott, Sergey Edunov, Alexei Baevski, Angela Fan, Sam Gross, Nathan Ng, David Grangier, and Michael Auli.  fairseq: A fast, extensible toolkit for sequence modeling.  In _Proc. of NAACL System Demonstrations_ , 2019.
  * Garofolo et al. [1993] John S. Garofolo, Lori F. Lamel, William M. Fisher, Jonathon G. Fiscus, David S. Pallett, and Nancy L. Dahlgren.  The DARPA TIMIT Acoustic-Phonetic Continuous Speech Corpus CDROM.  _Linguistic Data Consortium_ , 1993.
  * Goodfellow et al. [2014] Ian J. Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial networks.  _Proc. of NIPS_ , 2014.
  * Baevski et al. [2020] Alexei Baevski, Yuhao Zhou, Abdelrahman Mohamed, and Michael Auli.  wav2vec 2.0: A framework for self-supervised learning of speech representations.  In _Proc. of NeurIPS_ , 2020.
  * Mohri [1997] Mehryar Mohri.  Finite-state transducers in language and speech processing.  _Computational linguistics_ , 23(2):269–311, 1997\.
  * Mohri et al. [2002] Mehryar Mohri, Fernando Pereira, and Michael Riley.  Weighted finite-state transducers in speech recognition.  _Computer Speech & Language_, 16(1):69–88, 2002\.
  * Park et al. [2019] Daniel S. Park, William Chan, Yu Zhang, Chung-Cheng Chiu, Barret Zoph, Ekin D. Cubuk, and Quoc V. Le.  Specaugment: A simple data augmentation method for automatic speech recognition.  In _Proc. of Interspeech_ , 2019.
  * Gulrajani et al. [2017] Ishaan Gulrajani, Faruk Ahmed, Martin Arjovsky, Vincent Dumoulin, and Aaron Courville.  Improved training of wasserstein gans.  _Proc. of NIPS_ , 2017.
  * Ioffe and Szegedy [2015] Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In _International conference on machine learning_ , pages 448–456. PMLR, 2015.
  * Panayotov et al. [2015] Vassil Panayotov, Guoguo Chen, Daniel Povey, and Sanjeev Khudanpur.  Librispeech: an asr corpus based on public domain audio books.  In _Proc. of ICASSP_ , pages 5206–5210. IEEE, 2015.
  * Pratap et al. [2020] Vineel Pratap, Qiantong Xu, Anuroop Sriram, Gabriel Synnaeve, and Ronan Collobert.  Mls: A large-scale multilingual dataset for speech research.  In _Proc. of Interspeech_ , 2020.
  * Ardila et al. [2020] Rosana Ardila, Megan Branson, Kelly Davis, Michael Henretty, Michael Kohler, Josh Meyer, Reuben Morais, Lindsay Saunders, Francis M. Tyers, and Gregor Weber.  Common voice: A massively-multilingual speech corpus.  _Proc. of LREC_ , 2020.
  * Park and Kim [2019] Kyubyong Park and Jongseok Kim.  g2pe.  https://github.com/Kyubyong/g2p, 2019.
  * Ito and Johnson [2017] Keith Ito and Linda Johnson.  The lj speech dataset.  https://keithito.com/LJ-Speech-Dataset/, 2017.
  * Synnaeve et al. [2020] Gabriel Synnaeve, Qiantong Xu, Jacob Kahn, Tatiana Likhomanenko, Edouard Grave, Vineel Pratap, Anuroop Sriram, Vitaliy Liptchinsky, and Ronan Collobert.  End-to-end ASR: from Supervised to Semi-Supervised Learning with Modern Architectures.  _Proc. of ICML workshop on Self-supervision in Audio and Speech (SAS)_ , 2020.
  * Rivière et al. [2020] Morgane Rivière, Armand Joulin, Pierre-Emmanuel Mazaré, and Emmanuel Dupoux.  Unsupervised pretraining transfers well across languages.  In _Proc. of ICASSP_ , 2020.
  * Ren et al. [2019] Yi Ren, Xu Tan, Tao Qin, Sheng Zhao, Zhou Zhao, and Tie-Yan Liu.  Almost unsupervised text to speech and automatic speech recognition.  In _International Conference on Machine Learning_ , pages 5410–5419. PMLR, 2019.
  * Liu et al. [2020] Alexander H Liu, Tao Tu, Hung-yi Lee, and Lin-shan Lee.  Towards unsupervised speech recognition and synthesis with quantized speech representation learning.  In _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 7259–7263. IEEE, 2020.
  * Tan et al. [2020] Zheng-Hua Tan, Achintya Kumar Sarkar, and Najim Dehak.  rvad: An unsupervised segment-based robust voice activity detection method.  _Computer speech & language_, 59:1–21, 2020.
  * Conneau et al. [2020] Alexis Conneau, Alexei Baevski, Ronan Collobert, Abdelrahman Mohamed, and Michael Auli.  Unsupervised cross-lingual representation learning for speech recognition.  _arXiv_ , abs/2006.13979, 2020.
  * Can et al. [2018] Dogan Can, Victor R. Martinez, Pavlos Papadopoulos, and Shrikanth S. Narayanan.  Pykaldi: A python wrapper for kaldi.  In _Proc. of ICASSP_ , 2018.
  * Hsu et al. [2021] Wei-Ning Hsu, Yao-Hung Hubert Tsai, Benjamin Bolte, Ruslan Salakhutdinov, and Abdelrahman Mohamed.  Hubert: How much can a bad teacher benefit ASR pre-training?  In _Proc. of ICASSP_ , 2021.
  * Amodei et al. [2016] Dario Amodei, Sundaram Ananthanarayanan, Rishita Anubhai, Jingliang Bai, Eric Battenberg, Carl Case, Jared Casper, Bryan Catanzaro, Qiang Cheng, Guoliang Chen, et al.  Deep speech 2: End-to-end speech recognition in english and mandarin.  In _Proc. of ICML_ , 2016.
  * Zeghidour et al. [2018] Neil Zeghidour, Qiantong Xu, Vitaliy Liptchinsky, Nicolas Usunier, Gabriel Synnaeve, and Ronan Collobert.  Fully convolutional speech recognition.  _arXiv_ , abs/1812.06864, 2018.
  * Xu et al. [2018] H. Xu, K. Li, Y. Wang, J. Wang, S. Kang, X. Chen, D. Povey, and S. Khudanpur.  Neural network language modeling with letter-based features and importance sampling.  In _Proc. of ICASSP_ , 2018.
  * Xu et al. [2020] Qiantong Xu, Tatiana Likhomanenko, Jacob Kahn, Awni Hannun, Gabriel Synnaeve, and Ronan Collobert.  Iterative pseudo-labeling for speech recognition.  _Proc. of Interspeech_ , 2020.
  * Park et al. [2020] Daniel S. Park, Yu Zhang, Ye Jia, Wei Han, Chung-Cheng Chiu, Bo Li, Yonghui Wu, and Quoc V. Le.  Improved noisy student training for automatic speech recognition.  _Proc. of Interspeech_ , 2020.
  * Zhang et al. [2020] Yu Zhang, James Qin, Daniel S. Park, Wei Han, Chung-Cheng Chiu, Ruoming Pang, Quoc V. Le, and Yonghui Wu.  Pushing the limits of semi-supervised learning for automatic speech recognition.  _Proc. of NeurIPS SAS Workshop_ , 2020.
  * Fer et al. [2017] Radek Fer, Pavel Matějka, František Grézl, Oldřich Plchot, Karel Veselỳ, and Jan Honza Černockỳ.  Multilingually trained bottleneck features in spoken language recognition.  _Computer Speech & Language_, 46, 2017.
  * Lee et al. [2015] Chia-ying Lee, Timothy J. O’Donnell, and James R. Glass.  Unsupervised lexicon discovery from acoustic input.  _TACL_ , 2015.