Also at ]School of Engineering and Information Technology, University of New
South Wales, Canberra, Australia

Also at ]Centre for Audio, Acoustics and Vibration, University of Technology
Sydney, Sydney, Australia

# Willis coupling-induced acoustic radiation force and torque reversal

Shahrokh Sepehrirahnama<EMAIL_ADDRESS>Centre for Audio,
Acoustics and Vibration, University of Technology Sydney, Sydney, Australia
Sebastian Oberst<EMAIL_ADDRESS>[  Centre for Audio, Acoustics
and Vibration, University of Technology Sydney, Sydney, Australia  Yan Kei
Chiang  [  School of Engineering and Information Technology, University of New
South Wales, Canberra, Australia  David A. Powell  School of Engineering and
Information Technology, University of New South Wales, Canberra, Australia

###### Abstract

Acoustic meta-atoms serve as the building blocks of metamaterials, with linear
properties designed to achieve functions such as beam steering, cloaking and
focusing. They have also been used to shape the characteristics of incident
acoustic fields, which led to the manipulation of acoustic radiation force and
torque for development of acoustic tweezers with improved spatial resolution.
However, acoustic radiation force and torque also depend on the shape of the
object, which strongly affects its scattering properties. We show that by
designing linear properties of an object using metamaterial concepts, the
nonlinear acoustic effects of radiation force and torque can be controlled.
Trapped objects are typically small compared to the wavelength, and are
described as particles, inducing monopole and dipole scattering. We extend
such models to a polarizability tensor including Willis coupling terms, as a
measure of asymmetry, capturing the significance of geometrical features. We
apply our model to a three-dimensional, sub-wavelength meta-atom with maximal
Willis coupling, demonstrating that the force and the torque can be reversed
relative to an equivalent symmetrical particle. By considering shape asymmetry
in the acoustic radiation force and torque, Gorkov’s fundamental theory of
acoustophoresis is thereby extended. Asymmetrical shapes influence the
acoustic fields by shifting the stable trapping location, highlighting a
potential for tunable, shape-dependent particle sorting.

Willis Coupling, Acoustic Radiation Force, Acoustic Radiation Torque,
Bianisotropy, metamaterials, Ultrasound

††preprint: APS/123-QED

Acoustic radiation force and radiation torque are the physical quantities
underlying acoustophoresis - the manipulation of sub-wavelength objects by an
incident acoustic field [1, 2, 3, 4, 5, 6, 7, 8]. Applications such as
ultrasonic sorting, separation and levitation have been developed by applying
force and torque to objects such as biological cells [7, 9, 10, 11, 12, 13,
14, 15]. Considering their magnitude relative to gravitational force and fluid
drag, acoustic radiation force and torque are most suitable for practical
manipulation of sub-wavelength objects, with sizes ranging from sub-
micrometers to a few millimeters [16, 17, 18, 19, 20, 21, 22, 23, 24, 8].
Existing theoretical derivations of the acoustic radiation force and radiation
torque use simple shapes such as spheres, spheroids, cylinders and disks,
leading to simple expressions which show the influence of the object’s volume,
material and aspect ratio [25, 26, 27, 28, 29, 30, 31]. This simplifying
assumption neglects asymmetry, which is a potential degree of freedom to tune
acoustophoretic processes, for example manipulation of a heterogeneous mixture
of sub-wavelength objects, as illustrated in Fig. 1. Neglecting shape
asymmetry, Gorkov’s theory of acoustophoresis shows that small objects can be
trapped at pressure or velocity nodes, cf. Fig. 1(a), depending on their
dominant scattering mode [20, 32]. However, when shape asymmetry is accounted
for, the trapping locations may shift away from the pressure node, as shown in
Fig. 1(b).

Figure 1: Acoustic trapping of objects with (a) zero and (b) non-zero Willis
Coupling, corresponding to symmetric and asymmetric shapes, respectively. The
asymmetry in the shape affects the locations of an acoustic trap and the
stable orientation of the objects in it. A meta-atom based on a Helmholtz
resonator exhibiting such non-conventional acoustic trapping. P.N. and V.N.
stand for Pressure Node and Velocity Node of an incident acoustic field.
Stable trap refers to the location of zero radiation force with negative force
gradient.

Considering the key role of acoustic radiation force and torque in the
ultrasonic range, it is important to capture the effects of geometric
asymmetry on these force and torque fields. Recent advances in asymmetric
metamaterials indicate growing interest in designs based on engineered
geometries, either as a material of a single meta-atom for beam steering [33,
34, 35, 36] or used as a material made of an array of meta-gratings for beam
splitting [37]. Such asymmetric shapes are characterised by Willis coupling,
which describes the monopole and dipole moments induced by the incident
velocity and pressure fields, respectively [38]. This coupling has a
theoretical upper limit and approaches zero as the object approaches a mirror-
symmetric configuration [33]. By geometrical modifications such as adding
protruding holes and internal cavities, the Willis coupling and its effect on
the acoustic field can be controlled [38, 39, 40, 36]. In particular, a two-
dimensional C-shaped Helmholtz resonator can achieve the maximum Willis
coupling near its resonance, or other arbitrary values by having multiple
apertures of different sizes [35].

To investigate the influence of Willis coupling on the acoustic radiation
force and torque, we consider a three-dimensional version of the Helmholtz
resonator based meta-atom, as illustrated in Fig. 1(b). The cavity can be
considered as a secondary scatterer that is employed to tune the overall
scattering response of the meta-atom. Since our design targets the ultrasonic
range, narrow and long passages are avoided due to their role in thermo-
viscous losses [36]. The losses for airborne acoustic meta-materials with
narrow passages depend on the depth of the viscous boundary layer, which is
inversely proportional to the square root of the frequency and is of the order
of micrometers [41]. Here, thermo-viscous losses are neglected for the
theoretical investigation of shape-controlled meta-material behavior. This
assumption is validated in the Supplementary Notes [42].

The 3D meta-atom has an internal cavity, occupying 60$\%$ of the volume, and
one aperture of width $W=0.2a$ with $a$ being the outer radius of the
spherical shell, which acts as the normalization parameter, cf. Fig. 1(b).
These parameters are tuned to achieve maximum Willis coupling along a single
axis [35]. The values of $a$ are varied from $0.01\lambda$ to $0.16\lambda$,
corresponding to $ka\approx 0.06$ to $0.96$, respectively, to investigate the
effects of shape asymmetry in the Rayleight limit of $ka\ll 1$. Beyond $ka>1$,
higher order terms, e.g. quadrupole, octopole etc., are required for a more
accurate representation of scattered pressure and the formulation of acoustic
radiation force and torque. We assume that the meta-atom behaves as a sound-
hard and immovable object in the acoustic domain. This allows comparison with
existing radiation force theory, which predicts, for a plane standing wave,
sound-hard immovable spheres or spheroids are always pushed towards the
pressure nodes. These are the stable acoustic trap locations with zero force
and negative force gradient, as illustrated in Fig. 1. For studying radiation
torque, a sound-hard immovable prolate spheroid of sub-wavelength size is used
as a reference [43, 44].  The radiation force and torque applied to
axisymmetric shapes were studied analytically and numerically [45, 46, 47, 48,
49] in different beam types [50, 51, 52, 53, 54]; however, the locations of
acoustic traps were the same as in Gorkov’s theory. Here, we demonstrate the
correlation between shape asymmetry and the location of an acoustic trap due
to the force reversal phenomenon.

The incident field is taken as a plane standing wave with
$p_{i}=P_{a}\cos(kz)e^{-j\omega t}$, and the meta-atom is represented by its
monopole-dipole polarizability tensor
$\pmb{\alpha}=\begin{bmatrix}\alpha_{pp}&\pmb{\alpha}_{pv}^{T}\\\
\pmb{\alpha}_{vp}&\pmb{\alpha}_{vv}\end{bmatrix}$. The acoustic radiation
force and torque are derived as the sum of symmetric and asymmetric terms
[55], $\mathbf{F}=\mathbf{F}_{\text{sym}}+\mathbf{F}_{\text{asym}}$ and
$\mathbf{T}=\mathbf{T}_{\text{sym}}+\mathbf{T}_{\text{asym}}$ with

$\displaystyle\begin{split}\mathbf{F}_{\text{sym}}=&\frac{kE_{i}}{\rho_{f}}\Big{(}\frac{\Re{\big{[}\alpha_{pp}\big{]}}}{\kappa_{f}}\mathbf{e}_{z}-kc_{f}\Im{\big{[}\pmb{\alpha}_{vv}\mathbf{e}_{z}\big{]}}\cdot\mathbf{e}_{z}\mathbf{e}_{z}\Big{)}\sin(2kz),\\\
\mathbf{F}_{\text{asym}}=&\frac{kE_{i}}{\rho_{f}}c_{f}\Im{\big{[}2\pmb{\alpha}_{pv}\cdot\big{(}\mathbf{e}_{z}\mathbf{e}_{z}\big{)}\big{]}}\cos(2kz),\end{split}$
(1)
$\displaystyle\begin{split}\mathbf{T}_{\text{sym}}=&-\frac{kE_{i}}{\rho_{f}}c_{f}\Im{\big{[}2\pmb{\alpha}_{vv}\mathbf{e}_{z}\big{]}}\times\mathbf{e}_{z}\sin^{2}(kz),\\\
\mathbf{T}_{\text{asym}}=&\frac{kE_{i}}{\rho_{f}}\Big{(}\frac{1}{\kappa_{f}}\Re{\big{[}\pmb{\alpha}_{vp}\big{]}}\times\mathbf{e}_{z}\Big{)}\sin(2kz),\end{split}$
(2)

where subscript sym denotes the partial force and torque associated with
direct-polarization coefficients $\alpha_{pp}$ and $\pmb{\alpha}_{vv}$ [2,
55]. The terms with the subscript asym are induced by the Willis coupling
coefficients $\pmb{\alpha}_{pv}$ and $\pmb{\alpha}_{vp}$, which obey
$\pmb{\alpha}_{pv}=j\omega\rho_{f}\pmb{\alpha}_{vp}^{T}$ for reciprocal media.
These terms show the correlation between shape asymmetry and the two fields of
acoustic radiation force and torque. This is illustrated in Fig. 1(b), with
objects exhibiting strong Willis coupling being trapped in a different
location than their symmetric counterparts, cf. Fig. 1(a). Eqs. (1) and (2)
were derived using Gorkov’s far-field approach, from the radiated momentum and
the monopole-dipole approximation of the scattered pressure field [55], and
are applicable to objects with arbitrary shape in a plane standing wave.
Similar expressions for radiation force and torque for other types of incident
waves, e.g. Bessel, Gaussian and Vortex beams, can be derived from Eqs. (17)
and (24) in [55]. The polarizability coefficients are calculated numerically
using the novel approach in [55], which is based on the Boundary Element
Method and multipole translation and rotation theory. Details of these
numerical calculations and the verification study are provided in [55].
Similar to the torque expression for spheroids, Eq. (21) in [44],
$T_{\text{sym}}$ in Eq. (2), and Eq. (I.9) in [42], show the spatial
dependences of $\sin^{2}(kz)$ and $\sin(2\theta)$ that apply to objects of
arbitrary shape in a plane standing wave.

An incident pressure field of 10 mm wavelength in air is considered, with the
speed of sound $c_{f}=343.140$ $\text{ms}^{-1}$ and mean density
$\rho_{f}=1.204$ $\text{kgm}^{-3}$, to investigate force reversal in airborne
levitation within the long wavelength range $ka<1$.  The acoustic radiation
force and torque values are normalized as force contrast factor
$\mathbf{Q}=\mathbf{F}/\big{(}E_{i}\pi a^{2}\big{)}$ and torque contrast
factor $\mathbf{Z}=\mathbf{T}/\big{(}E_{i}\pi a^{3}\big{)}$. The energy
density of the incident wave is $E_{i}=p_{i}^{2}/(4\rho_{f}c_{f}^{2})$.
Although the effects of viscous losses can be captured by the polarizability
coefficients, they are neglected since the depth of the viscous boundary layer
is at least three times smaller than the smallest $W$ in this study (for
validation see the Supplementary Notes [42]).

Figure 2: Panel (a) shows the absolute value of the non-zero and normalized
polarizability coefficients $\pmb{\alpha}$, for the 3D meta-atom. Normalized
radiation force and torque with respect to the size index $ka$ for the
$\theta=-\pi/4$ incidence angle are shown in panels (b) and (c), respectively.
The magnitude of these quantities are denoted by $|\mathbf{Q}|$ and
$|\mathbf{Z}|$. Panels (d) and (e) are the normalized force and torque,
respectively, for the case of $ka_{\text{(II)}}=0.157$ with respect to
position $z/\lambda$ and incidence angle $\theta$. The arrows indicate the
directions of the force and torque in different regions of panels (d) and (e).
Panel (f) shows the loci of zero force and zero torque and locations of saddle
points, unstable and stable acoustic traps. Numerical force and torque results
are obtained using direct integration of stresses over a fictitious surface
enclosing the object.

The non-zero coefficients of the normalized polarizability tensor of the meta-
atom are shown in Fig. 2(a), with respect to the size index $ka$. The results
are scaled by $6\pi/(c_{f}^{2}k^{3})$, which is the maximum permissible value
of normalized Willis coupling coefficients [33]. Real and imaginary parts of
polarizability coefficients are shown in Fig. (2) of the Supplementary Notes
[42]. A resonant peak is found at $ka\approx 0.17$, and a trough, indicating a
zero root, is observed for $\alpha_{pp}$ at $ka\approx 0.12$, implying no
contribution from the incident pressure field to the monopole moment. The
resonant peak is only observed in the polarizability coefficients that are
associated with the aperture facing the $z$-direction. These two effects
originate from the internal cavity, changing the scattering strength of the
monopole mode.  The radiation force and torque obtained using Eqs. (1) and (2)
are shown with respect $ka$ in Fig. 2(b) and Fig.2(c), respectively. The meta-
atom is oriented at a fixed angle of $-\pi/4$ with respect to the positive
wave direction $\hat{z}$, as shown in the inset of Fig.2(c), to ensure non-
zero torque in the $y$-direction. The force has two equal components in $x$\-
and $z$-directions. Equations (1) and (2) are verified by comparing with the
force obtained by direct integration of radiation stresses over a fictitious
surface of radius $4\lambda$ enclosing the meta-atom [2, 27, 55]. The
contributions of the direct and Willis-coupling partial forces are plotted in
Fig. 2(b). The local minimum in $Q_{\text{sym}}$ at $ka\approx 0.16$ emerges
due to the root of $\alpha_{pp}$ at $ka_{\text{root}}\approx 0.12$, and
$\alpha_{vv}^{zz}$ being non-zero within this entire range of $ka$. The
radiation torque exhibits no local minimum since it is independent of the
monopole scattering term $\alpha_{pp}$, as shown in Eq. (2). The orientation
of the object with respect to the wave propagation direction changes the
acoustic radiation force and radiation torque, as shown analytically in [42].

Figures 2(d) and 2(e) show the acoustic radiation force and radiation torque
with respect to position $z$ and incidence angle $\theta$ for the specific
case of $ka_{(\text{II})}\approx 0.16$, at which the radiation force is
dominated by the Willis-coupling partial force. The direction of the force and
torque are indicated by the arrows in each region. In contrast to the case of
symmetric objects, the radiation force is no longer zero at the pressure node
$z/\lambda=0.25$; implying a significant shift in the stable trapping location
due to shape asymmetry. It is also observed that the force depends on angle
$\theta$, meaning that the stable acoustic trap occurs where both the force
and torque are zero, corresponding to the intersection of the zero contours in
Fig. 2(d) and (e). These zero force and torque contours are plotted in Fig.
2(f), leading to four intersection points. Only one of these points is stable
(circle marker), since both the force and torque gradients are negative.
Analysis of the force and torque gradients reveals that the other three points
are unstable. For the point marked by a square, any perturbation of position
or orientation causes the particle to escape. The points marked with diamonds
are saddle points, where particular combinations of torque and force
perturbation will lead to the particle escaping.

Figure 3: Variations of the normalized radiation force, panels (a)-(c), and
radiation torque, panels (d)-(f), over a half wavelength span between two
velocity nodes of a plane standing wave in the $z$ direction. The results are
shown for  $ka_{\text{(I)}}=0.094$ (Rayleigh Regime), $ka_{\text{(II)}}=0.157$
(minimum $F_{sym}$), and $ka_{\text{(III)}}=0.176$ (resonance peak), and the
orientation angle $\theta=-\pi/4$. The equivalent spheroid is obtained from
the polarizability tensor of the meta-atom by setting the Willis couplings to
zero. For the equivalent sphere, $\alpha^{zz}_{vv}$ is additionally set to
$\alpha^{xx}_{vv}$. The trap locations of the equivalent sphere and spheroid
are found from the analytical expressions in [2, 56] and [57, 44],
respectively.

We now investigate the force and torque distributions at different positions
$z/\lambda$ within the standing wave, as shown in Fig. 3.  A sphere and a
prolate spheroid with similar polarizability response to the meta-atom are
considered as reference cases, with their acoustic trap being at the pressure
node [2, 56, 57, 44, 55]. We quantify the role of Willis coupling as how far
the trap location of the meta-atom shifts away from the pressure node.
Although equivalent objects are selected based on the polarizability response
of the meta-atom, the prediction of their trap location is independent of
their size as long as $ka\ll 1$. Three different frequency regimes are
considered to show the drastic impact of shape asymmetry on the trapping
location. Figures 3(a) and 3(d) show the results for $ka_{(\text{I})}=0.094$,
corresponding to a low frequency limit where the Willis coupling is much
smaller than the direct polarizability terms. As expected, the Willis-coupling
partial force (asym in Fig. 3(a)) is negligible compared to the direct
polarization contribution, $F_{\text{sym}}$ (sym in Fig. 3(a)). Since
$F_{\text{sym}}$ indicates the net force applied to the equivalent sphere, its
trapping point is at the pressure node, where the force is zero with a
negative gradient (Fig. 3(a)), matching the results of the classical theory
for sound-hard spheres in a plane standing wave [16, 18, 17, 2, 20, 32, 58].
However, the trap location of the meta-atom determined by the sum of the two
partial forces, shown by the solid blue line, is slightly shifted from the
pressure node, as a result of $F_{\text{asym}}$ being dominant in that region.
In Fig. 3(d), we observe the Willis-coupling torque being dominant, resulting
in an additional zero-torque location at $z/\lambda\approx 0.21$, unlike the
reference prolate spheroid that only experiences zero-torque at the velocity
nodes $z/\lambda=0$ and $z/\lambda=0.5$ [57, 44, 55].

In Fig. 3(b) and 3(e), we consider the frequency $ka_{(\text{II})}=0.157$,
corresponding to a local minimum of the direct force term $F_{\text{{sym}}}$,
so that the total radiation force is determined by the Willis-coupling
contribution. This causes the stable trap to shift from the pressure node by
$z/\lambda\approx 0.125$, to where the symmetric contribution to the force is
maximum - illustrating that asymmetric particles can show behavior directly
contradicting the conventional Gorkov theory. Although $F_{\text{{sym}}}$ is
relatively small, we observe a direction reversal of the force vector in
comparison to $ka_{(\text{I})}$. Such force reversal is typically expected
only for $ka>1$ [59, 28, 52, 54], due to the resonant modes of the outer fluid
domain. However, our proposed meta-atom undergoes such force reversal due to
the internal cavity acting as a Helmholtz resonator and the maximized Willis
coupling associated with its single aperture. The net torque is determined by
the large contribution from the Willis-coupling partial torque, and the stable
angular equilibrium occurs at a point further to the left of the pressure node
than for $ka_{(\text{I})}$.

Close to the resonance peak, at $ka_{(\text{III})}=0.176$ , the total
radiation force and $F_{\text{{sym}}}$ shown in Fig. 3(c) exhibit the same
direction as predicted by Gorkov theory. We also observe a reversal of the
Willis-coupling partial force $F_{\text{{asym}}}$, resulting in a large
negative shift of more than $10a$ from the pressure node. The reversal of the
torque shown in Fig. 3(f) indicates its tendency to orient the aperture in the
opposite direction to that obtained for $ka_{(\text{I})}$ and
$ka_{(\text{II})}$. The stable angular equilibrium is still found at
$z/\lambda\approx 0.21$, due to the direction reversal of both partial
torques.

Our results provide physical insight into the relationship between shape
asymmetry and force and torque reversal, which we show here are also possible
for $ka<1$. Willis coupling in general describes structure-inherent asymmetry,
which can be related to geometrical anomalies of a trapped object. In addition
to anomalous scattering, meta-atoms with strong Willis coupling were shown to
experience different acoustic radiation forces and radiation torques,
indicating the potential to control these nonlinear acoustic phenomena through
Willis coupling. Applications such as particle sorting, could be tuned to
process objects with similar shape features by generating traps at locations
other than pressure/velocity nodes. An example is to switch frequency during a
sorting process [60], between resonant and off-resonant states, which can be
used as a propulsion technique to drive micro/nano-bots between target
locations in an acoustofluidic system. Further, single objects can be re-
oriented, not only in air but also in other fluids, by adjusting the
wavelength according to its asymmetric feature size to maximise the Willis
coupling. As these effects already take place within the sub-wavelength
regime, target-specific ultrasound manipulation techniques for the control of
biological cells, rafts of particles and small tissue ensembles or mesoscopic
asymmetrical structures such as insect appendages could be designed by mixing
in specifically engineered meta-objects. This control of field quantities
using metamaterial concepts for shape-dependent separation of heterogeneous
suspension mixtures or adapting self-assembly processes with tailored acoustic
beams is likely to provide a fresh approach to the field of acoustophoresis.

This research was supported by the Australian Research Council through
Discovery Projects DP200101708 and DP200100358.

## References

  * Bruus [2011] H. Bruus, Lab on a Chip 11, 3742 (2011).
  * Bruus [2012a] H. Bruus, Lab on a Chip 12, 1014 (2012a).
  * Lenshof  _et al._ [2012] A. Lenshof, C. Magnusson, and T. Laurell, Lab on a Chip 12, 1210 (2012).
  * Bruus [2012b] H. Bruus, Lab on a Chip 12, 1578 (2012b).
  * Dual  _et al._ [2012] J. Dual, P. Hahn, I. Leibacher, D. Möller, T. Schwarz, and J. Wang, Lab on a Chip 12, 4010 (2012).
  * Wiklund  _et al._ [2012] M. Wiklund, R. Green, and M. Ohlin, Lab on a Chip 12, 2438 (2012).
  * Hartono  _et al._ [2011] D. Hartono, Y. Liu, P. L. Tan, X. Y. S. Then, L.-Y. L. Yung, and K.-M. Lim, Lab on a Chip 11, 4072 (2011).
  * Mohapatra  _et al._ [2018] A. R. Mohapatra, S. Sepehrirahnama, and K.-M. Lim, Physical Review E 97, 053105 (2018).
  * Augustsson  _et al._ [2012] P. Augustsson, C. Magnusson, M. Nordin, H. Lilja, and T. Laurell, Analytical chemistry 84, 7954 (2012).
  * Bernassau  _et al._ [2014] A. Bernassau, P. Glynne-Jones, F. Gesellchen, M. Riehle, M. Hill, and D. Cumming, Ultrasonics 54, 268 (2014).
  * Antfolk  _et al._ [2015] M. Antfolk, C. Magnusson, P. Augustsson, H. Lilja, and T. Laurell, Analytical chemistry 87, 9322 (2015).
  * Marzo  _et al._ [2015] A. Marzo, S. A. Seah, B. W. Drinkwater, D. R. Sahoo, B. Long, and S. Subramanian, Nature communications 6, 8661 (2015).
  * Wijaya  _et al._ [2016] F. B. Wijaya, A. R. Mohapatra, S. Sepehrirahnama, and K.-M. Lim, Microfluidics and Nanofluidics 20, 69 (2016).
  * Marzo and Drinkwater [2019] A. Marzo and B. W. Drinkwater, PNAS 116, 84 (2019).
  * Polychronopoulos and Memoli [2020] S. Polychronopoulos and G. Memoli, Scientific reports 10, 1 (2020).
  * King [1934] L. V. King, Proceedings of the Royal Society of London. Series A-Mathematical and Physical Sciences 147, 212 (1934).
  * Yosioka and Kawasima [1955] K. Yosioka and Y. Kawasima, Acta Acustica united with Acustica 5, 167 (1955).
  * Gorkov [1962] L. P. Gorkov, Soviet Physics - Doklady 6, 773 (1962).
  * Doinikov [1994a] A. A. Doinikov, Journal of Fluid Mechanics 267, 1 (1994a).
  * Settnes and Bruus [2012] M. Settnes and H. Bruus, Physical Review E 85, 016327 (2012).
  * Silva and Bruus [2014] G. T. Silva and H. Bruus, Physical Review E 90, 063007 (2014).
  * Marston [2006] P. L. Marston, The Journal of the Acoustical Society of America 120, 3518 (2006).
  * Mitri [2015a] F. Mitri, Wave Motion 57, 231 (2015a).
  * Garcia-Sabaté  _et al._ [2014] A. Garcia-Sabaté, A. Castro, M. Hoyos, and R. González-Cinca, The Journal of the Acoustical Society of America 135, 1056 (2014).
  * Doinikov [1994b] A. A. Doinikov, Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 447, 447 (1994b).
  * Foresti  _et al._ [2012] D. Foresti, M. Nabavi, and D. Poulikakos, Journal of Fluid Mechanics 709, 581 (2012).
  * Wijaya and Lim [2015] F. B. Wijaya and K.-M. Lim, Acta Acustica united with Acustica 101, 531 (2015).
  * Mitri [2015b] F. Mitri, Journal of Applied Physics 118, 214903 (2015b).
  * Wei  _et al._ [2004] W. Wei, D. B. Thiessen, and P. L. Marston, The Journal of the Acoustical Society of America 116, 201 (2004).
  * Xie and Wei [2004] W. Xie and B. Wei, Physical Review E 70, 046611 (2004).
  * Garbin  _et al._ [2015] A. Garbin, I. Leibacher, P. Hahn, H. Le Ferrand, A. Studart, and J. Dual, The Journal of the Acoustical Society of America 138, 2759 (2015).
  * Sepehrirahnama and Lim [2020a] S. Sepehrirahnama and K.-M. Lim, Physical Review E 102, 043307 (2020a).
  * Quan  _et al._ [2018] L. Quan, Y. Ra’di, D. L. Sounas, and A. Alù, Physical Review Letters 120, 254301 (2018).
  * Jordaan  _et al._ [2018] J. Jordaan, S. Punzet, A. Melnikov, A. Sanches, S. Oberst, S. Marburg, and D. A. Powell, Applied Physics Letters 113, 224102 (2018).
  * Melnikov  _et al._ [2019] A. Melnikov, Y. K. Chiang, L. Quan, S. Oberst, A. Alù, S. Marburg, and D. Powell, Nature communications 10, 1 (2019).
  * Chiang  _et al._ [2020] Y. K. Chiang, S. Oberst, A. Melnikov, L. Quan, S. Marburg, A. Alù, and D. A. Powell, Physical Review Applied 13, 064067 (2020).
  * Ni  _et al._ [2019] H. Ni, X. Fang, Z. Hou, Y. Li, and B. Assouar, Physical Review B 100, 104104 (2019).
  * Sieck  _et al._ [2017] C. F. Sieck, A. Alù, and M. R. Haberman, Physical Review B 96, 104303 (2017).
  * Su and Norris [2018] X. Su and A. N. Norris, Physical Review B 98, 174305 (2018).
  * Melnikov  _et al._ [2020] A. Melnikov, M. Maeder, N. Friedrich, Y. Pozhanka, A. Wollmann, M. Scheffler, S. Oberst, D. Powell, and S. Marburg, The Journal of the Acoustical Society of America 147, 1491 (2020).
  * Jiang  _et al._ [2017] X. Jiang, Y. Li, and L. Zhang, The Journal of the Acoustical Society of America 141, EL363 (2017).
  * sup [2022] See supplemental material at [url will be inserted by publisher] for analytical derivation of the force and torque in terms of orientation angle, details of numerical modelling of polarizability coefficients, and verification study of the effects of viscous losses., (2022).
  * Leao-Neto  _et al._ [2021] J. P. Leao-Neto, M. Hoyos, J.-L. Aider, and G. T. Silva, The Journal of the Acoustical Society of America 149, 285 (2021).
  * Leão-Neto  _et al._ [2020] J. P. Leão-Neto, J. H. Lopes, and G. T. Silva, The Journal of the Acoustical Society of America 147, 2177 (2020).
  * Marston and Zhang [2017] P. L. Marston and L. Zhang, The Journal of the Acoustical Society of America 141, 3042 (2017).
  * Lima  _et al._ [2020] E. B. Lima, J. P. Leão-Neto, A. S. Marques, G. C. Silva, J. H. Lopes, and G. T. Silva, Physical Review Applied 13, 064048 (2020).
  * Jerome  _et al._ [2021] T. S. Jerome, Y. A. Ilinskii, E. A. Zabolotskaya, and M. F. Hamilton, The Journal of the Acoustical Society of America 149, 2081 (2021).
  * Lopes  _et al._ [2020] J. H. Lopes, E. B. Lima, J. P. Leão-Neto, and G. T. Silva, Physical Review E 101, 043102 (2020).
  * Lima and Glauber [2021] E. B. Lima and T. S. Glauber, The Journal of the Acoustical Society of America 150, 376 (2021).
  * Zhang and Marston [2011a] L. Zhang and P. L. Marston, Physical Review E 84, 035601 (2011a).
  * Zhang and Marston [2011b] L. Zhang and P. L. Marston, Physical Review E 84, 065601 (2011b).
  * Zhang [2018] L. Zhang, Physical Review Applied 10, 034039 (2018).
  * Fan and Zhang [2019] X.-D. Fan and L. Zhang, Physical Review Applied 11, 014055 (2019).
  * Fan and Zhang [2021] X.-D. Fan and L. Zhang, The Journal of the Acoustical Society of America 150, 102 (2021).
  * Sepehrirahnama  _et al._ [2021] S. Sepehrirahnama, S. Oberst, Y. Chiang, and D. Powell, Physical Review E 104, 065003 (2021).
  * Marston  _et al._ [2006] P. L. Marston, W. Wei, and D. B. Thiessen, in  _AIP Conference Proceedings_ , Vol. 838 (AIP, 2006) pp. 495–499.
  * Fan  _et al._ [2008] Z. Fan, D. Mei, K. Yang, and Z. Chen, The Journal of the Acoustical Society of America 124, 2727 (2008).
  * Sepehrirahnama and Lim [2020b] S. Sepehrirahnama and K.-M. Lim, Microfluidics and Nanofluidics 24, 1 (2020b).
  * Mitri [2009] F. Mitri, Ultrasonics 49, 794 (2009).
  * Liu and Lim [2011] Y. Liu and K.-M. Lim, Lab on a Chip 11, 3167 (2011).


2022

[1,2]Vahid Partovi Nia

[1]Huawei Noah’s Ark Lab, Montreal Research Centre, 7101 Park Avenue,
Montreal, H3N 1X9, Quebec, Canada

2]Department of Mathematics and Industrial Engineering, Polytechnique
Montreal, 2500 Chem. Polytechnique, Montreal, H3T 1J4, Quebec, Canada

3]Department of Mathematics and Statistics, McGill University, 805 Sherbrooke
Street West, Montreal, H3A 0B9, Quebec, Canada

