Title: QUATERNION RECURRENT NEURAL NETWORKS

Abstract: Recurrent neural networks (RNNs) are powerful architectures to model sequential data, due to their capability to learn short and long-term dependencies between the basic elements of a sequence. Nonetheless, popular tasks such as speech or images recognition, involve multi-dimensional input features that are characterized by strong internal dependencies between the dimensions of the input vector. We propose a novel quaternion recurrent neural network (QRNN), alongside with a quaternion long-short term memory neural network (QLSTM), that take into account both the external relations and these internal structural dependencies with the quaternion algebra. Similarly to capsules, quaternions allow the QRNN to code internal dependencies by composing and processing multidimensional features as single entities, while the recurrent operation reveals correlations between the elements composing the sequence. We show that both QRNN and QLSTM achieve better performances than RNN and LSTM in a realistic application of automatic speech recognition. Finally, we show that QRNN and QLSTM reduce by a maximum factor of 3.3x the number of free parameters needed, compared to real-valued RNNs and LSTMs to reach better results, leading to a more compact representation of the relevant information.

1. INTRODUCTION:

In the last few years, deep neural networks (DNN) have encountered a wide success in different domains due to their capability to learn highly complex input to output mapping. Among the different DNN-based models, the recurrent neural network (RNN) is well adapted to process sequential data. Indeed, RNNs build a vector of activations at each timestep to code latent relations between input vectors. Deep RNNs have been recently used to obtain hidden representations of speech unit sequences (Ravanelli et al., 2018a) or text word sequences (Conneau et al., 2018) , and to achieve state-of-the-art performances in many speech recognition tasks (Graves et al., 2013a; b; Amodei et al., 2016; Povey et al., 2016; Chiu et al., 2018) . However, many recent tasks based on multi-dimensional input features, such as pixels of an image, acoustic features, or orientations of 3D models, require to represent both external dependencies between different entities, and internal relations between the features that compose each entity. Moreover, RNN-based algorithms commonly require a huge number of parameters to represent sequential data in the hidden space.

Quaternions are hypercomplex numbers that contain a real and three separate imaginary components, perfectly fitting to 3 and 4 dimensional feature vectors, such as for image processing and robot kinematics (Sangwine, 1996; Pei & Cheng, 1999; Aspragathos & Dimitros, 1998) . The idea of bundling groups of numbers into separate entities is also exploited by the recent manifold and capsule networks (Chakraborty et al., 2018; Sabour et al., 2017) . Contrary to traditional homogeneous representations, capsule and quaternion networks bundle sets of features together. Thereby, quaternion numbers allow neural network based models to code latent inter-dependencies between groups of input features during the learning process with fewer parameters than RNNs, by taking advantage of the Hamilton product as the equivalent of the ordinary product, but between quaternions. Early applications of quaternion-valued backpropagation algorithms (Arena et al., 1994; 1997) have efficiently solved quaternion functions approximation tasks. More recently, neural networks of complex and hypercomplex numbers have received an increasing attention (Hirose & Yoshida, 2012; Tygert et al., 2016; Danihelka et al., 2016; Wisdom et al., 2016) , and some efforts have shown promising results in different applications. In particular, a deep quaternion network (Parcollet et al., 2016; 2017a; b) , a deep quaternion convolutional network (Gaudet & Maida, 2018; Parcollet et al., 2018) , or a deep complex convolutional network (Trabelsi et al., 2017) have been employed for challenging tasks such as images and language processing. However, these applications do not include recurrent neural networks with operations defined by the quaternion algebra.

This paper proposes to integrate local spectral features in a novel model called quaternion recurrent neural network1 (QRNN), and its gated extension called quaternion long-short term memory neural network (QLSTM). The model is proposed along with a well-adapted parameters initialization and turned out to learn both inter-and intra-dependencies between multidimensional input features and the basic elements of a sequence with drastically fewer parameters (Section 3), making the approach more suitable for low-resource applications. The effectiveness of the proposed QRNN and QLSTM is evaluated on the realistic TIMIT phoneme recognition task (Section 4.2) that shows that both QRNN and QLSTM obtain better performances than RNNs and LSTMs with a best observed phoneme error rate (PER) of 18.5% and 15.1% for QRNN and QLSTM, compared to 19.0% and 15.3% for RNN and LSTM. Moreover, these results are obtained alongside with a reduction of 3.3 times of the number of free parameters. Similar results are observed with the larger Wall Street Journal (WSJ) dataset, whose detailed performances are reported in the Appendix 6.1.1.

2. MOTIVATIONS:

A major challenge of current machine learning models is to well-represent in the latent space the astonishing amount of data available for recent tasks. For this purpose, a good model has to efficiently encode local relations within the input features, such as between the Red, Green, and Blue (R,G,B) channels of a single image pixel, as well as structural relations, such as those describing edges or shapes composed by groups of pixels. Moreover, in order to learn an adequate representation with the available set of training data and to avoid overfitting, it is convenient to conceive a neural architecture with the smallest number of parameters to be estimated. In the following, we detail the motivations to employ a quaternion-valued RNN instead of a real-valued one to code inter and intra features dependencies with fewer parameters.

As a first step, a better representation of multidimensional data has to be explored to naturally capture internal relations within the input features. For example, an efficient way to represent the information composing an image is to consider each pixel as being a whole entity of three strongly related elements, instead of a group of uni-dimensional elements that could be related to each other, as in traditional real-valued neural networks. Indeed, with a real-valued RNN, the latent relations between the RGB components of a given pixel are hardly coded in the latent space since the weight has to find out these relations among all the pixels composing the image. This problem is effectively solved by replacing real numbers with quaternion numbers. Indeed, quaternions are fourth dimensional and allow one to build and process entities made of up to four related features. The quaternion algebra and more precisely the Hamilton product allows quaternion neural network to capture these internal latent relations within the features encoded in a quaternion. It has been shown that QNNs are able to restore the spatial relations within 3D coordinates (Matsui et al., 2004) , and within color pixels (Isokawa et al., 2003) , while real-valued NN failed. This is easily explained by the fact that the quaternion-weight components are shared through multiple quaternion-input parts during the Hamilton product , creating relations within the elements. Indeed, Figure 1 shows that the multiple weights required to code latent relations within a feature are considered at the same level as for learning global relations between different features, while the quaternion weight w codes these internal relations within a unique quaternion Q out during the Hamilton product (right).