# Training Integer-Only Deep Recurrent Neural Networks

<EMAIL_ADDRESS>Eyyüb Sari<EMAIL_ADDRESS>Vanessa Courville
<EMAIL_ADDRESS>Masoud Asgharian<EMAIL_ADDRESS>*
[  [

###### Abstract

Recurrent neural networks (RNN) are the backbone of many text and speech
applications. These architectures are typically made up of several
computationally complex components such as; non-linear activation functions,
normalization, bi-directional dependence and attention. In order to maintain
good accuracy, these components are frequently run using full-precision
floating-point computation, making them slow, inefficient and difficult to
deploy on edge devices. In addition, the complex nature of these operations
makes them challenging to quantize using standard quantization methods without
a significant performance drop. We present a quantization-aware training
method for obtaining a highly accurate integer-only recurrent neural network
(iRNN). Our approach supports layer normalization, attention, and an adaptive
piecewise linear (PWL) approximation of activation functions, to serve a wide
range of state-of-the-art RNNs. The proposed method enables RNN-based language
models to run on edge devices with $2\times$ improvement in runtime, and
$4\times$ reduction in model size while maintaining similar accuracy as its
full-precision counterpart.

###### keywords:

Recurrent Neural Network, LSTM, Model compression, Quantization, NLP, ASR

##  1 Introduction

Recurrent Neural Network (RNN) rumelhart1986rnn  architectures such as Long-
Short Term Memory (LSTM) hochreiter1997lstm  or Gated Recurrent Units (GRU)
cho2014gru  are the foundation of language modeling applications. They can be
found in most large-scale systems such as neural machine translation
chen2018translstm ; wang2019accelerating , sentiment analysis
DBLP:journals/corr/abs-1801-07883 , image captioning You_2016_CVPR ; Mao-mRNN
and on-device systems such as Automatic Speech Recognition (ASR)
he2019streaming . Despite their pervasive use in academia and industry, RNN
architectures require more elaborated studies in order to improve inference
efficiency on-device.

Recently, Transformer-based models vaswani2017transformer ; devlin2019bert
have been proposed as an alternative to RNNs. Transformer architectures use an
encoder-decoder model, where the encoder maps an input data sequence to
learned vector representation, while the decoder uses the learned data to
produce a single output. Transformers differ from RNNs because sequential data
(such as words in a sentence) can be processed through the transformer’s
encoder concurrently, allowing for parallel training, not present in RNN
models. However, RNNs, especially LSTMs, are still powerful tools for building
highly accurate industry models. The combination of LSTM and transformer
architectures has provided impressive results recently; chen2018translstm ,
and wang2019accelerating  proposed an encoder-decoder neural machine
translator which harvests the inherent parallelism of transformers for its
encoder while utilizing LSTMs for a faster decoder. The LSTM decoder achieved
a faster inference runtime than a transformer decoder in these studies. It was
shown, for these cases, that Transformers may lose their parallelization
capabilities when generating outputs based on conditional inputs. Moreover,
they also suffer from quadratic computation complexity in terms of the
sequence length wang2019accelerating . kasai2021t2r  took transformers one
step further and fine-tuned a pre-trained transformer into RNNs, benefiting
the best of the two worlds.

When running these models on edge devices, such as mobile or IoT devices, the
physical hardware constraints must also be taken into consideration. In many
of these devices, the number of computing cores is limited to a handful of
processors, in which parallel-friendly transformer-based models lose their
advantage. There have been several studies on quantizing transformers to adapt
them for edge devices, but RNNs are largely ignored. Deploying RNN-based
chatbot, conversational agent, and ASR on edge devices with limited memory and
energy requires further computational improvements. The 8-bit integer neural
networks quantization jacob2017quantization  for convolutional architectures
(CNNs) is shown to be an almost free lunch approach to tackle the memory,
energy, and latency costs, with a negligible accuracy drop
krishnamoorthi2018quantizing .

Figure 1: Example of an integer-only LSTM cell (iLSTM). Layer normalization
changes to quantized integer friendly MadNorm (QMadNorm), full-precision
matrix multiplications change to integer matrix multiplications (QMatmul),
sigmoid, and tanh activation functions are replaced by their corresponding
piecewise linear (PWL) approximations.

Intuitively, quantizing RNNs is more challenging because the errors introduced
by quantization will propagate in two directions, i) to the successive layers,
like feedforward networks ii) across timesteps. Furthermore, RNN cells are
computationally more complex; they include several element-wise additions and
multiplications. They also have different activation functions that rely on
the exponential function, such as sigmoid and hyperbolic tangent (tanh).

Accurate fully-integer RNNs call for a new cell built using integer-friendly
operations. Our primary motivation is to enable integer-only inference of RNNs
on edge AI computing hardware with no floating-point units, so we constrained
the new LSTM cell to include only integer operations.

We propose a fully integer LSTM cell (iLSTM) in which its inference requires
integer-only computation units, see Figure 1. Our method can be applied to any
RNN architecture, but we focus on LSTM networks which are the most commonly
used RNNs.

Our contributions can be summarized as

  * •

providing a quantization-aware piece-wise linear approximation algorithm to
convert non-linear activation functions (e.g., sigmoid and tanh) to integer-
only equivalents

  * •

introducing an integer-friendly normalization layer based on mean absolute
deviation (MAD),

  * •

proposing integer-only attention,

  * •

combining these modules into an integer-only LSTM cell.

We also implement our method on a Huawei P20Pro smartphone, effectively
showing $2\times$ speedup and $4\times$ memory compression, hence providing
evidence that our method allows more RNN-based applications, such as ASR, to
run on edge devices.

##  2 Quantization

Quantization is a process whereby an input set is mapped to a lower resolution
discrete set, called the quantization set $\mathcal{Q}$. The mapping is
performed from floating-points to integers (e.g. float32 to int8) or from a
dense integer to another integer set with lower cardinality, e.g. int32 to
int8. We follow the Quantization-Aware Training (QAT) scheme described in
jacob2017quantization . Given $x\in[x_{\min},x_{\max}]$, we define the
quantization process as

$\displaystyle q_{x}=\mathrm{q}(x)$
$\displaystyle=\Big{\lfloor}\frac{x}{S_{x}}\Big{\rceil}+Z_{x}$ (1)
$\displaystyle r_{x}=\mathrm{r}(x)$ $\displaystyle=S_{x}(q_{x}-Z_{x})$ (2)
$\displaystyle S_{x}=\frac{x_{\max}-x_{\min}}{2^{b}-1},~{}~{}$ $\displaystyle
Z_{x}=\Big{\lfloor}\frac{-x_{\min}}{S_{x}}\Big{\rceil}$ (3)

where the input is clipped between $x_{\min}$ and $x_{\max}$ beforehand;
$\lfloor.\rceil$ is the round-to-nearest function; $S_{x}$ is the scaling
factor (also known as the step-size); $b$ is the bitwidth, e.g. $b=8$ for
8-bit quantization, $b=16$ for 16-bit quantization; $Z_{x}$ is the zero-point
corresponding to the quantized value of 0 (note that zero should always be
included in $[x_{\min},x_{\max}]$); $\mathrm{q}(x)$ quantizes $x$ to an
integer number and $\mathrm{r}(x)$ gives the floating-point value that
$\mathrm{q}(x)$ represents, i.e. $\mathrm{r}(x)\approx x$. We refer to
$\\{x_{\min},x_{\max},b,S_{x},Z_{x}\\}$ as quantization parameters of $x$.
Note that for inference, $S_{x}$ is expressed as a fixed-point integer number
rather than a floating-point number, allowing for integer-only arithmetic
computations jacob2017quantization .

The following example showcases how to apply Equations (1-3). Given a number
within the previously defined range i.e. $x\in[-1,1]$ and a target bit-width
of $b=8$, $S_{x}$ and $Z_{x}$ are computed as;