Then, while bigger neural networks allow better performances, quaternion neural networks make it possible to deal with the same signal dimension but with four times less neural parameters. Indeed, a 4-number quaternion weight linking two 4-number quaternion units only has 4 degrees of freedom, whereas a standard neural net parametrization has 4 × 4 = 16, i.e., a 4-fold saving in memory. Therefore, the natural multidimensional representation of quaternions alongside with their ability to drastically reduce the number of parameters indicate that hyper-complex numbers are a better fit than real numbers to create more efficient models in multidimensional spaces. Based on the success of previous deep quaternion convolutional neural networks and smaller quaternion feed-forward architectures (Kusamichi et al., 2004; Isokawa et al., 2009; Parcollet et al., 2017a) , this work proposes to adapt the representation of hyper-complex numbers to the capability of recurrent neural networks in a natural and efficient framework to multidimensional sequential tasks such as speech recognition.

Modern automatic speech recognition systems usually employ input sequences composed of multidimensional acoustic features, such as log Mel features, that are often enriched with their first, second and third time derivatives (Davis & Mermelstein, 1990; Furui, 1986) , to integrate contextual information. In standard RNNs, static features are simply concatenated with their derivatives to form a large input vector, without effectively considering that signal derivatives represent different views of the same input. Nonetheless, it is crucial to consider that time derivatives of the spectral energy in a given frequency band at a specific time frame represent a special state of a time-frame, and are linearly correlated (Tokuda et al., 2003) . Based on the above motivations and the results observed on previous works about quaternion neural networks, we hypothesize that quaternion RNNs naturally provide a more suitable representation of the input sequence, since these multiple views can be directly embedded in the multiple dimensions space of the quaternion, leading to better generalization.

3. QUATERNION RECURRENT NEURAL NETWORKS:

This Section describes the quaternion algebra (Section 3.1), the internal quaternion representation (Section 3.2), the backpropagation through time (BPTT) for quaternions (Section 3.3.2), and proposes an adapted weight initialization to quaternion-valued neurons (Section 3.4).

3.1. QUATERNION ALGEBRA:

The quaternion algebra H defines operations between quaternion numbers. A quaternion Q is an extension of a complex number defined in a four dimensional space as:

EQUATION

where r, x, y, and z are real numbers, and 1, i, j, and k are the quaternion unit basis. In a quaternion, r is the real part, while xi + yj + zk with i 2 = j 2 = k 2 = ijk = -1 is the imaginary part, or the vector part. Such a definition can be used to describe spatial rotations. The information embedded in the quaterion Q can be summarized into the following matrix of real numbers:

EQUATION

The conjugate Q * of Q is defined as:

Q * = r1 -xi -yj -zk.

(3) Then, a normalized or unit quaternion Q is expressed as:

EQUATION

)

Finally, the Hamilton product ⊗ between two quaternions Q 1 and Q 2 is computed as follows:

Q 1 ⊗ Q 2 =(r 1 r 2 -x 1 x 2 -y 1 y 2 -z 1 z 2 ) + (r 1 x 2 + x 1 r 2 + y 1 z 2 -z 1 y 2 )i+ (r 1 y 2 -x 1 z 2 + y 1 r 2 + z 1 x 2 )j + (r 1 z 2 + x 1 y 2 -y 1 x 2 + z 1 r 2 )k.

(5) The Hamilton product (a graphical view is depicted in Figure 1 ) is used in QRNNs to perform transformations of vectors representing quaternions, as well as scaling and interpolation between two rotations following a geodesic over a sphere in the R 3 space as shown in (Minemoto et al., 2017) .

3.2. QUATERNION REPRESENTATION:

The QRNN is an extension of the real-valued (Medsker & Jain, 2001) and complex-valued (Hu & Wang, 2012; Song & Yam, 1998) recurrent neural networks to hypercomplex numbers. In a quaternion dense layer, all parameters are quaternions, including inputs, outputs, weights, and biases. The quaternion algebra is ensured by manipulating matrices of real numbers (Gaudet & Maida, 2018) . Consequently, for each input vector of size N , output vector of size M , dimensions are split into four parts: the first one equals to r, the second is xi, the third one equals to yj, and the last one to zk to compose a quaternion Q = r1 + xi + yj + zk. The inference process of a fully-connected layer is defined in the real-valued space by the dot product between an input vector and a real-valued M × N weight matrix. In a QRNN, this operation is replaced with the Hamilton product (Eq. 5) with quaternion-valued matrices (i.e. each entry in the weight matrix is a quaternion). The computational complexity of quaternion-valued models is discussed in Appendix 6.1.2

3.3. LEARNING ALGORITHM:

The QRNN differs from the real-valued RNN in each learning sub-processes. Therefore, let x t be the input vector at timestep t, h t the hidden state, W hx , W hy and W hh the input, output and hidden states weight matrices respectively. The vector b h is the bias of the hidden state and p t , y t are the output and the expected target vectors. More details of the learning process and the parametrization are available on Appendix 6.2.

3.3.1. FORWARD PHASE:

Based on the forward propagation of the real-valued RNN (Medsker & Jain, 2001) , the QRNN forward equations are extended as follows:

EQUATION