$\displaystyle S_{x}$
$\displaystyle=\frac{x_{\max}-x_{\min}}{2^{b}-1}=\frac{2}{255}\approx 0.0078$
$\displaystyle Z_{x}$
$\displaystyle=\Big{\lfloor}\frac{-x_{\min}}{S_{x}}\Big{\rceil}=\Big{\lfloor}\frac{1}{S_{x}}\Big{\rceil}=128.$

Assuming $x=0.2$,

$\displaystyle\mathrm{q}(x)$
$\displaystyle=\Big{\lfloor}\frac{x}{S_{x}}\Big{\rceil}+Z_{x}=\Big{\lfloor}\frac{0.2}{0.0078}\Big{\rceil}+128=154,$

hence $q_{x}=154$. The corresponding floating-point representation of $q_{x}$,
$r_{x}$, is

$\displaystyle\mathrm{r}(x)$ $\displaystyle=S_{x}(q_{x}-Z_{x})$
$\displaystyle=0.0078(154-128)=0.2028.$

Therefore, after quantizing the floating-point value $x=0.2$, its dequantized
value $r_{x}$ is $0.2028$, so the error introduced is $\epsilon=\mid
x-r_{x}\mid=0.0028$. The purpose of QAT is to try to force the network to be
insensitive to noise introduced by these $\epsilon$ errors. The error
introduced by quantization is bounded by $\epsilon\leq\frac{S_{x}}{2}$. The
smaller the scaling factor, the better the precision. The solutions to make
the scaling factors small are either to have a smaller range
$[x_{\min},x_{\max}]$ or to use a bigger bitwidth $b$. In practice, for a
given bitwidth, we would like to have a large range, including all possible
values the network can produce while reducing the error as much as possible.

###  2.1 Integer-only arithmetic

In order to map floating-point operations to integer-only operators, all
element-wise multiplication and additions must be adapted to use Equations
(1-3). We describe this mapping below. jacob2017quantization  describes how to
do this for matrix-multiplications and other operators.

At inference, the inputs to any operation are already quantized. We define two
quantized inputs; $q_{a}$ and $q_{b}$ and a quantized output, $q_{c}$, along
with their respective scaling factors $S_{a}$, $S_{b}$, $S_{c}$ and zero-
points $Z_{a}$, $Z_{b}$, $Z_{c}$.

The multiplication operation does not require $q_{a}$ and $q_{b}$ to share
quantization parameters therefore, it is defined as

$q_{c}=\Big{\lfloor}\frac{S_{a}S_{b}}{S_{c}}\Big{(}q_{a}q_{b}-q_{a}Z_{b}-q_{b}Z_{a}+Z_{a}Z_{b}\Big{)}\Big{\rceil}+Z_{c}$
(4)

where the multiplication $q_{a}q_{b}$ is a $b$-bit multiplication (example
given in Appendix B.1).

Addition can take two forms based on the quantization parameters of $q_{a}$
and $q_{b}$. If they share the same quantization parameters (i.e.
$S_{b}=S_{a}$ and $Z_{b}=Z_{a}$), it can be computed as follows,

$\displaystyle
q_{c}=\Big{\lfloor}\frac{S_{a}}{S_{c}}(q_{a}+q_{b}-2Z_{a})\Big{\rceil}+Z_{c}$
(5)

where $q_{a}+q_{b}$ is a $b$-bit addition. However, if they do not share the
same quantization parameters,

$\displaystyle
q_{c}=\Big{\lfloor}\frac{S_{a}}{S_{c}}(q_{a}-Z_{a})+\frac{S_{b}}{S_{c}}(q_{b}-Z_{b})\Big{\rceil}+Z_{c}$
(6)

rescaling $q_{a}$ and $q_{b}$ is needed, which means the addition has to be
computed with more bitwidth, e.g. 32-bit (examples given in Appendix B.2).

###  2.2 Inference

For fully-integer inference, all variables must be expressed as integers.
However, the scaling-factor $S_{x}$ described in Equation (3) is a floating-
point value. To solve this, $S_{x}$ can be represented as fixed-point
integers, computed off-line (after training but before inference). This
enables the inference to run on edge devices which typically have limited
power and data format support courville2019deep .

A fixed-point number is an integer representation of a floating-point number,
and is described using what is known as Q-format. A fixed-point integer having
format $Q_{i.f}$ is said to have one sign bit, $i$ integral bits and $f$
fractional bits, bitwidth $1+i+f$, a resolution of $2^{-f}$ and span the range
$[-2^{-f},2^{i}-2^{-f}]$. For instance, $Q_{3.4}$ represents numbers in
$[-8,7.9375]$ with resolution $0.0625$ and $Q_{0.30}$ represents numbers in
$[-1,0.99999]$ with $\approx 9\times 10^{-10}$. Assuming we do not have to
worry about overflows, quantizing a floating-point number $M$ to a fixed-point
integer $M_{fx}$ is given by,

$\displaystyle M_{fx}=\lfloor 2^{f}M\rceil.$ (7)

Converting a fixed-point integer $M_{fx}$ to a floating-point number $M$ is
given by,

$\displaystyle M=2^{-f}M_{fx}$ (8)

An equivalent expression to (8) is

$\displaystyle M=M_{fx}\gg f$ (9)

where $\gg$ is the right shift operator. This can be incorporated into
floating-point expressions, often of the form $q_{y}=\lfloor Mq\rceil+Z_{y}$,
which can be expressed as

$\displaystyle q_{y}$
$\displaystyle=\Big{\lfloor}2^{-f}M_{fx}q\Big{\rceil}+Z_{y}$ (10)
$\displaystyle=M_{fx}q\gg f+(M_{fx}q\gg(f-1))~{}\&~{}0\textrm{x}1+Z_{y}.$ (11)

Note that $\&~{}0\textrm{x}1$ masks every bit in of the number except the
rightmost one, e.g. $2\&~{}0\textrm{x}1=0$ because 2 is $(10)_{2}$ but
$3\&~{}0\textrm{x}1=1$ because 3 is $(11)_{2}$. It returns 0 if your number is
even and 1 if your number is odd. We extracted the negative sign out of the
binary shifts, as shifting on negative number has different behaviour
depending on the hardware and software implementation. Note that in practice,
when we know a constant can only be positive, it will be expressed as an
unsigned fixed-point integer and therefore the sign bit can be allocated to
either the integral or fractional part. More details on fixed-point arithmetic
are provided in Appendix B.3.

##  3 Recurrent Neural Network

###  3.1 Related Work

In recent years, neural network quantization research jacob2017quantization ;
courbariaux2016bnn ; darabi2018bnnplus ; esser2020lsq  has enabled exploration
into the quantization of RNN architectures. ott2016rnnlimitprec , which
explores low-bit quantization of weights for RNNs, showed that binarizing
weights leads to a massive accuracy drop, but ternarizing them keeps the model
performance. courbariaux2016bnn  demonstrates that quantizing RNNs to
extremely low bits is challenging; they quantize weights and matrix product to
4-bit, but other operations such activation functions are computed in full-
precision.

Low-precision RNNs has also been a recent focal point in industry research.
he2016effectivernn  introduces Bit-RNN and improves 1-bit and 2-bit RNNs
quantization by constraining values within a fixed range carefully; they keep
activation computation and element-wise operations in full-precision.
kapur2017lowprecrnn  builds upon Bit-RNN and propose a low-bit RNN with
minimal performance drop, but they increase the number of neurons to
compensate for performance drop; they run activation functions in floating-
point precision as well. hou2019normalization  quantizes LSTM weights to 1-bit
and 2-bit and empirically show that low-bit quantized LSTMs suffer from
exploding gradients. It has been shown that gradient explosion can be
alleviated using normalization layers and leads to the successful training of
low bit weights ardakani2018learning . sari2020normalization  studied the
effect of normalization in low-bit networks theoretically and proved that low-
bit training without normalization operation is mathematically impossible;
their work demonstrates the fundamental importance of involving normalization
layers in quantized networks.

wu2016gnmt  is a pioneering work in LSTM quantization, which demonstrates
speed-up inference of large-scale LSTM models with limited performance drop by
partially quantizing RNN cells. Their proposed method is tailored towards
specific hardware. They use an 8-bit integer for matrix multiplications and a
16-bit integer for tanh, sigmoid, and element-wise operations but do not
quantize attention. bluche2020sonosquantlstm  propose an effective 8-bit
integer-only LSTM cell for Keyword Spotting application on microcontrollers.
They enforce weights and activations to be symmetric on fixed ranges $[-4,4]$
and $[-1,1]$. This prior assumption about the network’s behavior restricts
generalizing their approach for many RNN models. They propose a look-up table
of 256 slots to represent the quantized tanh and sigmoid activations. However,
the look-up table memory requirement explodes for a bigger bitwidth. Their
solution does not serve complex tasks such as automatic speech recognition due
to large look-up table memory consumption. While demonstrating strong results
on the Keyword Spotting task, their assumptions on quantization range and
bitwidth make their method task-specific.

To explain the methodology in the following subsections, we use the common
linear algebra notation and use plain symbols to denote scalar values, e.g.
$x\in\mathbb{R}$, bold lower-case letters to denote vectors, e.g.
$\textbf{x}\in\mathbb{R}^{n}$, and bold upper-case letters to denote matrices,
e.g. $\textbf{X}\in\mathbb{R}^{m\times n}$. The element-wise multiplication is
represented by $\odot$.

###  3.2 Long short term memory

We define an LSTM cell as

$\displaystyle\begin{pmatrix}\mathbf{i}_{t}\\\ \mathbf{f}_{t}\\\
\mathbf{j}_{t}\\\ \mathbf{o}_{t}\end{pmatrix}$
$\displaystyle=\mathbf{W}_{x}\mathbf{x}_{t}+\mathbf{W}_{h}\mathbf{h}_{t-1},$
(12) $\displaystyle\mathbf{c}_{t}$
$\displaystyle=\sigma(\mathbf{f}_{t})\odot\mathbf{c}_{t-1}+\sigma(\mathbf{i}_{t})\odot\text{tanh}(\mathbf{j}_{t}),$
(13) $\displaystyle\mathbf{h}_{t}$
$\displaystyle=\sigma(\mathbf{o}_{t})\odot\text{tanh}(\mathbf{c}_{t}),$ (14)

where $\sigma(\cdot)$ is the sigmoid function; $n$ is the input hidden units
dimension, and $m$ is the state hidden units dimension;
$\mathbf{x}_{t}\in\mathbb{R}^{n}$ is the input for the current timestep
$t\in\\{1,...,T\\}$; $\mathbf{h}_{t-1}\in\mathbb{R}^{m}$ is the hidden state
from the previous timestep and $\mathbf{h}_{0}$ is initialized with zeros;
$\mathbf{W}_{x}\in\mathbb{R}^{4m\times n}$ is the input to state weight
matrix; $\mathbf{W}_{h}\in\mathbb{R}^{4m\times m}$ is the state to state
weight matrix;
$\\{\mathbf{i}_{t},\mathbf{f}_{t},\mathbf{o}_{t}\\}\in\mathbb{R}^{m}$ are the
pre-activations to the {input, forget, output} gates;
$\mathbf{j}_{t}\in\mathbb{R}^{m}$ is the pre-activation to the cell candidate;
$\\{\mathbf{c}_{t},\mathbf{h}_{t}\\}\in\mathbb{R}^{m}$ are the cell state and
the hidden state for the current timestep, respectively. We omit the biases
for the sake of notation simplicity. For a bidirectional LSTM (BiLSTM) the
output hidden state at timestep $t$ is the concatenation of the forward hidden
state $\overrightarrow{\mathbf{h}}_{t}$ and the backward hidden state
$\overleftarrow{\mathbf{h}}_{t}$, $[\overrightarrow{\mathbf{h}}_{t}\mathchar
24635\relax\;\overleftarrow{\mathbf{h}}_{t}]$.

###  3.3 LayerNorm

Layer normalization ba2016layernorm  standardizes inputs across the hidden
units dimension with zero location and unit scale. Given hidden units
$\mathbf{x}\in\mathbb{R}^{H}$, LayerNorm is defined as

$\displaystyle\mu=\frac{1}{H}\sum_{i=1}^{H}x_{i}\quad\quad\hat{x}_{i}=x_{i}-\mu$
(15)
$\displaystyle\sigma_{\mathrm{std}}^{2}=\frac{1}{H}\sum_{i=1}^{H}\hat{x}_{i}^{2}\quad\quad\sigma_{\mathrm{std}}=\sqrt{\sigma_{\mathrm{std}}^{2}}$
(16)
$\displaystyle\mathrm{LN}(\mathbf{x})_{i}=y_{i}=\frac{\hat{x_{i}}}{\sigma_{\mathrm{std}}}$
(17)

$\forall i,j1\leq i,j\leq n\quad\text{and}\quad c_{ij}=\sum a_{ik}b_{kj}$ (18)

where $\mu$ (15) is the hidden unit mean, $\hat{x}_{i}$ (15) is the centred
hidden unit $x_{i}$, $\sigma_{\mathrm{std}}^{2}$ (16) is the hidden unit
variance, and $y_{i}$ (17) is the normalized hidden unit. In practice, one can
scale $y_{i}$ by a learnable parameter $\gamma$ or shift by a learnable
parameter $\beta$. The LayerNormLSTM cell is defined as in ba2016layernorm .

###  3.4 Attention

Attention is often used in encoder-decoder RNN architectures
bahdanau2015attention ; chorowski2015speechattention ; wu2016gnmt . We employ
Bahdanau attention, also called additive attention bahdanau2015attention . The
attention mechanism allows the decoder network to attend to the variable-
length output states from the encoder based on their relevance to the current
decoder timestep. At each of its timesteps, the decoder extracts information
from the encoder’s states and summarizes it as a context vector,

$\displaystyle\mathbf{s}_{t}$
$\displaystyle=\sum_{i=1}^{T_{\mathrm{enc}}}\alpha_{ti}\odot\mathbf{h}_{\mathrm{enc}_{i}}$
(19) $\displaystyle\alpha_{ti}$
$\displaystyle=\frac{\text{exp}(e_{ti})}{\sum_{j=1}^{T_{\mathrm{enc}}}\text{exp}(e_{tj})}$
(20) $\displaystyle e_{ti}$
$\displaystyle=\mathbf{v}^{\top}\text{tanh}(\mathbf{W}_{q}\mathbf{h}_{t-1}+\mathbf{W}_{k}\mathbf{h}_{enc_{i}})$
(21)

where $\mathbf{s}_{t}$ is the context at decoder timestep $t$ which is a
weighted sum of the encoder hidden states outputs
$\mathbf{h}_{\mathrm{enc}_{i}}\in\mathbb{R}^{m_{\mathrm{enc}}}$ along encoder
timesteps $i\in\\{1,...,T_{\mathrm{enc}}\\}$ ; $0<\alpha_{ti}<1$ are the
attention weights attributed to each encoder hidden states based on the
alignments $e_{ti}\in\mathbb{R}$; $m_{\mathrm{dec}}$ and $m_{\mathrm{enc}}$
are respectively the decoder and encoder hidden state dimension;
$\\{\mathbf{W}_{q}\in\mathbb{R}^{m_{\mathrm{att}}\times
m_{\mathrm{dec}}},\mathbf{W}_{k}\in\mathbb{R}^{m_{\mathrm{att}}\times
m_{\mathrm{enc}}}\\}$ are the weights matrices of output dimension
$m_{\mathrm{att}}$ respectively applied to the query $\mathbf{h}_{t-1}$ and
the keys $\mathbf{h}_{\mathrm{enc}_{i}}$;
$\mathbf{v}\in\mathbb{R}^{m_{\mathrm{att}}}$ is a learned weight vector. The
context vector is incorporated into the LSTM cell by modifying (12) to

$\displaystyle\begin{pmatrix}\mathbf{i}\\\ \mathbf{f}\\\ \mathbf{j}\\\
\mathbf{o}\end{pmatrix}=\mathbf{W}_{x}\mathbf{x}_{t}+\mathbf{W}_{h}\mathbf{h}_{t-1}+\mathbf{W}_{s}\mathbf{s}_{t}$
(22)

where $\mathbf{W}_{s}\in\mathbb{R}^{4m_{\mathrm{dec}}\times
m_{\mathrm{enc}}}$.

##  4 Methodology

In this section, we describe our task-agnostic quantization-aware training
method to enable integer-only RNN (iRNN).

###  4.1 Integer-only activation

First, we need to compute activation functions without relying on floating-
point operations to take the early step towards an integer-only RNN. At
inference, the non-linear activation is applied to the quantized input
$q_{x}$, operations using integer-only arithmetic and output the quantized
result $q_{y}$. Clearly, given the activation function $f$,
$q_{y}=q(f(q_{x}))$; as the input and the activation and output are both
quantized, we obtain a discrete mapping from $q_{x}$ to $q_{y}$. There are
several ways to formalize this operation. The first solution is a Look-Up
Table (LUT), where $q_{x}$ is the index and $q_{y}=\text{LUT}[q_{x}]$. For
8-bit quantization, a LUT would require 256 bytes of storage ($2^{8}$ 8-bit
elements), however 16-bit quantization would require up to 128 kilobytes of
memory ($2^{16}$ 16-bit elements). Therefore, this method does not scale well
and is not cache-friendly for large indexing bitwidths. The second solution
approximates the full-precision activation function using a fixed-point
integer Taylor approximation, but the amount of computations grows as the
approximation order grows.

We propose to replace all non-linear activation functions of an LSTM cell with
integer-only arithmetic by approximating them using a Quantization-Aware PWL
function which defines the PWL pieces at training time. This enables the
precision of the PWL to be adapted to the task at hand and provides highly
accurate data-dependent activation approximation with fewer pieces.

Figure 2:  $\tanh$ approximations using quantization-aware PWLs with 4 pieces
(left panel), 16 pieces (right panel) using (23). The dashed cyan curves are
the true tanh functions, while the solid orange curves are its approximation
from Algorithm 1. Red dots are the knots. The more we add pieces, the better
the approximation is. Our algorithm can prioritize sections of the function
with more curvatures; see Appendix A Figure 3 for more instances.

A PWL is defined as follows,

$\displaystyle
g(x)=\sum_{i=1}^{N}\mathbbm{1}_{[k_{i},k_{i+1})}\Big{(}a_{i}(x-k_{i})+b_{i}\Big{)},$
(23)

where $N$ is the number of linear pieces defined by $N+1$ knots (also known as
cutpoints or breakpoints); $\\{a_{i},k_{i},b_{i}=f(k_{i})\\}$ are the slope,
the knot, and the intercept of the $i^{\text{th}}$ piece respectively;
$\mathbbm{1}_{A}(x)=1$ is the indicator function on $A$. The more the linear
pieces, the better the activation approximation is (see Figure 2). A PWL is
suitable for simple fixed-point integer operations. It only relies on basic
arithmetic operations and is easy to parallelize because the computation of
each piece is independent. Therefore, the challenge is to select the knot
locations that provide the best PWL approximation to the original function
$f$. Note in this regime, we only approximate the activation function on the
subset corresponding quantized inputs and not the whole full-precision range.
In our proposed method if $x=k_{i}$ then $g(x)=g(k_{i})=b_{i}$, i.e. recovers
the exact output $f(k_{i})$. Hence, if the PWL has $2^{b}$ knots (i.e.
$2^{b}-1$ pieces), it is equivalent to a look-up table representing the
quantized activation function. Thus, we constraint the knots to be a subset of
the quantized inputs of the function we are approximating (i.e.
$\\{k_{i}\\}_{i=1}^{N+1}\subseteq\mathcal{Q}$).

Algorithm 1 Piece-wise linear approximation of activation $f$

We propose a recursive greedy algorithm, Algorithm 1, to locate the knots
during the quantization-aware PWL. The algorithm starts with $2^{b}-1$ pieces
and recursively removes one knot at a time until it reaches the specified
number of pieces. The absolute differences between adjacent slopes are
computed, and the shared knot from the pair of slopes that minimizes the
absolute difference is removed; see Appendix A, Figure 4. Applied only once
per training step, the generic nature of this algorithm enables it to be
applied to any nonlinear function. Note that the PWL is specific to a given
set of quantization parameters, i.e. the quantization parameters are kept
frozen after its creation.

At inference, the quantization-aware PWL is computed as follows

$\displaystyle
q_{y}=\Big{\lfloor}\sum_{i=1}^{N}\mathbbm{1}_{[q_{k_{i}},q_{k_{i+1}})}\Big{(}\frac{S_{x}a_{i}}{S_{y}}(q_{x}-q_{k_{i}})+\frac{b_{i}}{S_{y}}\Big{)}\Big{\rceil}+Z_{y},$
(24)

where the constants are expressed as fixed-point integers.

###  4.2 Integer-only normalization

Normalization greatly helps the convergence of quantized networks
hou2019normalization ; sari2020normalization . There are many location and
scale measures to define a normalization operation. The commonly used measure
of dispersion is the standard deviation to define normalization, which is
imprecise and costly to compute on integer-only hardware. However, the mean
absolute deviation (MAD) is integer-friendly and defined as

$d=\frac{1}{H}\sum_{i=1}^{H}\left\lvert x_{i}-\mu\right\rvert$ (25)

While the mean minimizes the standard deviation, the median minimizes MAD. We
suggest measuring deviation from the mean for two reasons: i) the median is
computationally more expensive ii) the absolute deviation from the mean is
closer to the standard deviation. For Gaussian data, the MAD is $\approx
0.8\sigma_{\mathrm{std}}$ so that it might be exchanged with standard
deviation. We propose to LayerNorm in LSTM with MAD instead of standard
deviation and refer to it as MadNorm., where (17) is replaced by

$y_{i}=\frac{\hat{x}_{i}}{d}.$ (26)

MadNorm involves simpler operations, as there is no need to square and no need
to take the square root while taking absolute value instead of these two
operations is much cheaper. The values $\\{\mu,\hat{x}_{i},d,y_{i}\\}$ are
8-bit quantized and computed as follows

$\displaystyle q_{\mu}$
$\displaystyle=\Big{\lfloor}\frac{S_{x}}{S_{\mu}N}\Big{(}\sum_{i=1}^{N}q_{x_{i}}-NZ_{x}\Big{)}\Big{\rceil}+Z_{\mu},$
(27) $\displaystyle q_{\hat{x}_{i}}$
$\displaystyle=\Big{\lfloor}\frac{S_{x}}{S_{\hat{x}}}(q_{x_{i}}-Z_{x})-\frac{S_{\mu}}{S_{\hat{x}}}(q_{\mu}-Z_{\mu})\Big{\rceil}+Z_{\hat{x}},$
(28) $\displaystyle q_{d}$
$\displaystyle=\Big{\lfloor}\frac{S_{\hat{x}}}{S_{d}N}\sum_{i=1}^{N}\left\lvert
q_{\hat{x}_{i}}-Z_{\hat{x}}\right\rvert\Big{\rceil}+Z_{d},$ (29)
$\displaystyle q_{y_{i}}$
$\displaystyle=\Big{\lfloor}\frac{\frac{S_{\hat{x}}}{S_{y}S_{d}}(q_{\hat{x_{i}}}-Z_{\hat{x}})}{\mathrm{\max}(q_{d},1)}\Big{\rceil}+Z_{y},$
(30)

where all floating-point constants can be expressed as fixed-point integer
numbers, allowing for integer-only arithmetic computations. Note that (27-30)
are only examples of ways to perform integer-only arithmetic for MadNorm, and
may change depending on the software implementation and the target hardware.
We propose to quantize $\\{\mathbf{v},\mathbf{W}_{q},\mathbf{W}_{k}\\}$ to
8-bit. The vectors $\mathbf{h}_{t-1}$ and $\mathbf{h}_{\mathrm{enc}_{i}}$ are
quantized thanks to the previous timestep and/or layer. The matrix
multiplications in (21) are performed in 8-bit, and their results are
quantized to 8-bit, each with its quantization parameters. Since those matrix
multiplications do not share the same quantization parameters, the sum (21)
requires proper rescaling, and the result is quantized to 16-bit. We found
that 8-bit quantization adds too much noise, thus preventing the encoder-
decoder model from working correctly. The tanh function in (21) is computed
using quantization-aware PWL, and its outputs are quantized to 8-bit. The
alignments $e_{ti}$ are quantized to 16-bit (21). The exponential function in
$\alpha_{ti}$ is computed using a quantization-aware PWL, with its outputs
quantized to 8-bit. We found that quantizing the softmax denominator (20) to
8-bit introduces too much noise and destroys attention. 8-bit attention does
not offer enough flexibility and prevents fine-grained decoder attention to
the encoder. We left the denominator in a 32-bit integer value. The context
vector $\mathbf{s}_{t}$ is quantized to 8-bit. Note that in practice we shift
the inputs of softmax for numerical stability reasons $e_{ti}\leftarrow
e_{ti}-\max\limits_{j}e_{tj}$.

######  Theorem 1 (Scale convergence).

Suppose $\\{X_{i}\\}_{i=1}^{\infty}$ is a sequence of integrable pairwise
independent random variables defined on the same probability space
$(\Omega,\mathcal{F},P)$ with $\mu=\mathrm{I\\!E}(X_{i})$. Then $D_{n}={1\over
n}\sum_{i=1}^{n}\left\lvert X_{i}-\mu\right\rvert$ converges almost surely to
$\tilde{\sigma}=\mathrm{I\\!E}(\left\lvert X-\mu\right\rvert).$

Proof: Integrability assumption implies existence of
$\mu=\mathrm{I\\!E}(X_{i})<\mathrm{I\\!E}(\left\lvert
X_{i}\right\rvert)<\infty$ and hence existence of $\mathrm{I\\!E}\left\lvert
X_{i}-\mu\right\rvert<\mathrm{I\\!E}\left\lvert
X_{i}\right\rvert+\left\lvert\mu\right\rvert<\infty$. The desired result then
follows from the Strong Law of Large Numbers to $Y_{i}=\left\lvert
X_{i}-\mu\right\rvert.\ \blacksquare$