where α is a quaternion split activation function (Xu et al., 2017; Tripathi, 2016) defined as:

EQUATION

with f corresponding to any standard activation function. The split approach is preferred in this work due to better prior investigations, better stability (i.e. pure quaternion activation functions contain singularities), and simpler computations. The output vector p t is computed as:

EQUATION

where β is any split activation function. Finally, the objective function is a classical loss applied component-wise (e.g., mean squared error, negative log-likelihood).

3.3.2. QUATERNION BACKPROPAGATION THROUGH TIME:

The backpropagation through time (BPTT) for quaternion numbers (QBPTT) is an extension of the standard quaternion backpropagation (Nitta, 1995) , and its full derivation is available in Appendix 6.3. The gradient with respect to the loss E t is expressed for each weight matrix as ∆ t hy = ∂Et ∂W hy , ∆ t hh = ∂Et ∂W hh , ∆ t hx = ∂Et ∂W hx , for the bias vector as ∆ t b = ∂Et ∂B h , and is generalized to ∆ t = ∂Et ∂W with:

EQUATION

)

Each term of the above relation is then computed by applying the chain rule. Indeed, and conversaly to real-valued backpropagation, QBPTT must defines the dynamic of the loss w.r.t to each component of the quaternion neural parameters. As a use-case for the equations, the mean squared error at a timestep t and named E t is used as the loss function. Moreover, let λ be a fixed learning rate. First, the weight matrix W hy is only seen in the equations of p t . It is therefore straightforward to update each weight of W hy at timestep t following:

EQUATION

where h * t is the conjugate of h t . Then, the weight matrices W hh , W hx and biases b h are arguments of h t with h t-1 involved, and the update equations are derived as:

EQUATION

with,

EQUATION

and,

EQUATION

with h preact n and p preact n the pre-activation values of h n and p n respectively.

3.4. PARAMETER INITIALIZATION:

A well-designed parameter initialization scheme strongly impacts the efficiency of a DNN. An appropriate initialization, in fact, improves DNN convergence, reduces the risk of exploding or vanishing gradient, and often leads to a substantial performance improvement (Glorot & Bengio, 2010) . It has been shown that the backpropagation through time algorithm of RNNs is degraded by an inappropriated parameter initialization (Sutskever et al., 2013) . Moreover, an hyper-complex parameter cannot be simply initialized randomly and component-wise, due to the interactions between components. Therefore, this Section proposes a procedure reported in Algorithm 1 to initialize a matrix W of quaternion-valued weights. The proposed initialization equations are derived from the polar form of a weight w of W :

EQUATION

and, w r = ϕ cos(θ), w i = ϕ q imagi sin(θ), w j = ϕ q imagj sin(θ), w k = ϕ q imagk sin(θ). (15)

The angle θ is randomly generated in the interval [-π, π]. The quaternion q imag is defined as purely normalized imaginary, and is expressed as q imag = 0 + xi + yj + zk. The imaginary components xi, yj, and zk are sampled from an uniform distribution in [0, 1] to obtain q imag , which is then normalized (following Eq. 4) to obtain q imag . The parameter ϕ is a random number generated with respect to well-known initialization criterions (such as Glorot or He algorithms) (Glorot & Bengio, 2010; He et al., 2015) . However, the equations derived in (Glorot & Bengio, 2010; He et al., 2015) are defined for real-valued weight matrices. Therefore, the variance of W has to be investigated in the quaternion space to obtain ϕ (the full demonstration is provided in Appendix 6.2). The variance of W is:

EQUATION

)

Algorithm 1 Quaternion-valued weight initialization

1: procedure QINIT(W, n in , n out ) 2: σ ← 1 √ 2(nin+nout)

w.r.t to Glorot criterion and Eq. 18

3:

for w in W do 4:

θ ← rand(-π, π)

5:

ϕ ← rand(-σ, σ)

6:

x, y, z ← rand(0, 1)

7:

q imag ← Quaternion(0, x, y, z)

8:

q imag ← qimag √ x 2 +y 2 +z 2 9:

w r ← ϕ × cos(θ) See Eq. 15 10:

w i ← ϕ × q imagi × sin(θ)

11:

w j ← ϕ × q imagj × sin(θ) 12: w k ← ϕ × q imag k × sin(θ) 13:

w ← Quaternion(w r , w i , w j , w k )

Indeed, the weight distribution is normalized. The value of V ar(W ) = E(|W |2 ), instead, is not trivial in the case of quaternion-valued matrices. Indeed, W follows a Chi-distribution with four degrees of freedom (DOFs). Consequently, V ar(W ) is expressed and computed as follows:

EQUATION

)

The Glorot (Glorot & Bengio, 2010) and He (He et al., 2015) criterions are extended to quaternion as:

EQUATION

with n in and n out the number of neurons of the input and output layers respectively. Finally, ϕ can be sampled from [-σ, σ] to complete the weight initialization of Eq. 15.

4. EXPERIMENTS:

This Section details the acoustic features extraction (Section 4.1), the experimental setups and the results obtained with QRNNs, QLSTMs, RNNs and LSTMs on the TIMIT speech recognition tasks (Section 4.2). The results reported in bold on tables are obtained with the best configurations of the neural networks observed with the validation set.

4.1. QUATERNION ACOUSTIC FEATURES:

The raw audio is first splitted every 10ms with a window of 25ms. Then 40-dimensional log Mel-filterbank coefficients with first, second, and third order derivatives are extracted using the pytorch-kaldi 2 (Ravanelli et al., 2018b) toolkit and the Kaldi s5 recipes (Povey et al., 2011 ). An acoustic quaternion Q(f, t) associated with a frequency f and a time-frame t is formed as follows:

EQUATION

Q(f, t) represents multiple views of a frequency f at time frame t, consisting of the energy e(f, t) in the filter band at frequency f , its first time derivative describing a slope view, its second time derivative describing a concavity view, and the third derivative describing the rate of change of the second derivative. Quaternions are used to learn the spatial relations that exist between the 3 described different views that characterize a same frequency (Tokuda et al., 2003) . Thus, the quaternion input vector length is 160/4 = 40. Decoding is based on Kaldi (Povey et al., 2011) and weighted finite state transducers (WFST) (Mohri et al., 2002) that integrate acoustic, lexicon and language model probabilities into a single HMM-based search graph.

4.2. THE TIMIT CORPUS:

The training process is based on the standard 3, 696 sentences uttered by 462 speakers, while testing is conducted on 192 sentences uttered by 24 speakers of the TIMIT (Garofolo et al., 1993) dataset.

A validation set composed of 400 sentences uttered by 50 speakers is used for hyper-parameter tuning. The models are compared on a fixed number of layers M = 4 and by varying the number of neurons N from 256 to 2, 048, and 64 to 512 for the RNN and QRNN respectively. Indeed, it is worth underlying that the number of hidden neurons in the quaternion and real spaces do not handle the same amount of real-number values. Indeed, 256 quaternion neurons output are 256 × 4 = 1024 real values. Tanh activations are used across all the layers except for the output layer that is based on a softmax function. Models are optimized with RMSPROP with vanilla hyper-parameters and an initial learning rate of 8 • 10 -4 . The learning rate is progressively annealed using a halving factor of 0.5 that is applied when no performance improvement on the validation set is observed. The models are trained during 25 epochs. All the models converged to a minimum loss, due to the annealed learning rate. A dropout rate of 0.2 is applied over all the hidden layers (Srivastava et al., 2014) except the output one. The negative log-likelihood loss function is used as an objective function.

All the experiments are repeated 5 times (5-folds) with different seeds and are averaged to limit any variation due to the random initialization. The results on the TIMIT task are reported in Table 1 . The best PER in realistic conditions (w.r.t to the best validation PER) is 18.5% and 19.0% on the test set for QRNN and RNN models respectively, highlighting an absolute improvement of 0.5% obtained with QRNN. These results compare favorably with the best results obtained so far with architectures that do not integrate access control in multiple memory layers (Ravanelli et al., 2018a ). In the latter, a PER of 18.3% is reported on the TIMIT test set with batch-normalized RNNs . Moreover, a remarkable advantage of QRNNs is a drastic reduction (with a factor of 2.5×) of the parameters needed to achieve these results. Indeed, such PERs are obtained with models that employ the same internal dimensionality corresponding to 1, 024 real-valued neurons and 256 quaternion-valued ones, resulting in a number of parameters of 3.8M for QRNN against the 9.4M used in the real-valued RNN. It is also worth noting that QRNNs consistently need fewer parameters than equivalently sized RNNs, with an average reduction factor of 2.26 times. This is easily explained by considering the content of the quaternion algebra. Indeed, for a fully-connected layer with 2, 048 input values and 2, 048 hidden units, a real-valued RNN has 2, 048 2 ≈ 4.2M parameters, while to maintain equal input and output dimensions the quaternion equivalent has 512 quaternions inputs and 512 quaternion hidden units. Therefore, the number of parameters for the quaternion-valued model is 512 2 × 4 ≈ 1M. Such a complexity reduction turns out to produce better results and has other advantages such as a smaller memory footprint while saving models on budget memory systems. This characteristic makes our QRNN model particularly suitable for speech recognition conducted on low computational power devices like smartphones (Chen et al., 2014) . QRNNs and RNNs accuracies vary accordingly to the architecture with better PER on bigger and wider topologies. Therefore, while good PER are observed with a higher number of parameters, smaller architectures performed at 23.9% and 23.4%, with 1M and 0.6M parameters for the RNN and the QRNN respectively. Such PER are due to a too small number of parameters to solve the task.

4.3. QUATERNION LONG-SHORT TERM MEMORY NEURAL NETWORKS:

We propose to extend the QRNN to state-of-the-art models such as long-short term memory neural networks (LSTM), to support and improve the results already observed with the QRNN compared to the RNN in more realistic conditions. LSTM (Hochreiter & Schmidhuber, 1997) neural networks Published as a conference paper at ICLR 2019 were introduced to solve the problems of long-term dependencies learning and vanishing or exploding gradient observed with long sequences. Based on the equations of the forward propagation and back propagation through time of QRNN described in Section 3.3.1, and Section 3.3.2, one can easily derive the equations of a quaternion-valued LSTM. Gates are defined with quaternion numbers following the proposal of Danihelka et al. (2016) . Therefore, the gate action is characterized by an independent modification of each component of the quaternion-valued signal following a componentwise product with the quaternion-valued gate potential. Let f t ,i t , o t , c t , and h t be the forget, input, output gates, cell states and the hidden state of a LSTM cell at time-step t:

EQUATION

EQUATION

EQUATION

EQUATION

where W are rectangular input weight matrices, R are square recurrent weight matrices, and b are bias vectors. α is the split activation function and × denotes a component-wise product between two quaternions. Both QLSTM and LSTM are bidirectional and trained on the same conditions than for the QRNN and RNN experiments. The results on the TIMIT corpus reported on Table 2 support the initial intuitions and the previously established trends. We first point out that the best PER observed is 15.1% and 15.3% on the test set for QLSTMs and LSTM models respectively with an absolute improvement of 0.2% obtained with QLSTM using 3.3 times fewer parameters compared to LSTM. These results are among the top of the line results (Graves et al., 2013b; Ravanelli et al., 2018a) and prove that the proposed quaternion approach can be used in state-of-the-art models. A deeper investigation of QLSTMs performances with the larger Wall Street Journal (WSJ) dataset can be found in Appendix 6.1.1.

5. CONCLUSION:

Summary. This paper proposes to process sequences of multidimensional features (such as acoustic data) with a novel quaternion recurrent neural network (QRNN) and quaternion long-short term memory neural network (QLSTM). The experiments conducted on the TIMIT phoneme recognition task show that QRNNs and QLSTMs are more effective to learn a compact representation of multidimensional information by outperforming RNNs and LSTMs with 2 to 3 times less free parameters. Therefore, our initial intuition that the quaternion algebra offers a better and more compact representation for multidimensional features, alongside with a better learning capability of feature internal dependencies through the Hamilton product, have been demonstrated.

Future Work. Future investigations will develop other multi-view features that contribute to decrease ambiguities in representing phonemes in the quaternion space. In this extent, a recent approach based on a quaternion Fourier transform to create quaternion-valued signal has to be investigated. Finally, other high-dimensional neural networks such as manifold and Clifford networks remain mostly unexplored and can benefit from further research.

6. APPENDIX 6.1 WALL STREET JOURNAL EXPERIMENTS AND COMPUTATIONAL COMPLEXITY:

This Section proposes to validate the scaling of the proposed QLSTMs to a bigger and more realistic corpus, with a speech recognition task on the Wall Street Journal (WSJ) dataset. Finally, it discuses the impact of the quaternion algebra in term of computational compexity. It is important to notice that reported results on Table 3 compare favorably with equivalent architectures (Graves et al., 2013a) (WER of 11.7% on 'test-dev93'), and are competitive with state-of-the-art and much more complex models based on better engineered features (Chan & Lane, 2015) (WER of 3.8% with the 81 hours of training data, and on 'test-eval92'). According to Table 3 , QLSTMs outperform LSTM in all the training conditions (14 hours and 81 hours) and with respect to both the validation and testing sets. Moreover, QLSTMs still need 2.9 times less neural parameters than LSTMs to achieve such performances. This experiment demonstrates that QLSTMs scale well to larger and more realistic speech datasets and are still more efficient than real-valued LSTMs.

6.1.2. NOTES ON COMPUTATIONAL COMPLEXITY:

A computational complexity of O(n 2 ) with n the number of hidden states has been reported by Morchid (2018) for real-valued LSTMs. QLSTMs just involve 4 times larger matrices during computations. Therefore, the computational complexity remains unchanged and equals to O(n 2 ). Nonetheless, and due to the Hamilton product, a single forward propagation between two quaternion neurons uses 28 operations, compared to a single one for two real-valued neurons, implying a longer training time (up to 3 times slower). However, such worst speed performances could easily be alleviated with a proper engineered cuDNN kernel for the Hamilton product, that would helps QNNs to be more efficient than real-valued ones. A well-adapted CUDA kernel would allow QNNs to perform more computations, with fewer parameters, and therefore less memory copy operations from the CPU to the GPU.

6.2. PARAMETERS INITIALIZATION:

Let us recall that a generated quaternion weight w from a weight matrix W has a polar form defined as:

EQUATION

with q imag = 0 + xi + yj + zk a purely imaginary and normalized quaternion. Therefore, w can be computed following: w r = ϕ cos(θ), w i = ϕ q imagi sin(θ), w j = ϕ q imagj sin(θ), w k = ϕ q imagk sin(θ).

(26) However, ϕ represents a randomly generated variable with respect to the variance of the quaternion weight and the selected initialization criterion. The initialization process follows (Glorot & Bengio, 2010) and (He et al., 2015) to derive the variance of the quaternion-valued weight parameters. Indeed, the variance of W has to be investigated:

EQUATION

) [E(|W |)]

2 is equals to 0 since the weight distribution is symmetric around 0. Nonetheless, the value of V ar(W ) = E(|W | 2 ) is not trivial in the case of quaternion-valued matrices. Indeed, W follows a Chi-distribution with four degrees of freedom (DOFs) and E(|W | 2 ) is expressed and computed as follows:

EQUATION

With f (x) is the probability density function with four DOFs. A four-dimensional vector X = {A, B, C, D} is considered to evaluate the density function f (x). X has components that are normally distributed, centered at zero, and independent. Then, A, B, C and D have density functions:

EQUATION

)

The four-dimensional vector X has a length L defined as

L = √ A 2 + B 2 + C 2 + D 2

with a cumulative distribution function F L (x; σ) in the 4-sphere (n-sphere with n = 4) S x :

EQUATION

where

S x = {(a, b, c, d) : √ a 2 + b 2 + c 2 + d 2 <

x} and dS x = da db dc dd. The polar representations of the coordinates of X in a 4-dimensional space are defined to compute dS x :

a = ρ cos θ, b = ρ sin θ cos φ, c = ρ sin θ sin φ cos ψ, d = ρ sin θ sin φ sin ψ, where ρ is the magnitude (ρ = √ a 2 + b 2 + c 2 + d 2

) and θ, φ, and ψ are the phases with 0 ≤ θ ≤ π, 0 ≤ φ ≤ π and 0 ≤ ψ ≤ 2π. Then, dS x is evaluated with the Jacobian J f of f defined as: (31)

J f = ∂(a, b, c, d) ∂(ρ, θ, φ, ψ) = da

Therefore, by the Jacobian J f , we have the polar form:

da db dc dd = ρ 3 sin 2 θ sin φ dρ dθ dφ dψ.

(32)

Then, writing Eq.( 30) in polar coordinates, we obtain:

F L (x, σ) = 1 √ 2πσ 2 4 x 0 e -a 2 /2σ 2 e -b 2 /2σ 2 e -c 2 /2σ 2 e -d 2 /2σ 2 dS x = 1 4π 2 σ 4 2π 0 π 0 π 0 x 0 e -ρ 2 /2σ 2 ρ 3 sin 2 θ sin φ dρ dθ dφ dψ = 1 4π 2 σ 4 2π 0 dψ π 0 sin φ dφ π 0 sin 2 θ dθ x 0 ρ 3 e -ρ 2 /2σ 2 dρ = 1 4π 2 σ 4 2π2 θ 2 - sin 2θ 4 π 0 x 0 ρ 3 e -ρ 2 /2σ 2 dρ = 1 4π 2 σ 4 4π π 2 x 0 ρ 3 e -ρ 2 /2σ 2 dρ,

Then,

EQUATION

The probability density function for X is the derivative of its cumulative distribution function, which by the fundamental theorem of calculus is:

EQUATION

)

The expectation of the squared magnitude becomes:

E(|W | 2 ) = ∞ 0 x 2 f (x) dx = ∞ 0 x 2 1 2σ 4 x 3 e -x 2 /2σ 2 dx = 1 2σ 4 ∞ 0

x 5 e -x 2 /2σ 2 dx.

With integration by parts we obtain:

EQUATION

)

The expectation E(|W | 2 ) is the sum of two terms. The first one:

-x 4 e -x 2 /2σ 2 ∞ 0 = lim x→+∞ -x 4 e -x 2 /2σ 2 -lim x→+0 x 4 e -x 2 /2σ 2 = lim x→+∞ -x 4 e -x 2 /2σ 2 ,

Based on the L'Hôpital's rule, the undetermined limit becomes:

lim x→+∞ -x 4 e -x 2 /2σ 2 = -lim x→+∞ x 4 e x 2 /2σ 2 = . . . = -lim x→+∞ 24

(1/σ 2 )(P (x)e x 2 /2σ 2 ) (36) = 0.

With P (x) is polynomial and has a limit to +∞. The second term is calculated in a same way (integration by parts) and E(|W | 2 ) becomes from Eq.( 35):

E(|W | 2 ) = 1 2σ 2 ∞ 0 4x 3 e -x 2 /2σ 2 dx = 2 σ 2 x 2 σ 2 e -x 2 /2σ 2 ∞ 0 + ∞ 0 σ 2 2xe -x 2 /2σ 2 dx .

(37) The limit of first term is equals to 0 with the same method than in Eq.( 36). Therefore, the expectation is:

E(|W | 2 ) = 4 ∞ 0 xe -x 2 /2σ 2 dx = 4σ 2 .

(38) And finally the variance is:

V ar(|W |) = 4σ 2 .

(39)

Let us recall the forward equations and parameters needed to derive the complete quaternion backpropagation through time (QBPTT) algorithm.

6.3.1. RECALL OF THE FORWARD PHASE:

Let x t be the input vector at timestep t, h t the hidden state, W hh , W xh and W hy the hidden state, input and output weight matrices respectively. Finally b h is the biases vector of the hidden states and p t , y t are the output and the expected target vector.

EQUATION

) and α is the quaternion split activation function (Xu et al., 2017) of a quaternion Q defined as:

EQUATION

) and f corresponding to any standard activation function. The output vector p t can be computed as:

EQUATION

) and β any split activation function. Finally, the objective function is a real-valued loss function applied component-wise. The gradient with respect to the MSE loss is expressed for each weight matrix as ∂Et ∂W hy , ∂Et ∂W hh , ∂Et ∂W hx , and for the bias vector as ∂Et ∂B h . In the real-valued space, the dynamic of the loss is only investigated based on all previously connected neurons. In this extent, the QBPTT differs from BPTT due to the fact that the loss must also be derived with respect to each component of a quaternion neural parameter, making it bi-level. This could act as a regularizer during the training process. 

6.3.2. OUTPUT WEIGHT MATRIX:

EQUATION

By regrouping in a matrix form the h t components from these equations, one can define:

EQUATION

Therefore,

EQUATION

6.3.3. HIDDEN WEIGHT MATRIX:

Conversely to W hy the weight matrix W hh is an argument of h t with h t-1 involved. The recursive backpropagation can thus be derived as:

EQUATION

And,

EQUATION

with N the number of timesteps that compose the sequence. As for W hy we start with (64)

The same operations are performed for i,j,k in Eq. 68 and ∂Et ∂W hh can finally be expressed as: 

EQUATION

Nonetheless, since biases are not connected to any inputs or hidden states, the matrix of derivatives defined in Eq. 59 becomes a matrix of 1. Consequently ∂Et ∂B h can be summarized as:

EQUATION

:

https://github.com/Orkis-Research/Pytorch-Quaternion-Neural-Networks

pytorch-kaldi is available at https://github.com/mravanelli/pytorch-kaldi

Figures and Tables:

Figure 1: Illustration of the input features (Q in ) latent relations learning ability of a quaternion-valued layer (right) due to the quaternion weight sharing of the Hamilton product (Eq. 5), compared to a standard real-valued layer (left).

φ ρ sin θ cos φ -ρ sin θ sin φ 0 sin θ sin φ cos ψ ρ cos θ sin φ cos ψ ρ sin θ cos φ cos ψ -ρ sin θ sin φ sin ψ sin θ sin φ sin ψ ρ cos θ sin φ sin ψ ρ sin θ cos φ sin ψ ρ sin θ sin φ cos ψ . And, J f = ρ 3 sin 2 θ sin φ.

Phoneme error rate (PER%) of QRNN and RNN models on the development and test sets of the TIMIT dataset. "Params" stands for the total number of trainable parameters.