One may prove the central limit theorem by replacing absolute integrability
with square integrability, and exchanging pairwise independence with mutual
independence. Convergence to the population scale $\tilde{\sigma}$ in Theorem
1 paves the way to show that our MadNorm enjoys a concentration inequality
similar to LayerNorm.

######  Theorem 2 (Concentration inequality).

Suppose $X$ defined on $(\Omega,\mathcal{F},P)$ is an integrable random
variable with mean $\mu$. Then for a positive $k$,

$P\left(\left\lvert X-\mu\over\tilde{\sigma}\right\rvert<k\right)\geq
1-{1\over k}$

Proof: The result follows from Markov’s inequality, see Sari_iRNN_ICPRAM2022
for details. $\blacksquare$

Theorem 2 assures that independent of the distribution of data, this MadNorm
brings the mass of the distribution around the origin. This is somehow
expected from any normalization method. Therefore, it is not surprising that
LayerNorm also has a similar property. Therefore, LayerNorm and MadNorm assure
that the tail probability far from the origin is negligible.

There is a slight difference between the concentration inequality of LayerNorm
and MadNorm. The LayerNorm provides a tighter bound, i.e., the bound in
Theorem 2 changes from $1-{1\over k}$ to $1-{1\over k^{2}}$, but it also
requires more assumptions like the square integrability of $X$.

###  4.3 Integer-only attention

Attention plays a crucial role in modern encoder-decoder architectures. The
decoder relies on attention to extract information from the encoder and
provide predictions. Attention is the bridge between the encoder and the
decoder. Careless quantization of attention breaks apart the decoder due to
quantization noise.

###  4.4 Integer-only LSTM

A standard LSTM cell comprises matrix multiplications, element-wise additions,
element-wise multiplications, tanh, and sigmoid activations (12 \- 14). We
quantize the weights matrices $\mathbf{W}_{x}$ and $\mathbf{W}_{h}$ to 8-bit.
The inputs $\mathbf{x}_{t}$ and hidden states $\mathbf{h}_{t-1}$ are already
8-bit quantized from the previous layer and from the previous timestep. The
cell states $\mathbf{c}_{t}$ are theoretically unbounded (13); therefore the
amount of quantization noise potentially destroys the information carried by
$\mathbf{c}_{t}$, if it spans a large range. When performing QAT on some pre-
trained models, it is advised to quantize $\mathbf{c}_{t}$ to 16-bit.
Therefore, $\mathbf{c}_{t}$ is 8-bit quantized unless stated otherwise but can
be quantized to 16-bit if necessary. Matrix multiplications in (12) are
performed with 8-bit arithmetic, and their outputs are quantized to 8-bit
based on their respective quantization parameters. The sum between the two
matrix multiplications outputs in (12) requires proper rescaling because they
do not share the same quantization parameters.

The sum results are quantized to 8-bit; however, 16-bit quantization might be
necessary for complex tasks. The sigmoid and tanh activations in (13) and (14)
are replaced with their own quantization-aware PWL, and their output is always
quantized to 8-bit. The element-wise multiplications operations are
distributive, and sharing quantization parameters is not required. In (13),
the element-wise multiplications are quantized to 8-bit, but can be quantized
to 16-bit if $\mathbf{c}_{t}$ is quantized to 16-bit as well; the element-wise
additions are quantized based on $\mathbf{c}_{t}$’s bitwidth (i.e. 8-bit or
16-bit).

The element-wise multiplications between sigmoid and tanh in (14) are always
quantized to 8-bit because $\mathbf{h}_{t}$ are always quantized to 8-bit. We
obtain an integer-only arithmetic LSTM cell following this recipe, see Figure.
1. For LSTM cells with LayerNorm quantized MadNorm layers are used instead of
LayerNorm. Appendix A.1 provides details about the quantization of other types
of layers in an LSTM model.

##  5 Experimental results

We evaluate our proposed method, iRNN, on language modeling and automatic
speech recognition. We also implemented our approach on a smartphone to
benchmark inference speedup, see Section 5.4.

###  5.1 Language modeling on PTB

LayerNorm LSTM | val | test
---|---|---
Full-precision | $98.58\pm 0.35$ | $94.84\pm 0.21$
PWL4 | $101.40\pm 0.70$ | $98.11\pm 0.75$
PWL8 | $98.14\pm 0.11$ | $95.03\pm 0.16$
PWL16 | $98.09\pm 0.06$ | $94.92\pm 0.05$
PWL32 | $\mathbf{97.97}\pm 0.01$ | $\mathbf{94.81}\pm 0.02$
Table 1: Word-level perplexities on PTB with a LayerNorm LSTM and quantized
models with a different number of PWL pieces. LayerNorm is replaced with
MadNorm for the quantized models (iRNN). Best results are averaged across 3
runs $\pm$ standard deviation.

As a proof of concept, we perform several experiments on full-precision and
entirely 8-bit quantized models on the Penn TreeBank (PTB) dataset
marcus1993ptb . We report perplexity per word as a performance metric. For the
quantized models, the LayerNorm is replaced with MadNorm. We do not train
full-precision models with MadNorm to make our method comparable with common
full-precision architectures. The quantized PWL models are initialized from
the best full-precision checkpoint and train for another 100 epochs.
Experimental setup details are provided in Appendix A.2.1. We can draw two
conclusions from the results presented in Table 1, i) replacing LayerNorm with
MadNorm does not destroy model performance, ii) using eight linear pieces is
enough to retain the performance of the model, but adding more linear pieces
improves the performance. Furthermore, we could obtain even superior results
in the quantized model compared to the full-precision model because of the
regularization introduced by quantization errors.

###  5.2 Language modeling on WikiText2

Mogrifier LSTM | val | test
---|---|---
Full-precision | $60.27\pm 0.34$ | $58.02\pm 0.34$
PWL8 | $60.91\pm 0.04$ |  $58.54\pm$ 0.07
PWL16 | $60.65\pm 0.09$ | $58.21\pm 0.08$
PWL32 | $\mathbf{60.37}\pm 0.03$ | $\mathbf{57.93}\pm 0.07$
Table 2: Word-level perplexities on WikiText2 with Mogrifier LSTM and
quantized models with different number of PWL pieces. Best results are
averaged across 3 runs $\pm$ standard deviations.

We evaluated our proposed method on the WikiText2 dataset merity2016wikitext
with a state-of-the-art RNN, Mogrifier LSTM melis2020mogrifier . The original
code111https://github.com/deepmind/lamb was written in TensorFlow, and we
reimplemented our version in PyTorch by staying as close as possible to the
TensorFlow version. We follow the experimental setup from the
authors222https://github.com/deepmind/lamb/blob/254a0b0e330c44e00cf535f98e9538d6e735750b/lamb/experiment/mogrifier/config/c51c838b33a5+_tune_wikitext-2_35m_lstm_mos2_fm_d2_arms/trial_747/config
as we found it critical to get similar results. We use a two-layer Mogrifier
LSTM. The setup and hyper-parameters for the experiments can be found in
Appendix A.2.2 to save some space. We present our results averaged over 3 runs
in Table 2. We use the best full-precision model, which scores $59.93$
perplexity to initialize the quantized models. Our method can produce an 8-bit
quantized integer-only Mogrifier LSTM with similar performance to the full-
precision model with only about $0.3$ perplexity increase for the quantized
model with a PWL of 32 pieces and a maximum of about $0.9$ perplexity increase
with several pieces as low as 8. Interestingly, a pattern emerged by doubling
the number of pieces, as we get a decrease in perplexity by about $0.3$. We
also perform a thorough ablation study of our method in Table 5. Surprisingly,
we found that stochastic weight averaging for quantized models exhibits the
same behavior as for full-precision models and improved performance thanks to
regularization. While experiments on the PTB dataset demonstrated the
potential of our method, these experiments on WikiText2 show that our proposed
method can stay on par with state-of-the-art RNN models. This helps emphasize
that our fully-integer iRNN model is capable of maintaining comparable
accuracy to the full-precision model.

###  5.3 ASR on LibriSpeech

ESPRESSO LSTM | set | clean | other
---|---|---|---
Full-precision | dev | 2.99 | 8.77
iRNN PWL96* | dev | 3.73 | 10.02
Full-precision | test | 3.37 | 9.49
iRNN PWL96* | test | 4.11 | 10.71
Table 3: WER% on LibriSpeech with ESPRESSO LSTM (Encoder-Decoder LSTM with
Attention) with LM shallow fusion. *(160 pieces were used for the exponential
function)

ASR is a critical edge AI application and is challenging due to the nature of
the task. Voice is diverse as the human voice may vary in pitch, accent,
pronunciation style, voice volume, etc. While we showed our method is working
for a competitive language modeling task, one can argue that ASR is a more
practical and, at the same time, more difficult task for edge and IoT
applications. Therefore, we experiment on an ASR task based on the setup of
wang2019espresso  and their ESPRESSO
framework333https://github.com/freewym/espresso. We used an LSTM-based
Attention Encoder-Decoder (ESPRESSO LSTM) trained on the strong ASR
LibriSpeech dataset panayotov2015librispeech . Experiments setup and hyper-
parameters are provided in Appendix A.2.3. We initialize the quantized model
from the pre-trained full-precision ESPRESSO LSTM. Our early experiments
confirmed that quantizing the model to 8-bit would not give comparable
results. After investigation, we noticed it was mainly due to two reasons, i)
the cell states $\mathbf{c}_{t}$ had large ranges (e.g., $[-17,15]$), and ii)
the attention mechanism was not letting the decoder attend the encoder outputs
accurately. Therefore, we quantize the pre-activation gates (12), the element-
wise multiplications in (13) and cell states $\mathbf{c}_{t}$ to 16-bit. The
attention is quantized following our described integer-only attention method.
Everything else is quantized to 8-bit following our described method. The
quantized model has a similar performance to the full-precision model, with a
maximum of $1.25$ WER% drop (Table 3). We believe allowing the model to train
longer would reduce the gap.

###  5.4 Inference measurements

LSTM | quantized | forward pass (ms) | iteration | speedup
---|---|---|---|---
| activation |  | in second |
Full-precision | X | 130 | 7.6 | $1.00\times$
iRNN w/ | PWL32 | 84 | 11.8 | $1.54\times$
iRNN | PWL8 | 61 | 14.9 | $\mathbf{1.95\times}$
iRNN | X | 127 | 7.8 | $1.02\times$
Table 4: Inference measurements on a Huawei P20 Pro smartphone, based on a
custom fork from PyTorch 1.7.1. The model is one LSTM cell with a state size
of 400. Milliseconds (ms) and iterations per second (iter/s) are averaged
across 100 runs.

We implemented an 8-bit quantized integer-only LSTM with PWL model based on a
custom PyTorch paszke2019pytorch  fork from 1.7.1. We implemented an integer-
only PWL kernel using NEON intrinsics. We benchmark the models on a Huawei P20
Pro smartphone, using the $\mathrm{speed\\_benchmark\\_torch}$
tool444https://github.com/pytorch/pytorch/blob/1.7/binaries/speed_benchmark_torch.cc.
We warm up each model for five runs, measure the inference time a hundred
times, and report an average. The sequence length used is 128, and the batch
size is one. We benchmark our iRNN LSTM model using PWLs with 32 pieces and
eight pieces, which achieve up to $2\times$ speedup. We also evaluate our iRNN
with full-precision computations (iRNN w/o QAct) for the activation where no
speedup was observed for this state size. We believe it is due to round-trip
conversions between floating points and integers (Table 4). There is much room
for improvements to achieve even greater speedup, such as writing a C++
integer-only LSTM cell, fusing operations, and better PWL kernel
implementation.

Unimplemented Operations | val | test
---|---|---
Quantized Weights (Full-precision) | $60.27\pm 0.34$ | $58.02\pm 0.34$
PWL | $60.40\pm 0.05$ | $57.90\pm 0.01$
Quantized Activations | $60.40\pm 0.03$ | $57.95\pm 0.003$
Quantized Element-wise ops | $60.08\pm 0.10$ | $57.61\pm 0.23$
Quantized Matmul | $60.10\pm 0.05$ | $57.64\pm 0.10$
Table 5: Ablation study on quantized Mogrifier LSTM training on WikiText2. iRNN without PWL is the quantized model using LUT instead of PWL to compute the activation function. Best results are averaged across three runs, and standard deviations are reported. Full-precision model | val | test
---|---|---
LayerNorm LSTM | $98.58\pm 0.35$ | $94.84\pm 0.21$
MadNorm LSTM | $\mathbf{97.20}\pm 0.47$ | $\mathbf{93.63}\pm 0.74$
Table 6: Word-level perplexities on PTB for a full-precision LSTM with
LayerNorm and a full-precision model with MadNorm. Best results are averaged
across three runs, and standard deviations are reported.

##  6 Conclusion

Speech and text applications are increasingly being deployed in IoT edge
applications (i.e. ASR in autonomous vehicles or mobile devices, neural
machine translation on portable translator devices, or even in many of today’s
smart home appliances). As RNN models form the basis of most these
applications, it is vital that efficient, hardware-friendly models are
designed to ensure usability on the edge. This requires that models use low-
precision computation which can improve device latency, power usage and
memory.