Table content: <table><tr><td>Models</td><td>Neurons</td><td>Dev.</td><td>Test</td><td>Params</td></tr><tr><td/><td>256</td><td>22.4</td><td>23.4</td><td>1M</td></tr><tr><td>RNN</td><td>512 1,024</td><td>19.6 17.9</td><td>20.4 19.0</td><td>2.8M 9.4M</td></tr><tr><td/><td>2,048</td><td>20.0</td><td>20.7</td><td>33.4M</td></tr><tr><td/><td>64</td><td>23.6</td><td>23.9</td><td>0.6M</td></tr><tr><td>QRNN</td><td>128 256</td><td>19.2 17.4</td><td>20.1 18.5</td><td>1.4M 3.8M</td></tr><tr><td/><td>512</td><td>17.5</td><td>18.7</td><td>11.2M</td></tr></table>

Phoneme error rate (PER%) of QLSTM and LSTM models on the development and test sets of the TIMIT dataset. "Params" stands for the total number of trainable parameters.

Table content: <table><tr><td>Models</td><td>Neurons</td><td>Dev.</td><td>Test</td><td>Params</td></tr><tr><td/><td>256</td><td>14.9</td><td>16.5</td><td>3.6M</td></tr><tr><td>LSTM</td><td>512 1,024</td><td>14.2 14.4</td><td>16.1 15.3</td><td>12.6M 46.2M</td></tr><tr><td/><td>2,048</td><td>14.0</td><td>15.9</td><td>176.3M</td></tr><tr><td/><td>64</td><td>15.5</td><td>17.0</td><td>1.6M</td></tr><tr><td>QLSTM</td><td>128 256</td><td>14.1 14.0</td><td>16.0 15.1</td><td>4.6M 14.4M</td></tr><tr><td/><td>512</td><td>14.2</td><td>15.1</td><td>49.9M</td></tr></table>

6.1.1 SPEECH RECOGNITION WITH THE WALL STREET JOURNAL CORPUSWe propose to evaluate both QLSTMs and LSTMs with a larger and more realistic corpus to validate the scaling of the observed TIMIT results (Section 4.2). Acoustic input features are described in Section 4.1, and extracted on both the 14 hour subset 'train-si84', and the full 81 hour dataset 'train-si284' of the Wall Street Journal (WSJ) corpus. The 'test-dev93' development set is employed for validation, while 'test-eval92' composes the testing set. Models architectures are fixed with respect to the best results observed with the TIMIT corpus (Section 4.2). Therefore, both QLSTMs and LSTMs contain four bidirectional layers of internal dimension of size 1, 024. Then, an additional layer of internal size 1, 024 is added before the output layer. The only change on the training procedure compared to the TIMIT experiments concerns the model optimizer, which is set to Adam(Kingma & Ba, 2014) instead of RMSPROP. Results are from a 3-folds average.

Table content: <table/>

Word error rates (WER %) obtained with both training set (WSJ14h and WSJ81h) of the Wall Street Journal corpus. 'test-dev93' and 'test-eval92' are used as validation and testing set respectively. L expresses the number of recurrent layers.

Table content: <table><tr><td>Models</td><td colspan="4">WSJ14 Dev. WSJ14 Test WSJ81 Dev. WSJ81 Test</td><td>Params</td></tr><tr><td>LSTM</td><td>11.2</td><td>7.2</td><td>7.4</td><td>4.5</td><td>53.7M</td></tr><tr><td>QLSTM</td><td>10.9</td><td>6.9</td><td>7.2</td><td>4.3</td><td>18.7M</td></tr></table>

The weight matrix W hy is used only in the computation of p t . It is therefore straightforward to compute ∂Et ∂W hy :

Table content: <table><tr><td/><td/><td/><td/><td/><td colspan="3">∂E t ∂W hy</td><td>=</td><td/><td colspan="2">∂E t ∂W r hy</td><td colspan="2">+ i</td><td>∂E t ∂W i hy</td><td>+ j</td><td>∂E t ∂W j hy</td><td>+ k</td><td>∂E t hy ∂W k</td><td>.</td><td>(45)</td></tr><tr><td colspan="15">Each quaternion component is then derived following the chain rule:</td></tr><tr><td colspan="2">∂E t ∂W r hy</td><td>=</td><td colspan="2">∂E t ∂p r t</td><td colspan="3">∂p r t ∂W r hy</td><td>+</td><td colspan="2">∂E t ∂p i t</td><td colspan="3">∂p i t ∂W r hy</td><td>+</td><td>∂E t ∂p j t</td><td>∂p j t ∂W r hy</td><td>+</td><td>∂E t ∂p k t</td><td>∂p k t ∂W r hy</td><td>(46)</td></tr><tr><td/><td/><td colspan="13">= (p r t -y r t ) × h r t + (p i t -y i t ) × h i t + (p j t -y j t ) × h j t + (p k t -y k t ) × h k t .</td></tr><tr><td>∂E t ∂W i hy</td><td>=</td><td colspan="2">∂E t ∂p r t</td><td colspan="2">∂p r t ∂W i hy</td><td>+</td><td colspan="3">∂E t ∂p i t</td><td colspan="2">∂p i t ∂W i hy</td><td>+</td><td colspan="2">∂E t ∂p j t</td><td>∂p j t ∂W i hy</td><td>+</td><td>∂E t ∂p k t</td><td>∂p k t ∂W i hy</td><td>(47)</td></tr><tr><td/><td colspan="14">= (p r t -y r t ) × -h i t + (p i t -y i t ) × h r t + (p j t -y j t ) × h k t + (p k t -y k t ) × -h j t .</td></tr><tr><td>∂E t ∂W j hy</td><td>=</td><td colspan="2">∂E t ∂p r t</td><td colspan="2">∂p r t ∂W j hy</td><td>+</td><td colspan="3">∂E t ∂p i t</td><td colspan="2">∂p i t ∂W j hy</td><td>+</td><td colspan="2">∂E t ∂p j t</td><td>∂p j t ∂W j hy</td><td>+</td><td>∂E t ∂p k t</td><td>∂p k t ∂W j hy</td><td>(48)</td></tr><tr><td/><td colspan="14">= (p r t -y r t ) × -h j t + (p i t -y i t ) × -h k t + (p j t -y j t ) × h r t + (p k t -y k t ) × h i t .</td></tr><tr><td>∂E t ∂W k hy</td><td>=</td><td colspan="2">∂E t ∂p r t</td><td colspan="2">∂p r t ∂W k hy</td><td>+</td><td colspan="3">∂E t ∂p i t</td><td colspan="2">∂p i t ∂W k hy</td><td>+</td><td colspan="2">∂E t ∂p j t</td><td>∂p j t ∂W k hy</td><td>+</td><td>∂E t ∂p k t</td><td>∂p k t ∂W k hy</td></tr></table>

,∂W hx is computed in the exact same manner as ∂Et ∂W hh .

Table content: <table><tr><td>And,</td><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td colspan="2">∂E t ∂B h</td><td colspan="2">=</td><td colspan="2">t m=0</td><td colspan="2">∂E m ∂B r h</td><td colspan="2">+ i</td><td colspan="2">∂E m ∂B r h</td><td>+ j</td><td>∂E m ∂B i h</td><td>+ k</td><td>∂E m h ∂B k</td><td>.</td></tr><tr><td>δ n =</td><td/><td/><td colspan="10">W  *  hh ⊗ δ n+1 × α (h preact n W  *  hy ⊗ (p n -y n ) × β (p preact ) n</td><td>if n = t ) else.</td><td>(66)</td></tr><tr><td colspan="4">6.3.4 INPUT WEIGHT MATRIX</td><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td/><td/><td/><td/><td/><td/><td colspan="3">∂E ∂W hx</td><td colspan="2">=</td><td>N t=0</td><td>∂E t ∂W hx</td><td>.</td><td>(67)</td></tr><tr><td>And,</td><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td>∂E t ∂W hx</td><td colspan="2">=</td><td colspan="2">t m=0</td><td colspan="3">∂E m ∂W r hx</td><td colspan="2">+ i</td><td colspan="3">∂E m ∂W r hx</td><td>+ j</td><td>∂E m ∂W i hx</td><td>+ k</td><td>∂E m hx ∂W k</td><td>.</td><td>(68)</td></tr><tr><td colspan="5">Therefore ∂Et ∂W hx is easily extent as:</td><td/><td/><td/><td/><td/><td/><td/></tr><tr><td/><td/><td/><td/><td colspan="3">∂E t ∂W hx</td><td>=</td><td colspan="5">t m=0 ( n=m+1 t</td><td>δ n ) ⊗ x  *  t .</td><td>(69)</td></tr><tr><td>6.3.5 HIDDEN BIASES</td><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td/><td/><td/><td/><td/><td/><td/><td colspan="2">∂E ∂B h</td><td colspan="2">=</td><td>N t=0</td><td>∂E t ∂B h</td><td>.</td><td>(70)</td></tr></table>

Peer review (iclr2019):
 After the discussion with authors, I am happy to recommend acceptance. ———————————————————— 1. In “Consequently, for each input vector of size N, output vector of size M, dimensions are split into four parts: the first one equals to r, the second is xi, the third one equals to yj, and the last one to zk to compose a quaternion Q = r1 + xi + yj + zk”, are you splitting dimension M or M.times N? And if you split M .times N (I believe that’s what you are doing), in which order you are splitting (row major right?) Please explain. 2. I did not understand why authors didn’t go in the negative direction of the gradient in Eq. (10-11)? 3. In section 3.4, authors mentioned “Moreover, an hyper-complex parameter cannot be simply initialized randomly and component-wise, due to the interactions between components.” which I strongly agree. But in Eq. (7) and (9) why the update rules and activation function are applied component wise? 4. I really like the elegance in the parameter initialization. Couple of minor things here: (1) It’s better to mention in Eq. (16) why E(|W|) is 0 because of symmetry. (2) Reference should be 6.1 instead of 5.1. 5. Another reasonable baseline will be using a complex network like (https://openreview.net/forum?id=H1T2hmZAb) and use the first two terms in Eq. (19) for representation. This will also possibly justify the usefulness of using higher order partials. 6. The authors mentioned multiple times about the achieved state-of-the-art results without giving any citation. As a reader not well versed in the acoustic domain, it will be nice to see some references to cross-validate the claim made. General Comments: 1. I understand the necessity of defining RNN/ LSTM model in the space of quaternions. But unit quaternions can be identified with other spaces where convolution is defined recently, e.g., with S^3 (https://arxiv.org/abs/1809.06211). I can see that this paper is contemporary, but at least can authors comment on the applicability of this general method in their case? Given that in NIPS’18 the following paper talked about RNN model on non-Euclidean spaces (https://arxiv.org/pdf/1805.11204.pdf), one can extend these ideas to develop an RNN model in the space of quaternions. Authors should look into it rigorously as future directions? But at least please comments on the applicability. 2. The experimental results section is somewhat weak, the overall claim of using fewer parameters and achieving comparable results is only validated on TIMIT data. More experimentation is necessary. 3. In terms of technical novelty, though quaternion algebra is well-known, I like the parameter initialization algorithm. I can see the merit of this in ML/ vision community. Pros: 1. Nice well grounded methodological development on well-known algebra. (simple but elegant, so that*s good). 2. Nicely written and all the maths check out (that*s good). 3. Experimental result on TIMIT dataset shows usefulness in terms of using fewer parameters (but still can achieve SOA results). Cons: 1. See my comments above. I expect the authors to rebut/ address the aforementioned comments. Overall though simple but nice (and necessary) development of RNN/ LSTM framework in the space of quaternions. 2. Lacks extensive experimental validation. My reason for my rating is mainly because of (1) lack of experimental validation. (2) being aware of the recent development of general RNN model on non-Euclidean spaces, I want some comments in this direction (see detailed comment and reference above).