In this text, we propose a task-agnostic and flexible methodology which
enables integer-only RNNs that can tackle this problem. We are the first to
offer an approach to quantize all existing operations in modern RNNs,
supporting normalization and attention with a little performance drop. We
evaluated our approach on high-performance LSTM-based models for language
modeling and ASR, which have distinct architectures and variable computational
requirements. Our experiments showed promising results on popular
architectures and datasets namely; language modeling on WikiText2 with
Mogrifier LSTM where our integer-only Mogrifier LSTM achieved similar
performance to the full-precision model with only 0.3 perplexity increase with
a PWL of 32 pieces, as well as ASR on LibriSpeech with ESPRESSO LSTM where our
quantized solution only saw a $1.25$ WER% drop compared to the full-precision
model. These experiments sucessfully show that RNNs can be fully quantized
while achieving similar accuracy as their full-precision counterpart.

To emphasize the benefit of our integer-only RNN models on end-user devices,
we benchmark our iRNN LSTM on a Huawei P20 Pro smartphone, where we obtain up
to $2\times$ inference speedup and $4\times$ memory reduction. This
implementation showed that even a simple integer-only model can see
performance benefit compared to a full-precision while maintaining good
accuracy.

The method proposed in this text, highlights that it is now possible to deploy
a wide range of accurate and high performing RNN-based applications on the
edge. This opens the door for more improved usability and accessibility to
advanced speech and text processing applications for the end-user.

## Appendix A

###  A.1 Details on LSTM-based models

For BiLSTM cells, nothing stated in section Integer-only LSTM network is
changed except that we enforce the forward LSTM hidden state
$\overrightarrow{\mathbf{h}}_{t}$ and the backward LSTM hidden state
$\overleftarrow{\mathbf{h}}_{t}$ to share the same quantization parameters so
that they can be concatenated as a vector. If the model has embedding layers,
they are quantized to 8-bit as we found they were not sensitive to
quantization. If the model has residual connections (e.g., between LSTM
cells), they are quantized to 8-bit integers. In encoder-decoder models, the
attention layers would be quantized using the method described in Section 4.3.
The last fully-connected layer weights of the model are 8-bit quantized to
allow for 8-bit matrix multiplication. We do not quantize the outputs and let
them remain 32-bit integers as often this is where it is considered that the
model has done its job and that some postprocessing is performed (e.g., beam
search).

###  A.2 Experimental details

| parameters | standard | weight | activation | Total
---|---|---|---|---|---
|  | training | quant. | quant. | GPU time
LM PTB | 4M | 1h | 13h | 10h | 24h
LM WikiText2 | 44M | 37h | 87h | 50h | 174h
ASR LibriSpeech | 174M | 144h | 36h | 12h | 192h
Table 7: Details about the model parameters and training time for each
different experiment setup: language modeling on PTB (LM PTB), language
modeling on WikiText2 (LM WikiText2), automatic speech recognition on
LibriSpeech (ASR LibriSpeech). Parameters are given in the order of millions.
We report the total number of V100 GPU hours as well as the GPU hours used for
each different training phase, i) pre-trained model (standard training), ii)
quantized model without PWL (weight quant.), iii) quantized model with PWL
activation quantization (activation quant.) and the total GPU time.

In this section, we provide further details of our experimental setups. The
number of parameters and training time of the model are reported in Table 7.

####  A.2.1 LayerNorm LSTM on PTB

Preprocessing of the dataset was performed mikolov2012subword . The vocabulary
size is 10K. We report the best _perplexity_ per word on the validation set
and test set for a language model of embedding size 200 with one LayerNormLSTM
cell of state size 200. The lower the perplexity, the better the model
performs. These experiments focus on the relative increase of perplexity
between the full-precision models and their 8-bit quantized counterparts. We
did not aim to reproduce state-of-the-art performance on PTB and went with a
naïve set of hyper-parameters. The full-precision network is trained for 100
epochs with batch size 20 and BPTT werbos1990bptt  window size of 35. We used
the SGD optimizer with weight decay of $10^{-5}$ and learning rate 20, which
is divided by 4 when the loss plateaus for more than two epochs without a
relative decrease of $10^{-4}$ in perplexity. We use gradient clipping of
0.25. We initialize the quantized models from the best full-precision
checkpoint and train for another 100 epochs. We did not enable quantization to
gather range statistics to compute the quantization parameters for the first
five epochs.

####  A.2.2 Mogrifier LSTM on WikiText2

We describe the experimental setup for Mogrifier LSTM on WikiText2. Note that
we follow the setup of melis2020mogrifier  where they do not use dynamic
evaluation krause2018dyneval  nor Monte Carlo dropout gal2016mcdropout . The
vocabulary size is 33279. We use a two-layer Mogrifier LSTM with embedding
dimension 272, state dimension 1366, and capped input gates. We use six
modulation rounds per Mogrifier layer with low-rank dimension 48. We use 2
Mixture-of-Softmax layers yang2018mos . The input and output embedding are
tied. We use a batch size of 64 and a BPTT window size of 70. We train the
full-precision Mogrifier LSTM for 340 epochs, after which we enable Stochastic
Weight Averaging (SWA) izmailov2018swa  for 70 epochs. For the optimizer we
used Adam kingma2014adam  with a learning rate of $\approx 3\times 10^{-3}$,
$\beta_{1}=0$, $\beta_{2}=0.999$ and weight decay $\approx 1.8\times 10^{-4}$.
We clip gradients’ norm to 10. We use the same hyper-parameters for the
quantized models from which we initialize with a pre-trained full-precision
and continue to train for 200 epochs. During the first two epochs, we do not
perform QAT, but we gather min and max statistics in the network to have a
correct starting estimate of the quantization parameters. After that, we
enable 8-bit QAT on every component of the Mogrifier LSTM: weights, matrix
multiplications, element-wise operations, and activations. Then we replace
activation functions in the model with quantization-aware PWLs and continue
training for 100 epochs. We perform complete ablation on our method to study
the effect of each component. Quantizing the weights and matrix
multiplications cover about $0.1$ of the perplexity increase. There is a clear
performance drop after adding quantization of element-wise operations with an
increase in the perplexity of about $0.3$. This is due to element-wise
operations in the cell and hidden state computations affecting the flow of
information across timesteps and the residual connections across layers.
Adding quantization of the activation does not impact the performance of the
network.

####  A.2.3 ESPRESSO LSTM on LibriSpeech

The encoder comprises 4 CNN-BatchNorm-ReLU blocks followed by 4 BiLSTM layers
with 1024 units. The decoder consists of 3 LSTM layers of units 1024, with
Bahdanau’s attention on hidden states of the encoder and residual connections
between each layer. The dataset preprocessing is precisely the same as in
wang2019espresso . We train the model for 30 epochs on one V100 GPU,
approximately six days to complete. We use a batch size of 24 while limiting
the maximum number of tokens in a mini-batch to 26000. Adam is used with a
starting learning rate of $0.001$, which is divided by 2 when the validation
set metric plateaus without a relative decrease of $10^{-4}$ in performance.
Cross-entropy with uniform label smoothing $\alpha=0.1$ szegedy2016inceptionv3
is used as a loss function. The model predictions are weighted at evaluation
time using a pre-trained full-precision 4-layer LSTM language model (shallow
fusion). We consider this language model as an external component to the
ESPRESSO LSTM; we do not quantize it due to the lack of resources. In our
language modeling experiments, we have already shown that quantized language
models retain their performance. We refer the reader to wang2019espresso  and
training
script555https://github.com/freewym/espresso/blob/master/examples/asr_librispeech/run.sh
for a complete description of the experimental setup.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 3: Quantization-aware PWL for several functions: exponential function
approximation over $x\in[-10,0]$ with a) 4, b) 8, c) 16, d) 32 pieces; cosine
over $x\in[-\pi,\pi]$ with e) 4, f) 8 pieces; GeLU hendrycks2016gelu  over
$x\in[-2,2]$ with g) 4, h) 5 pieces.

We initialize the quantized model from the pre-trained full-precision ESPRESSO
LSTM. Due to the lack of resources, we have trained the quantized model for
only four epochs. The quantized model is trained on 6 V100 GPUs where each
epoch takes two days, so a total of 48 GPU days. The batch size is set to 8
mini-batch per GPU with a maximum of 8600 tokens. We made these changes since,
otherwise, the GPU would run out of VRAM due to the added fake quantization
operations. For the first half of the first epoch, we gathered statistics for
quantization parameters then we enabled QAT. The activation functions are
swapped with quantization-aware PWL in the last epoch. The number of pieces
for the quantization-aware PWLs is 96, except for the exponential function in
the attention, which is 160, as we found out it was necessary to have more
pieces because of its curvature. The number of pieces used is higher than that
in the language modeling experiments. However, the difference is that the
inputs to the activation functions are 16-bit rather than 8-bit, although the
outputs are still quantized to 8-bit. It means we need more pieces to capture
the input resolution better. Note that it would not be feasible to use a
16-bit Look-Up Table to compute the activation functions due to the size and
cache misses, whereas using 96 pieces allows for a 170x reduction in memory
consumption compared to LUT.

Figure 4: Example of an iteration from our proposed quantization-aware PWL
Algorithm 1. The algorithm proceeds to reduce the number of pieces by merging
two similar adjacents pieces. In this figure, the slopes $S_{12}$ and $S_{23}$
are the most similar pieces; therefore, the knot $k_{2}$ is removed.

## Appendix B

The following section provides some examples of integer-only arithmetic and
more details on fixed-point scaling.

###  B.1 Multiplication

To illustrate how integer-only multiplication is achieved, we define an
example utilizing (4). Defining $u\in[u_{\min},u_{\max}]=[-1,1]$ and
$w\in[w_{\min},w_{\max}]=[0,5]$, the multiplication between two numbers from
those ranges will fall into $[z_{\min},z_{\max}]=[-5,5]$. From (3), for 8-bit
quantization, we have $S_{u}\approx 0.0078$, $Z_{u}=128$, $S_{w}\approx
0.0196,Z_{w}=0,S_{z}\approx 0.0392,Z_{z}=128$. Given $u=-0.8$ and $w=2.3$, we
have $q_{u}=25$ and $q_{w}=117$. Therefore, following (4),

$\displaystyle q_{z}$
$\displaystyle=\Big{\lfloor}\frac{S_{u}S_{w}}{S_{z}}\Big{(}q_{u}q_{w}-q_{u}Z_{w}-q_{w}Z_{u}+Z_{u}Z_{w}\Big{)}\Big{\rceil}+Z_{z}$
$\displaystyle=\Big{\lfloor}\frac{0.0078\times 0.0196}{0.0392}\Big{(}25\times
117-25\times 0-117\times 128+128\times 0\Big{)}\Big{\rceil}+128$
$\displaystyle=81.$

Using (2), the floating-point representation of $q_{z}$ is $r_{z}=-1.8424$
which is close to $uv=-1.8399$. Note that we lost precision at two levels, the
first time when quantizing $u$ and $v$, then the second time when quantizing
$z$, the multiplication output.

###  B.2 Addition

As mentioned in Section 2.1, addition with quantized numbers can take two
forms. The first form is when the two numbers to be added share the same
scaling factor and zero-point. For instance, given $x_{1}=-0.3,x_{2}=0.7$ from
$[-1,1]$, and $S_{x}=0.0078,Z_{x}=128$, we have $q_{x_{1}}=90$ and
$q_{x_{2}}=218$. The result value $y$ will fall into the range $[-2,2]$,
therefore $S_{y}\approx 0.0157$ and $Z_{y}=128$. Then, because they share the
same quantization parameters, following (5),

$\displaystyle q_{y}$
$\displaystyle=\Big{\lfloor}\frac{1}{S_{y}}\Big{(}S_{x}(q_{x_{1}}-Z_{x})+S_{x}(q_{x_{2}}-Z_{x})\Big{)}\Big{\rceil}+Z_{y}$
$\displaystyle=\Big{\lfloor}\frac{S_{x}}{S_{y}}\Big{(}q_{x_{1}}+q_{x_{2}}-2Z_{x}\Big{)}\Big{\rceil}+Z_{y}$
$\displaystyle=\Big{\lfloor}\frac{0.0078}{0.0157}\Big{(}90+218-256\Big{)}\Big{\rceil}+128$
$\displaystyle=154.$

We have $r_{y}=0.4082$, while $x_{1}+x_{2}=0.3999$. The second form is when
the two numbers do not share the same scaling factor and zero-point. Define
$a\in[a_{\min},a_{\max}]=[-1,1]$ and $b\in[b_{\min},b_{\max}]=[0,5]$, the
addition between two numbers from those ranges will fall into
$[c_{\min},c_{\max}]=[-1,6]$. We get $S_{a}\approx 0.0078$, $Z_{a}=128$,
$S_{b}\approx 0.0196,Z_{b}=0,S_{c}\approx 0.0274,Z_{c}=36$. For $a=-0.9$,
$b=3.9$, we have $q_{a}=13$ and $q_{b}=199$. The quantized addition result
$q_{c}$, following (6), is,

$\displaystyle q_{c}$
$\displaystyle=\Big{\lfloor}\frac{S_{a}}{S_{c}}(q_{a}-Z_{a})+\frac{S_{b}}{S_{c}}(q_{b}-Z_{b})\Big{\rceil}+Z_{c}$
$\displaystyle=\Big{\lfloor}\frac{0.0078}{0.0274}(13-128)+\frac{0.0196}{0.0274}199\Big{\rceil}+36$
$\displaystyle=146$

and $r_{y}=3.0140$ while $a+b=3.0$.

###  B.3 Fixed point arithmetic

Even with the most careful rounding, fixed-point values represented with a
scaling factor S may have an error of up to $\pm 0.5$ in the stored integer,
that is, $\pm 0.5S$ in the value. Therefore, smaller scaling factors generally
produce more accurate results. On the other hand, a smaller scaling factor
means a smaller range of values stored in a given program variable. The
maximum fixed-point value that can be stored in a variable is the largest
integer value that can be stored into it, multiplied by the scaling factor;
and similarly for the minimum value. For example, Table 1 gives the implied
scaling factor $S$, the minimum and maximum representable values. The accuracy
$\delta=S/2$ of values can be represented in 16-bit signed binary fixed-point
format, depending on the number f of implied fraction bits.

scaling | precision | signed (low, high) | unsigned (low, high)
---|---|---|---
$2^{1}$ | $2.0$ | $(-256,254)$ | $(0,510)$
$2^{0}$ | $1.0$ | $(-128,127)$ | $(0,255)$
$2^{-1}$ | $0.5$ | $(-64,63.5)$ | $(0,127.5)$
$2^{-2}$ | $0.25$ | $(-32,31.75)$ | $(0,63.75)$
$2^{-3}$ | $0.125$ | $(-16,15.875)$ | $(0,31.875)$
$2^{-4}$ | $0.0625$ | $(-8,7.9375)$ | $(0,15.9375)$
$2^{-5}$ | $0.03125$ | $(-4,3.96875)$ | $(0,7.96875)$
$2^{-6}$ | $0.015625$ | $(-2,1.984375)$ | $(0,3.984375)$
$2^{-7}$ | $0.0078125$ | $(-1,0.9921875)$ | $(0,1.9921875)$
$2^{-8}$ | $0.00390625$ | $(-0.5,0.49609375)$ | $(0,0.99609375)$
Table 1: Fixed point format with common signed scaling and bias format for 8
bit mantissa.

To convert a number from a floating-point to a fixed-point, one may divide it
by the scaling factor S, then round the result to the nearest integer. Care
must be taken to ensure that the result fits in the destination variable or
register. Depending on the scaling factor and storage size, and the range of
input numbers, the conversion may not entail any rounding. To convert a fixed-
point number to floating-point, in contrast, one may convert the integer to
floating-point and then multiply it by the scaling factor S. This conversion
may entail rounding if the integer’s absolute value is greater than 224 (for
binary single-precision IEEE floating point) or of 253 (for double-precision).
In addition, overflow or underflow may occur if $\left\lvert S\right\rvert$ is
very large or small. However, most computers with binary arithmetic have fast
bit shift instructions that can multiply or divide an integer by any power of
2, particularly an arithmetic shift instruction. These instructions can be
used to quickly change scaling factors that are powers of 2 while preserving
the sign of the number.

## References

  *   * (1) Rumelhart, D., Hinton, G.E., Williams, R.J.: Learning internal representations by error propagation. (1986)
  * (2) Hochreiter, S., Schmidhuber, J.: Long short-term memory. Neural Comput. 9(8), 1735–1780 (1997). https://doi.org/10.1162/neco.1997.9.8.1735
  * (3) Cho, K., Van Merriënboer, B., Bahdanau, D., Bengio, Y.: On the properties of neural machine translation: Encoder-decoder approaches. arXiv preprint arXiv:1409.1259 (2014)
  * (4) Chen, M.X., Firat, O., Bapna, A., Johnson, M., Macherey, W., Foster, G., Jones, L., Schuster, M., Shazeer, N., Parmar, N., Vaswani, A., Uszkoreit, J., Kaiser, L., Chen, Z., Wu, Y., Hughes, M.: The best of both worlds: Combining recent advances in neural machine translation. In: Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), pp. 76–86. Association for Computational Linguistics, Melbourne, Australia (2018). https://doi.org/10.18653/v1/P18-1008. https://www.aclweb.org/anthology/P18-1008
  * (5) Wang, C., Wu, S., Liu, S.: Accelerating transformer decoding via a hybrid of self-attention and recurrent neural network. arXiv preprint arXiv:1909.02279 (2019)
  * (6) Zhang, L., Wang, S., Liu, B.: Deep learning for sentiment analysis : A survey. CoRR abs/1801.07883 (2018) 1801.07883
  * (7) You, Q., Jin, H., Wang, Z., Fang, C., Luo, J.: Image captioning with semantic attention. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR) (2016)
  * (8) Mao, J., Xu, W., Yang, Y., Wang, J., Huang, Z., Yuille, A.: Deep captioning with multimodal recurrent neural networks (m-rnn) (2014). https://doi.org/10.48550/ARXIV.1412.6632
  * (9) He, Y., Sainath, T.N., Prabhavalkar, R., McGraw, I., Alvarez, R., Zhao, D., Rybach, D., Kannan, A., Wu, Y., Pang, R., et al.: Streaming end-to-end speech recognition for mobile devices. In: ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 6381–6385 (2019). IEEE
  * (10) Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A.N., Kaiser, L.u., Polosukhin, I.: Attention is all you need. In: Guyon, I., Luxburg, U.V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., Garnett, R. (eds.) Advances in Neural Information Processing Systems 30, pp. 5998–6008. Curran Associates, Inc., ??? (2017). http://papers.nips.cc/paper/7181-attention-is-all-you-need.pdf
  * (11) Devlin, J., Chang, M.-W., Lee, K., Toutanova, K.: BERT: Pre-training of deep bidirectional transformers for language understanding. In: Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers), pp. 4171–4186. Association for Computational Linguistics, Minneapolis, Minnesota (2019). https://doi.org/10.18653/v1/N19-1423. https://www.aclweb.org/anthology/N19-1423
  * (12) Kasai, J., Peng, H., Zhang, Y., Yogatama, D., Ilharco, G., Pappas, N., Mao, Y., Chen, W., Smith, N.A.: Finetuning pretrained transformers into rnns. arXiv preprint arXiv:2103.13076 (2021)
  * (13) Jacob, B., Kligys, S., Chen, B., Zhu, M., Tang, M., Howard, A.G., Adam, H., Kalenichenko, D.: Quantization and training of neural networks for efficient integer-arithmetic-only inference. 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2704–2713 (2017)
  * (14) Krishnamoorthi, R.: Quantizing deep convolutional networks for efficient inference: A whitepaper. arXiv preprint arXiv:1806.08342 (2018)
  * (15) Courville, V., Nia, V.P.: Deep learning inference frameworks for arm cpu. Journal of Computational Vision and Imaging Systems 5(1), 3–3 (2019)
  * (16) Hubara, I., Courbariaux, M., Soudry, D., El-Yaniv, R., Bengio, Y.: Quantized neural networks: Training neural networks with low precision weights and activations. Journal of Machine Learning Research 18(187), 1–30 (2018)
  * (17) Darabi, S., Belbahri, M., Courbariaux, M., Nia, V.P.: BNN+: improved binary network training. CoRR abs/1812.11800 (2018) arXiv:1812.11800
  * (18) Esser, S.K., McKinstry, J.L., Bablani, D., Appuswamy, R., Modha, D.S.: Learned step size quantization. arXiv preprint arXiv:1902.08153 (2019)
  * (19) Ott, J., Lin, Z., Zhang, Y., Liu, S.-C., Bengio, Y.: Recurrent neural networks with limited numerical precision. arXiv preprint arXiv:1608.06902 (2016)
  * (20) He, Q., Wen, H., Zhou, S., Wu, Y., Yao, C., Zhou, X., Zou, Y.: Effective quantization methods for recurrent neural networks. arXiv preprint arXiv:1611.10176 (2016)
  * (21) Kapur, S., Mishra, A., Marr, D.: Low precision rnns: Quantizing rnns without losing accuracy. arXiv preprint arXiv:1710.07706 (2017)
  * (22) Hou, L., Zhu, J., Kwok, J.T.-Y., Gao, F., Qin, T., Liu, T.-y.: Normalization helps training of quantized lstm (2019)
  * (23) Ardakani, A., Ji, Z., Smithson, S.C., Meyer, B.H., Gross, W.J.: Learning recurrent binary/ternary weights. arXiv preprint arXiv:1809.11086 (2018)
  * (24) Sari, E., Partovi Nia, V.: Batch normalization in quantized networks. In: Proceedings of the Edge Intelligence Workshop, pp. 6–9 (2020). https://www.gerad.ca/en/papers/G-2020-23
  * (25) Wu, Y., Schuster, M., Chen, Z., Le, Q.V., Norouzi, M., Macherey, W., Krikun, M., Cao, Y., Gao, Q., Macherey, K., et al.: Google’s neural machine translation system: Bridging the gap between human and machine translation. arXiv preprint arXiv:1609.08144 (2016)
  * (26) Bluche, T., Primet, M., Gisselbrecht, T.: Small-footprint open-vocabulary keyword spotting with quantized lstm networks. arXiv preprint arXiv:2002.10851 (2020)
  * (27) Ba, J., Kiros, J.R., Hinton, G.E.: Layer normalization. ArXiv abs/1607.06450 (2016)
  * (28) Bahdanau, D., Cho, K., Bengio, Y.: Neural machine translation by jointly learning to align and translate. In: Bengio, Y., LeCun, Y. (eds.) 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings (2015). http://arxiv.org/abs/1409.0473
  * (29) Chorowski, J., Bahdanau, D., Serdyuk, D., Cho, K., Bengio, Y.: Attention-based models for speech recognition. arXiv preprint arXiv:1506.07503 (2015)
  * (30) Sari., E., Courville., V., Partovi Nia., V.: iRNN: Integer-only Recurrent Neural Network. In: Proceedings of the 11th International Conference on Pattern Recognition Applications and Methods - ICPRAM, pp. 110–121 (2022)
  * (31) Marcus, M.P., Santorini, B., Marcinkiewicz, M.A.: Building a large annotated corpus of English: The Penn Treebank. Computational Linguistics 19(2), 313–330 (1993)
  * (32) Merity, S., Xiong, C., Bradbury, J., Socher, R.: Pointer sentinel mixture models. CoRR abs/1609.07843 (2016) arXiv:1609.07843
  * (33) Melis, G., Kočiský, T., Blunsom, P.: Mogrifier lstm. In: International Conference on Learning Representations (2020). https://openreview.net/forum?id=SJe5P6EYvS
  * (34) Wang, Y., Chen, T., Xu, H., Ding, S., Lv, H., Shao, Y., Peng, N., Xie, L., Watanabe, S., Khudanpur, S.: Espresso: A fast end-to-end neural speech recognition toolkit. In: 2019 IEEE Automatic Speech Recognition and Understanding Workshop (ASRU), pp. 136–143 (2019). https://doi.org/10.1109/ASRU46091.2019.9003968
  * (35) Panayotov, V., Chen, G., Povey, D., Khudanpur, S.: Librispeech: an asr corpus based on public domain audio books. In: Acoustics, Speech and Signal Processing (ICASSP), 2015 IEEE International Conference On, pp. 5206–5210 (2015). IEEE
  * (36) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al.: Pytorch: An imperative style, high-performance deep learning library. Advances in neural information processing systems 32 (2019)
  * (37) Mikolov, T., Sutskever, I., Deoras, A., Le, H.-S., Kombrink, S., Cernocky, J.: Subword language modeling with neural networks. preprint (http://www. fit. vutbr. cz/imikolov/rnnlm/char. pdf) 8, 67 (2012)
  * (38) Werbos, P.J.: Backpropagation through time: what it does and how to do it. Proceedings of the IEEE 78(10), 1550–1560 (1990)
  * (39) Krause, B., Kahembwe, E., Murray, I., Renals, S.: Dynamic evaluation of neural sequence models. In: International Conference on Machine Learning, pp. 2766–2775 (2018). PMLR
  * (40) Gal, Y., Ghahramani, Z.: A theoretically grounded application of dropout in recurrent neural networks. Advances in neural information processing systems 29 (2016)
  * (41) Yang, Z., Dai, Z., Salakhutdinov, R., Cohen, W.W.: Breaking the Softmax Bottleneck: A High-Rank RNN Language Model (2018)
  * (42) Izmailov, P., Podoprikhin, D., Garipov, T., Vetrov, D., Wilson, A.G.: Averaging weights leads to wider optima and better generalization. arXiv preprint arXiv:1803.05407 (2018)
  * (43) Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. CoRR abs/1412.6980 (2014)
  * (44) Szegedy, C., Vanhoucke, V., Ioffe, S., Shlens, J., Wojna, Z.: Rethinking the inception architecture for computer vision. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 2818–2826 (2016)
  * (45) Hendrycks, D., Gimpel, K.: Gaussian error linear units (gelus). arXiv preprint arXiv:1606.08415 (2016)