Orthogonal Gated Recurrent Unit with Neumann-Cayley Transformation

Edison Mucllari∗ Vasily Zadorozhnyy∗ Cole Pospisil Duc Nguyen Qiang Ye

Department of Mathematics, University of Kentucky, Lexington, Kentucky, USA
{edison.mucllari, vasily.zadorozhnyy, Cole.Pospisil, ducnguyen<EMAIL_ADDRESS>
###### Abstract

In recent years, using orthogonal matrices has been shown to be a promising
approach in improving Recurrent Neural Networks (RNNs) with training,
stability, and convergence, particularly, to control gradients. While Gated
Recurrent Unit (GRU) and Long Short Term Memory (LSTM) architectures address
the vanishing gradient problem by using a variety of gates and memory cells,
they are still prone to the exploding gradient problem. In this work, we
analyze the gradients in GRU and propose the usage of orthogonal matrices to
prevent exploding gradient problems and enhance long-term memory. We study
where to use orthogonal matrices and we propose a Neumann series-based Scaled
Cayley transformation for training orthogonal matrices in GRU, which we call
Neumann-Cayley Orthogonal GRU, or simply NC-GRU. We present detailed
experiments of our model on several synthetic and real-world tasks, which show
that NC-GRU significantly outperforms GRU as well as several other RNNs.

##  1 Introduction

One of the preferred neural network models for working with sequential data is
the Recurrent Neural Network (RNN) (Rumelhart et al.,, 1986; Hopfield,, 1982).
RNNs can efficiently model sequential data through a hidden sequence of
states. However, training vanilla RNNs have obstacles (Rumelhart et al.,,
1986; Hopfield,, 1982), one of which is their susceptibility to vanishing and
exploding gradients (Bengio et al.,, 1993). In the case of vanishing
gradients, the optimization algorithm faces difficulty continuing to learn due
to very small changes in the model parameters. In the case of exploding
gradients, the training could overstep local minima, potentially causing
instabilities such as divergence or oscillatory behavior. It may also lead to
computational overflows.

There have been several works studying how to solve these problems. For
example gates have been introduced into the RNN architecture: Long Short-Term
Memory (LSTM) (Hochreiter and Schmidhuber, 1997a, ) and Gated Recurrent Units
(GRU) (Cho et al.,, 2014). They can pass long term information and help to
overcome vanishing gradients. In practice, GRU and LSTM models are still prone
to the problem of exploding gradients.

More recently, several RNN models have been proposed using unitary or
orthogonal matrices for the recurrent weights (Helfrich et al.,, 2018; Dorado-
Rojas et al.,, 2020; Vorontsov et al., 2017a, ; Mhammedi et al.,, 2017; Jing
et al.,, 2017; Arjovsky et al.,, 2016; Wisdom et al.,, 2016; Jing et al.,,
2019; Maduranga et al.,, 2018) along with methods to preserve those
properties. The introduction of such weights into RNN models brought new
development into the RNN community. One of the main reasons behind this is a
theoretical explanation of why the performance is better when using unitary or
orthogonal weights (Arjovsky et al.,, 2016). The key step in all of these
methods is preserving orthogonal or unitary properties at every iteration of
training. There have been several different techniques for updating the
recurrent weights to preserve either orthogonal or unitary properties,
including for example multiplicative updates (Wisdom et al.,, 2016), Givens
rotations (Jing et al.,, 2017), Householder reflections (Mhammedi et al.,,
2017), Cayley transforms (Helfrich et al.,, 2018; Maduranga et al.,, 2018;
Helfrich and Ye,, 2019; Lezcano-Casado and Martínez-Rubio,, 2019), and other
similar ideas that have shown effective usage of orthogonal or unitary
matrices (Saxe et al.,, 2014; Arjovsky et al.,, 2016; Henaff et al.,, 2017;
Hyland and Gunnar,, 2017; Tagare,, 2011; Vorontsov et al., 2017b, ).

In this work, we study the benefits of applying orthogonal matrices to one of
the most widely used RNN models, Gated Recurrent Unit (GRU) (Cho et al.,,
2014), both theoretically and experimentally. We analyze the gradients of the
GRU loss and based on this analysis, we propose the usage of orthogonal
matrices in several hidden state weights of the model. We introduce a Neumann
series-based Scaled Cayley transforms for training the orthogonal weights. Our
method utilizes a reliable and stable method of Scaled Cayley transforms which
was studied and used in (Helfrich et al.,, 2018; Maduranga et al.,, 2018) for
training orthogonal weights for RNNs. In addition, we propose a Neumann series
approximation of matrix inverse inside the Cayley transform. Such substitution
not only performs either on the same or even better level as the traditional
inverse (see Section 4 for experiments), but also decreases computation time
which might come of particular assistance when working with larger models. We
call our method _Neumann-Cayley Orthogonal Gated Recurrent Unit_ or _NC-GRU_
for simplicity. Our experiments show that the proposed method is more stable
with faster convergence and produces better results on several synthetic and
real-world tasks.

The rest of the paper has the following layout. In Section 2, we discuss some
previous works that are most relevant for this paper. In Section 3, we
introduce Neumann-Cayley Orthogonal Gated Recurrent Unit (NC-GRU) method
together with backpropagation method and gradient analysis of GRU (Cho et
al.,, 2014) and NC-GRU. In Section 4, we present experiments of our model on
several synthetic and real-world problems: Adding task (Hochreiter and
Schmidhuber, 1997b, ), Copying task (Hochreiter and Schmidhuber, 1997b, ),
Parenthesis task (Jing et al.,, 2019; Foerster et al.,, 2016), and Denoise
task (Jing et al.,, 2019) as well as character and word level Penn TreeBank
prediction tasks (Marcus et al., 1993b, ). Then, in Section 5, we provide
several experiments to justify the use of orthogonal matrices, and Neumann
series as well as which hidden states benefit from them. Finally, in Section
6, we summarize our proposed method and contribution.

##  2 Related Work

Many models have been designed to improve classical RNN (Rumelhart et al.,,
1986; Hopfield,, 1982) for sequential data. They include, for example, the
establishment of gates (Hochreiter and Schmidhuber, 1997a, ; Cho et al.,,
2014), normalization methods (Ioffe and Szegedy,, 2015; Cooijmans et al.,,
2017; Wu and He,, 2018; Ulyanov et al.,, 2017; Salimans and Kingma,, 2016; Ba
et al.,, 2016; Xu et al.,, 2019), and the introduction of unitary and
orthogonal matrices into RNN structure (Arjovsky et al.,, 2016; Jing et al.,,
2017; Mhammedi et al.,, 2017; Vorontsov et al., 2017b, ; Jing et al.,, 2019;
Dorado-Rojas et al.,, 2020), etc. In this section, we will discuss some of the
works most relevant to our proposed method.

Unitary RNNs (uRNNs) (Arjovsky et al.,, 2016) presented an architecture that
learns a unitary weight matrix. The construction of the recurrent weight
matrix consists of a composition of diagonal matrices, reflection matrices in
the complex domain, and Fourier transform. The uRNN model presented in (Wisdom
et al.,, 2016) is based on the constrained optimization over the space of all
unitary matrices rather than a product of parameterized matrices. EUNN (Jing
et al.,, 2017) is another work that utilizes the product of unitary matrices.
The recurrent matrix in this architecture is parameterized with products of
Givens Rotations. Also, the representation capacity of the unitary space is
fully tunable and ranges from a subspace of unitary matrices to the entire
unitary space.

Work by (Mhammedi et al.,, 2017) proposed orthogonal RNNs, or simply oRNNs,
which involve the application of Householder reflections. Such
parameterization of the transition matrix allows efficient training and
maintains the orthogonality of the recurrent weights while training. The
method introduced in (Vorontsov et al., 2017b, ), proposes a weight matrix
factorization by bounding the matrix norms. Moreover, it controls the degree
of gradient expansion during backpropagation. Besides that, this technique
enforces orthogonality as well. GORU (Jing et al.,, 2019) presented a RNN
which is an extension of unitary RNN with a gating mechanism. Unitarity is
preserved by using the ideas from EURNN (Jing et al.,, 2017). The results
compared to GRU (Cho et al.,, 2014) are mixed, depending on the task. More
recently, (Dorado-Rojas et al.,, 2020) presented an embedding of a linear time
- invariant system that contains Laguerre polynomials in the model.

###  2.1 Gated Recurrent Unit (GRU)

In our work, we have studied in depth the Gated Recurrent Unit (GRU)
architecture which was proposed in (Cho et al.,, 2014) as an alternative to
the well-known LSTM (Hochreiter and Schmidhuber, 1997a, ) cell. The structure
of a GRU cell is

$\displaystyle r_{t}$
$\displaystyle=\sigma\left(W_{r}x_{t}+U_{r}h_{t-1}+b_{r}\right)$ (1)
$\displaystyle u_{t}$
$\displaystyle=\sigma\left(W_{u}x_{t}+U_{u}h_{t-1}+b_{u}\right)$
$\displaystyle c_{t}$
$\displaystyle=\Phi\left(W_{c}x_{t}+U_{c}\left(r_{t}\odot
h_{t-1}\right)+b_{c}\right)$ $\displaystyle h_{t}$
$\displaystyle=\left(1-u_{t}\right)\odot h_{t-1}+u_{t}\odot c_{t}$

where $W_{r}$, $W_{u}$, and $W_{c}$ are input weights in $\mathbb{R}^{n\times
m}$, $U_{r}$, $U_{u}$, and $U_{c}$ are recurrent weights in
$\mathbb{R}^{n\times n}$, and $b_{r}$, $b_{u}$ and $b_{c}$ are the bias
parameters in $\mathbb{R}^{n}$. Here, $m$ represents the dimension of the
input data and $n$ represents the dimension of the hidden state. In (1), the
activation functions $\sigma$ and $\Phi$ are sigmoid and hyperbolic tangent
function ($\tanh$) respectively, and $\odot$ is the Hadamard product. In
addition, the initial hidden state $h_{0}$ is set to zero.

The main difference of GRU from LSTM is the implementation of the long-term
memory not as a separate channel but inside the hidden state $h_{t}$ itself.
GRU has a single gate $u_{t}$ that controls both forget and input gates. For
example, if the output of $u_{t}$ is $1$ then the forget gate is open implying
that the input gate is closed. Similarly, if $u_{t}$ is $0$, the forget gate
is closed and input gate is open. This structure allows GRUs to discard random
or insignificant information but at the same time grasping the important
details.

###  2.2 Cayley Transform Orthogonal RNN

We have implemented and studied the effects of orthogonal matrices inside the
GRU cell based on the Cayley transforms (Tagare,, 2011). Some initial work to
use Cayley transform for orthogonal weights in RNNs was presented in (Helfrich
et al.,, 2018) together with scoRNN model. This model introduced a skew-
symmetric matrix $A$ which is used to define an orthogonal matrix $W$ via the
Scaled Cayley transform

$W=\left(I+A\right)^{-1}\left(I-A\right)D,$ (2)

where matrix $D$ is a diagonal matrix of ones and negative ones, which scales
the traditional Cayley transform (Tagare,, 2011). It is proved in (Kahan,,
2006) that the matrix $D$, with a suitable choice on the number of negative
ones, can avoid a potential problem of eigenvalue(s) of $A$ being negative
one, making matrix $I+A$ non-invertible. The number of negative ones in matrix
$D$ can be considered as a tunable hyperparameter. Further, it guarantees that
the skew-symmetric matrix $A$ that generates the orthogonal matrix will be
bounded.

(Helfrich et al.,, 2018) presents the following process to train scoRNN model
using Scaled Cayley transforms:

$\displaystyle A^{(k+1)}$
$\displaystyle=A^{(k)}-\lambda\nabla_{A}L\left(U_{sco}\left(A^{(k)}\right)\right)$
(3) $\displaystyle U_{sco}^{(k+1)}$
$\displaystyle=\left(I+A^{(k+1)}\right)^{-1}\left(I-A^{(k+1)}\right)D$ (4)

where $\nabla_{A}L\left(U_{sco}\left(A\right)\right)$ is computed using

$\nabla_{A}L\left(U_{sco}\left(A\right)\right)=V^{T}-V,\,$ (5)

with

$V=\left(I+A\right)^{-T}\nabla_{U_{sco}}L\left(U_{sco}\left(A\right)\right)\left(D+U_{sco}^{T}\right)$
(6)

in which $\nabla_{U_{sco}}L\left(U_{sco}\left(A\right)\right)$ is computed via
standard backpropogation methods.

Despite the fact that rounding errors may accumulate over a number of repeated
matrix multiplications, orthogonality in scoRNN (Helfrich et al.,, 2018) is
maintained to the machine precision. This property helps to achieve
significant improvements over other orthogonal/unitary RNNs for long sequences
on several benchmark tasks, see the Experiments section in (Helfrich et al.,,
2018) for more details.

##  3 Efficient Orthogonal Gated Recurrent Unit

We now present an efficient orthogonal GRU. The proofs of all theoretical
results presented in this section are given in section 8.

###  3.1 Gradient Analysis of Hidden States in GRU

Gradient behavior plays a very important role in model training, convergence,
stability, and most importantly performance. However when it comes to
backpropagation through time for the GRU model from (1), the gradients of the
loss function $\mathcal{L}$ with respect to intermediate hidden states,
weights, and biases can be found from the respective gradients of the final
hidden state $h_{T}$, which is simplified to finding the gradient of $h_{t}$
with respect to $h_{t-1}$ for $t$ between 1 and $T$. Namely,

$\dfrac{\partial{\cal L}}{\partial h_{i}}=\dfrac{\partial{\cal L}}{\partial
h_{T}}\prod_{t=i+1}^{T}\dfrac{\partial h_{t}}{\partial h_{t-1}}.$ (7)

Thus, to analyze the gradients, we consider the gradient of the hidden state
$h_{t}$ with respect to the hidden state $h_{t-1}$ as well as its upper bound
in the following theorem.

######  Theorem 1.

Let $h_{t-1}$ and $h_{t}$ be two consecutive hidden states from the GRU model
stated in (1). Then

$\left\lvert\left\lvert\dfrac{\partial h_{t}}{\partial
h_{t-1}}\right\rvert\right\rvert_{2}\leq\alpha+\beta\left\lvert\left\lvert
U_{c}\right\rvert\right\rvert_{2}$ (8)

where

$\displaystyle\alpha$
$\displaystyle=\delta_{u}\left(\max_{i}\left\\{[h_{t-1}]_{i}\right\\}+\max_{i}\left\\{[c_{t}]_{i}\right\\}\right)\left\lvert\left\lvert
U_{u}\right\rvert\right\rvert_{2}$ (9)
$\displaystyle\quad\quad\quad+\max_{i}\left\\{(1-[u_{t}]_{i})\right\\}$

and

$\displaystyle\beta$
$\displaystyle=\max_{i}\left\\{[u_{t}]_{i}\right\\}\left(\delta_{r}\left\lvert\left\lvert
U_{r}\right\rvert\right\rvert_{2}\max_{i}\left\\{[h_{t-1}]_{i}\right\\}\right.$
(10)
$\displaystyle\quad\quad\quad\left.+\max_{i}\left\\{[r_{t}]_{i}\right\\}\right),$

with constants $\delta_{u}$ and $\delta_{r}$ defined as follows:

$\delta_{u}=\max_{i}\left\\{\left[u_{t}\right]_{i}\left(1-\left[u_{t}\right]_{i}\right)\right\\}$
(11)

and

$\delta_{r}=\max_{i}\left\\{\left[r_{t}\right]_{i}\left(1-\left[r_{t}\right]_{i}\right)\right\\}.$
(12)

The following corollary provides some simple upper bounds for $\alpha,\beta$
obtained in the above theorem.

######  Corollary 1.

For the hyperbolic tangent activation function in (1) (i.e. $\Phi=\,$`tanh`),
we have $\delta_{u},\delta_{r}\leq\frac{1}{4}$, $[h_{t}]_{i}\leq 1$ for any
$i$ and $t$ as well as

$\alpha\leq\dfrac{1}{2}\left\lvert\left\lvert
U_{u}\right\rvert\right\rvert_{2}+1\quad\text{and}\quad\beta\leq\dfrac{1}{4}\left\lvert\left\lvert
U_{r}\right\rvert\right\rvert_{2}+1.$ (13)

These bounds may be pessimistic because the gate elements may be expected to
be close to 0 or 1. As a consequence, the following corollary presents the
relationship between the constants $\alpha$ and $\beta$ presented in Theorem 1
when GRU’s gates are nearly closed or opened. Below we use a notation
$x\lesssim y$ to denote that $x$ is bounded by a quantity that is
approximately equal to $y$.

######  Corollary 2.

When elements of GRU gates $u_{t}$ and $r_{t}$ are nearly either 0 or 1, then
constants $\alpha$ and $\beta$ from Theorem 1 satisfy the following
inequality:

$\alpha+\beta\lesssim 2.$ (14)

Moreover if $u_{t}$ and $r_{t}$ are nearly either the zero vector or the
vector of all ones, then

$\alpha+\beta\lesssim 1.$ (15)

###  3.2 Neumann-Cayley Orthogonal Transformation

Based on Theorem 1 and Corollary 2, we propose the usage of orthogonal weights
in the hidden parameters of GRU to obtain a better conditioned gradients. As
discussed in section 2, there have been different techniques proposed and used
to initialize and preserve orthogonal weights while training, for example
Givens rotations (Jing et al.,, 2017), Householder reflections (Mhammedi et
al.,, 2017) etc. In this work, we implement a version of the Scaled-Cayley
transformation method discussed in 2.2 with one key difference. The Scaled
Cayley transform method requires calculation of the inverse of $I+A^{(k)}$ to
update the orthogonal matrix $U^{(k)}$ in (4). When it comes to the
computation of this inverse, classical numerical methods such as using LU-
decomposition or solving the Least Squares problem can be implemented, and
these methods work well when the dimension of the matrix is small. However, if
the dimension of the matrix is large then these methods are very expensive
from both memory and computational time perspectives. Moreover, classical
methods might overflow and not converge at all. We propose to solve this
possible complication by using the Neumann Series method to approximate the
inverse of $I+A^{(k)}$.

To derive the Neumann series approximation for the inverse of $I+A^{(k)}$ in
(4) we consider the following:

$\displaystyle\left(I+A^{(k)}\right)^{-1}$
$\displaystyle=\left(I+A^{(k-1)}-\delta A^{(k)}\right)^{-1}$ (16) (17) (18)

where $\delta
A^{(k)}:=optimizer_{A}\left(\nabla_{A}\mathcal{L}={V^{(k)}}^{T}-V^{(k)}\right)$,
here $optimizer_{A}$ includes a learning rate $\lambda$ inside of it. Note
that the equality between Equations (17) and (18) relies on the assumption
that $\left\lvert\left\lvert\left(I+A^{(k-1)}\right)^{-1}\delta
A^{(k)}\right\rvert\right\rvert<1$ for some operator norm
$\left\lvert\left\lvert\cdot\right\rvert\right\rvert$, see (Demmel,, 1997) for
more details. We have conducted an ablation study that shows empirical
evidence that this condition is indeed satisfied, see section 5.3 for more
details.

In our experiments, we have considered the first and the second-order Neumann
series approximations, estimating the series in (18) with two ($i=0,1$) and
three ($i=0,1,2$) terms respectively. As expected, the performance of the
model is slightly better when using second-order approximation although it
comes with a marginal increase in computational time, see 5.1 for the ablation
study regarding the accuracy and computational time of such approximations.
Mathematically speaking, if we are using the second-order approximation, the
error is of order $\mathcal{O}\left(\left(\left(I+A^{(k-1)}\right)^{-1}\delta
A^{(k)}\right)^{3}\right)$. Even though this error is quite small, there is a
chance that the errors from this approximation can accumulate and cause a loss
of orthogonality. To avoid this issue we recommend resetting orthogonality by
computing the matrix inverse explicitly using a factorization method at the
beginning of each epoch. However, it might be necessary to do it more often
particularly in the earlier training (e.g. every 100 iterations) due to more
fluctuations in the gradients.

1Given: $D$, $A^{(0)}$, $U^{(0)}$,
$\nabla_{U}\mathcal{L}\left(U^{(0)}\left(A^{(0)}\right)\right)$,
$optimizer_{A}$

2Define: $\tilde{A}^{(0)}:=\left(I+A^{(0)}\right)^{-1}$

3

4for  _$k=1,2,\ldots$_ do

5
$V^{(k)}:=\text{$\tilde{A}^{(k-1)}$}^{T}\nabla_{U}\mathcal{L}\left(U^{(k-1)}\left(A^{(k-1)}\right)\right)\left(D+{U^{(k-1)}}^{T}\right)$

6 $\delta
A^{(k)}:=optimizer_{A}\left(\nabla_{A}\mathcal{L}={V^{(k)}}^{T}-V^{(k)}\right)$

7 $A^{(k)}:=A^{(k-1)}-\delta A^{(k)}$

8 $\tilde{A}^{(k)}:=\left(I+\tilde{A}^{(k-1)}{\delta
A^{(k)}}+\left(\tilde{A}^{(k-1)}{\delta
A^{(k)}}\right)^{2}\right)\tilde{A}^{(k-1)}$

9 $U^{(k)}:=\tilde{A}^{(k)}\left(I-A^{(k)}\right)D$

10 end for

Algorithm 1 Update Rule for Orthogonal Weight $U$

We present Algorithm 1 that outlines the Neumann-Cayley Orthogonal
Transformation method for training weight $A$ as well as updating the
corresponding orthogonal weight $U$. It is important to note that during the
initialization step, the weight $A^{(0)}$ is defined to be skew-symmetric
using the same initialization technique as in (Helfrich et al.,, 2018) which
is based on the idea from (Henaff et al.,, 2017). Then, we apply the Cayley
transform in $A^{(0)}$ to obtain $U^{(0)}$. Another peculiar detail that we
want to point out is that Algorithm 1 includes $optimizer_{A}$, which is the
standard optimizer such as SGD, RMSProp (Tieleman and Hinton,, 2012), or Adam
(Kingma and Ba,, 2014), etc., that takes
$\nabla_{A}\mathcal{L}={V^{(k)}}^{T}-V^{(k)}$ as an input. Moreover, the skew-
symmetric property of the weight $A$ and its gradient are preserved under the
such optimizer.

###  3.3 Neumann-Cayley Orthogonal GRU (NC-GRU)

Finally, we introduce a Neumann-Cayley Orthogonal GRU (NC-GRU) model that
utilizes the proposed Neumann-Cayley Orthogonal Transform.

The structure of a NC-GRU cell is shown below:

$\displaystyle r_{t}$
$\displaystyle=\sigma\left(W_{r}x_{t}+U_{r}(A_{r})h_{t-1}+b_{r}\right)$ (19)
$\displaystyle u_{t}$
$\displaystyle=\sigma\left(W_{u}x_{t}+U_{u}h_{t-1}+b_{u}\right)$
$\displaystyle c_{t}$
$\displaystyle=\Phi\left(W_{c}x_{t}+U_{c}(A_{c})\left(r_{t}\odot
h_{t-1}\right)\right)$ $\displaystyle h_{t}$
$\displaystyle=\left(1-u_{t}\right)\odot h_{t-1}+u_{t}\odot c_{t}$

here $\sigma$ represents sigmoid function, $\odot$ \- Hadamard product, and
$\Phi$ \- modReLU function defined in (Arjovsky et al.,, 2016) as

$\Phi(x):=\text{modReLU}(x):=\text{sgn}(x)\cdot\text{ReLU}\left(\left\lvert
x\right\rvert+b\right)$ (20)

with $b$ as a trainable bias.

We derive from most of the experiments that the best performance is achieved
using orthogonality in $U_{c}$ and $U_{r}$ hidden weights. In addition, we
present an ablation study in 5.2 about the usage and performance of orthogonal
weights throughout the GRU model. Similar to the GRU Cell, $W_{r}$, $W_{u}$,
$W_{c}$, $U_{u}$, $b_{r}$, $b_{u}$, $b$, are trainable parameters as well as
$U_{r}$ and $U_{c}$ together with their associated weights $A_{r}$ and $A_{c}$
respectively. Moreover, all of them except $U_{r}$ with $A_{r}$ and $U_{c}$
with $A_{c}$ are trained using standard backpropagation algorithms such as
Stochastic Gradient Descent (SGD), RMSProp (Tieleman and Hinton,, 2012), or
Adam (Kingma and Ba,, 2014) in a similar manner as in GRU (Cho et al.,, 2014),
however weights $U_{r}$, $A_{r}$, $U_{c}$, and $A_{c}$ are trained using
Algorithm 1.

As we mentioned before, orthogonal weights lead to a better conditioned
gradient, the following Corollary summarizes this result.

######  Corollary 3.

Let $h_{t-1}$ and $h_{t}$ be two consecutive hidden states from the NC-GRU
model defined in (19). Then $\left\lvert\left\lvert
U_{r}\right\rvert\right\rvert_{2}=\left\lvert\left\lvert
U_{c}\right\rvert\right\rvert_{2}=1$ and if the element of the gates $u_{t}$
and $r_{t}$ are nearly either 0 or 1, then the following inequality is
satisfied:

$\left\lvert\left\lvert\dfrac{\partial h_{t}}{\partial
h_{t-1}}\right\rvert\right\rvert_{2}\lesssim 2.$ (21)

Furthermore, if $u_{t}$ and $r_{t}$ are nearly either zero vector or vector of
all ones,

$\left\lvert\left\lvert\dfrac{\partial h_{t}}{\partial
h_{t-1}}\right\rvert\right\rvert_{2}\lesssim 1.$ (22)

(a) $T=100$

(b) $T=200$

Figure 1: Parenthesis Task Results

##  4 Experiments

We have performed a variety of experiments to demonstrate the robustness and
efficiency of our proposed NC-GRU method. To this end, we apply NC-GRU to 4
commonly used synthetic tasks: Parenthesis, Denoise, Adding, and Copying
Tasks. In addition, we have considered non-synthetic experiments, Language
Modeling with the character and word level tasks for the Penn Treebank (PTB)
dataset (Marcus et al., 1993a, ).

All the experiments were trained with equal numbers of trainable parameters
also known as parameter-matching architecture, see the implementation details
section of each experiment for more details. Codes for these experiments are
available online at github.com/Edison1994/NC-GRU.

###  4.1 Parenthesis Task

This experiment derives from the descriptions in (Jing et al.,, 2019) and
(Foerster et al.,, 2016). This task tests the ability of the network to
remember the number of unmatched parentheses contained in our input data. The
input data consists of 10 pairs of different types of parentheses combined
with some noise data in between and it is given as a one-hot encoding vector
of length $T$. As stated in (Jing et al.,, 2019), there are not more than $10$
types of parentheses. The output data is given as a one hot-encoding vector
too, counting the number of the unpaired parentheses in the corresponding
input data. The goal of our model is to forget the noise data and absorb
information from the long-term dependencies related to the parentheses. This
synthetic data requires the model to develop a memory and to be able to select
the most relevant information.

For this experiment, we compared 5 different models/architectures: LSTM
(Hochreiter and Schmidhuber, 1997a, ), GRU (Cho et al.,, 2014), scoRNN
(Helfrich et al.,, 2018), GORU (Jing et al.,, 2019), and NC-GRU (ours).

_Implementation Details:_ All of the models consisted of a single layer net
with the following hidden dimensions for each model: LSTM (Hochreiter and
Schmidhuber, 1997a, ) \- 42, GRU (Cho et al.,, 2014) \- 50, scoRNN (Helfrich
et al.,, 2018) \- 110, GORU (Jing et al.,, 2019) \- 64, and NC-GRU - 56.
Furthermore, we trained all of the models for 200 epochs with a batch size of
16, and the number of negative ones for the $D$ matrix in scoRNN (Helfrich et
al.,, 2018) and NC-GRU models were 20 and 40, respectively. All models were
trained using Adam optimizer (Kingma and Ba,, 2014) with a learning rate of
$10^{-3}$ including $A$ associated weights in scoRNN (Helfrich et al.,, 2018)
and NC-GRU models. We conducted experiments using input length of $T=100$, see
Figure 1a, and $T=200$, see Figure 1b.

We present two versions of NC-GRU model, the first one only has orthogonality
in $U_{c}$ weight (NC-GRU($U_{c}$)) however the second model utilizes
orthogonality in both weights $U_{r}$ and $U_{c}$ (NC-GRU($U_{r},U_{c}$)).
Both models were trained using Neumann series method with a reset at every 50
iterations.

_Results:_ We observed and recorded the behavior of the five models on the
Parenthesis task when the input length is set to $100$ and $200$. Our results
showed that both versions of NC-GRU models outperform GRU (Cho et al.,, 2014),
LSTM (Hochreiter and Schmidhuber, 1997a, ), scoRNN (Helfrich et al.,, 2018),
and GORU (Jing et al.,, 2019) models with a significant gap, see Figure 1. On
this task, NC-GRU($U_{r},U_{c}$) model shows a better performance than NC-
GRU($U_{c}$), however both of them perform better than the rest of the
compared models. The minimum value of the loss function attained during the
training are presented in Table 1.

Table 1: Parenthesis Task Results: Minimum attained loss values (the smaller the better). All results are based on our tests. NC-GRU$(U_{c})$ denotes the NC-GRU model (19) where Neumann-Cayley method was only applied to the weight $U_{c}$, similarly NC-GRU$(U_{r},U_{c})$ represents the NC-GRU model (19) where both $U_{r}$ and $U_{c}$ weights were updated using the Neumann-Cayley method. | Loss $\times 10^{-3}$ $\downarrow$
---|---
Sequence Length |  $T=100$ |  $T=200$
LSTM |  $0.62$ |  $6.93$
GRU |  $29.94$ |  $1.526$
scoRNN |  $24.85$ |  $18.957$
GORU |  $1.199$ |  $1.66$
NC-GRU$(U_{c})$ (ours) |  $0.286$ |  $0.227$
NC-GRU$(U_{r},U_{c})$ (ours) |  $0.149$ |  $0.167$

###  4.2 Denoise Task

The Denoise Task (Jing et al.,, 2019) is another synthetic problem that
requires filtering out the noise out of a noisy sequence. This problem
requires the forgetting ability of the network as well as learning long-term
dependencies coming from the data (Jing et al.,, 2019). The input sequence of
length $T$ contains 10 randomly located data points and the other $T-10$
points are considered noise data. These 10 points are selected from a
dictionary $\\{a_{i}\\}_{i=0}^{n+1}$, where the first $n$ elements are data
points, and the other two are the $``noise"$ and the $``marker"$ respectively.
The output data consists of the list of the data points from the input, and it
should be outputted as soon as it receives the $``marker"$. The model task to
filter out the noise part and output the random 10 data points chosen from the
input.

Similar to the Parenthesis Task, we compared our NC-GRU models to LSTM
(Hochreiter and Schmidhuber, 1997a, ), GRU (Cho et al.,, 2014), scoRNN
(Helfrich et al.,, 2018), and GORU (Jing et al.,, 2019) models.

_Implementation Details:_ We implemented one NC-GRU cell with a hidden size of
118 and the number of negative ones in the $D$ matrix to 50. The hidden size
for the LSTM (Hochreiter and Schmidhuber, 1997a, ) was 90, GRU (Cho et al.,,
2014) – 100, scoRNN (Helfrich et al.,, 2018) – 200, and GORU (Jing et al.,,
2019) – 128. We implemented Adam optimizer (Kingma and Ba,, 2014) with
learning rate of $10^{-3}$ to train all the aforementioned models including
scoRNN (Helfrich et al.,, 2018) and NC-GRU $A$ weights. We trained all the
models for 10,000 iterations with a batch size of 128. Similar to the
Parethesis Task 4.1, we implemented the Neumann series method of approximation
the $\left(I+A^{(k)}\right)^{-1}$ when training the NC-GRU model with the
reset option to be applied every 50 iterations.

_Results:_ Based on our experiments, NC-GRU($U_{c}$) and NC-GRU($U_{r},U_{c}$)
models significantly outperformed LSTM (Hochreiter and Schmidhuber, 1997a, ),
GRU (Cho et al.,, 2014), scoRNN (Helfrich et al.,, 2018), and GORU (Jing et
al.,, 2019) models on the denoise data, see Figure 2. Similar to what we
observed from the parenthesis task, usage of two orthogonal weights leads to
better results. In addition, Table 2 provides a comparison of attained minimum
loss for each of the models.

Table 2: Denoise Task Results: Minimum attained loss values (the smaller the better). All results are based on our tests. | Loss $\times 10^{-2}$ $\downarrow$
---|---
Sequence Length |  $T=200$ |  $T=400$
LSTM |  $10.08$ |  $5.427$
GRU |  $9.799$ |  $5.2338$
scoRNN |  $8.258$ |  $3.812$
GORU |  $4.453$ |  $2.29$
NC-GRU$(U_{c})$ (ours) |  $1.639$ |  $0.878$
NC-GRU$(U_{c},U_{r})$ (ours) |  $1.633$ |  $0.774$

(a) $T=200$ (b) $T=400$

Figure 2: Denoise Task Results

###  4.3 Adding Problem

The adding problem is the third synthetic task that we considered and it was
proposed in (Hochreiter and Schmidhuber, 1997b, ) for recurrent networks. Our
implementation of this problem is a variation of the original problem. The
input in the network is a 2-dimensional sequence of length $T$. In the first
dimension, we have a sequence of all zeros except for two randomly placed
ones, one in the first half of the sequence and one in the second one. In the
second dimension, we have a sequence of randomly selected numbers that are
chosen uniformly from the interval $[0,1)$. The goal of the adding task is to
take the second dimension numbers from the positions that correspond to the
ones in the first dimension and then output their sum. The size of the
training and testing/evaluation sets are $100,000$ and $10,000$, respectively.

We compared the performance of our NC-GRU models to LSTM (Hochreiter and
Schmidhuber, 1997a, ), GRU (Cho et al.,, 2014), scoRNN (Helfrich et al.,,
2018), and GORU (Jing et al.,, 2019) models.

_Implementation Details:_ In this task, we worked with a single layer cell
with hidden dimension set to be 68 for LSTM (Hochreiter and Schmidhuber,
1997a, ), 70 for GRU (Cho et al.,, 2014), 190 for scoRNN (Helfrich et al.,,
2018), 128 for GORU (Jing et al.,, 2019), and 80 for NC-GRU. The number of
negative ones for the $D$ matrix inside scoRNN and NC-GRU was set to 95 and
43, respectively. We used the Adam optimizer (Kingma and Ba,, 2014) with a
learning rate of $10^{-3}$ in all of the experiments, including the training
of the $A$ weight in scoRNN (Helfrich et al.,, 2018) and NC-GRU models.

All of the models were trained for 10 epochs with a batch size of 50, and
evaluation every 100 iterations. Neumann series were applied during the
training of NC-GRU models, with reset option implemented every 50 iterations.

_Results:_ Figures 3a and 3b present the performances of all the interested
methods on the tests for sequences of length 200 and 400, respectively. In
these experiment, our NC-GRU($U_{c}$) model where the Neumann-Cayley method
applied only to the $U_{c}$ weight showed the best performance out of all the
compared models including our NC-GRU($U_{r},U_{c}$) model which produced
comparable or marginally better results than others. Minimum attained
validation loss values presented in Table 3 for both experiments.

Table 3: Adding Task Results: Minimum attained loss values (the smaller the better). All results are based on our tests. | Loss $\times 10^{-5}$ $\downarrow$
---|---
Sequence Length |  $T=200$ |  $T=400$
LSTM |  $2.2829$ |  $269.7$
GRU |  $1.776$ |  $1.2976$
scoRNN |  $51.49$ |  $16,295.05$
GORU |  $2.758$ |  $2.3067$
NC-GRU$(U_{c})$ (ours) |  $1.022$ |  $0.81$
NC-GRU$(U_{c},U_{r})$ (ours) |  $5.0$ |  $4.93$

(a) $T=200$

(b) $T=400$

Figure 3: Adding Problem Results

###  4.4 Copying Problem

The copying problem was proposed in (Hochreiter and Schmidhuber, 1997b, ) as a
synthetic task for testing Recurrent Neural Networks. The setup of this
problem consists of a string of 10 random digits that are sampled uniformly
from the integers 1 through 8 and then fed into the recurrent model. These 10
digits are followed by a string of $T$ zeros and a digit 9 which marks the
start of a string of 9 zeros. Therefore, the total length of the fed string is
$T+20$. The objective of the task is to output the initial sequence of 10
random digits beginning at the marker location, copying the first ten elements
of the sequence in order. For the evaluation of the model, the cross-entropy
loss is used with an expected cross-entropy baseline of
$\frac{10\log{(8)}}{T+20}$ which represents the random selection of digits 1-8
after the 9.

Similar to the other 3 synthetic experiments, we have compared the performance
of our NC-GRU models to LSTM (Hochreiter and Schmidhuber, 1997a, ), GRU (Cho
et al.,, 2014), scoRNN (Helfrich et al.,, 2018), and GORU (Jing et al.,, 2019)
models.

_Implementation Details:_ We carried out our experiments on a single layer
models with hidden sizes for each model to match the number of trainable
parameters were 68 for LSTM (Hochreiter and Schmidhuber, 1997a, ), 78 for GRU
(Cho et al.,, 2014), 190 for scoRNN (Helfrich et al.,, 2018), 100 for GORU
(Jing et al.,, 2019), and 96 for NC-GRU. All the models were trained using
Adam (Kingma and Ba,, 2014) optimizer with the learning rate of $10^{-3}$
except $A$ matrix in scoRNN (Helfrich et al.,, 2018) and NC-GRU models where
Adam (Kingma and Ba,, 2014) optimizer learning rate was set to $10^{-4}$. The
number of training iterations and a batch size were set to 10,000 and 50,
respectively.

The number of negative ones for scaling matrix $D$ was set to 95 and 80 for
scoRNN (Helfrich et al.,, 2018) and NC-GRU models respectively. All of the
models were evaluated every 50 iterations. While training the model using the
NC-GRU models, we applied the Neumann series in a similar way as in the other
experiments with the reset option at every 20 iterations.

_Results:_ Figures 4a and 4b present the performance of aforementioned models
for the Copying Problem with string sizes of $T=1,000$ and $T=2,000$,
respectively. Moreover Table 4 shows minimum validation loss values. Both NC-
GRU models perform on the same level while outperforming LSTM (Hochreiter and
Schmidhuber, 1997a, ), GRU (Cho et al.,, 2014), and GORU (Jing et al.,, 2019)
models by a noticeable margin. However, for this task scoRNN (Helfrich et
al.,, 2018) performed better than all of the other models.

Table 4: Copying Task Results: Minimum attained loss values (the smaller the better). All results are based on our tests. | Loss $\times 10^{-2}$ $\downarrow$
---|---
Sequence Length |  $T=1,000$ |  $T=2,000$
LSTM |  $2.040$ |  $1.03$
GRU |  $2.039$ |  $0.89$
scoRNN |  $0.002$ |  $0.03$
GORU |  $1.488$ |  $0.718$
NC-GRU$(U_{c})$ (ours) |  $0.904$ |  $0.44$
NC-GRU$(U_{c},U_{r})$ (ours) |  $0.882$ |  $0.496$

(a) $T=1,000$

(b) $T=2,000$

Figure 4: Copying Problem Results

Language modeling is one of many natural language processing tasks. It is the
development of probabilistic models that are capable of predicting the next
character or word in a sequence using information that has preceded it. We
have conducted two experiments on the language modeling with the character and
word level tasks for the Penn Treebank (PTB) dataset (Marcus et al., 1993a, )
both of which were based on experiments and codes from the AWD-LSTM model
(Merity et al.,, 2018).

###  4.5 Character Level Penn Treebank

For this task, models were tested on their suitability for language modeling
tasks using the character level Penn Treebank dataset (Marcus et al., 1993b,
). This dataset is a collection of English-language Wall Street Journal
articles. The dataset consists of a vocabulary of 10,000 words with other
words replaced as $<unk>$, resulting in approximately 6 million characters
that are divided into 5.1 million, 400 thousand, and 450 thousand character
sets for training, validation, and testing, respectively with a character
alphabet size of 50. The goal of the character-level Language Modeling task is
to predict the next character given the preceding sequence of characters.

We compared the performance of our NC-GRU model to LSTM (Hochreiter and
Schmidhuber, 1997a, ), GRU (Cho et al.,, 2014), and scoRNN (Helfrich et al.,,
2018) models.

_Implementation Details:_ For this task, we have considered a three layer
models, where hidden dimensions for NC-GRU($U_{c}$) model were set to (430,
1000, 430). All the dropout coefficient were set to $0.15$. The learning rate
of $5\times 10^{-4}$ together with the Adam (Kingma and Ba,, 2014) optimizer
were used to train whole model including the matrix $A$. The number of
negative ones for matrix $D$ for every layer was set to $50$.

LSTM (Hochreiter and Schmidhuber, 1997a, ), GRU (Cho et al.,, 2014), and
scoRNN (Helfrich et al.,, 2018) model were trained in the similar manner with
hidden dimensions being (350, 880, 350), (415, 950, 415), and (500, 2000,
500), respectively. Batch size for all the experiments was set to be $32$ and
the Back Propagation Through Time (bptt) window of 100 was used for all the
models.

_Results:_ We evaluate all models using the bits-per-character (bpc) metric.
As shown in Table 5, NC-GRU showed a better performance than scoRNN, GRU, and
even LSTM models.

Table 5: Character Level PTB Results: Evaluated bits-per-character (bpc) for every model (the smaller the better). All results are based on our tests. |  bpc $\downarrow$
---|---
LSTM |  $1.4$
GRU |  $1.455$
scoRNN |  $1.6$
NC-GRU($U_{c}$) (ours) |  $1.385$

###  4.6 Word Level Penn Treebank

Finally, we also tested our proposed Neumann-Cayley method on the word level
Penn Treebank dataset (Marcus et al., 1993a, ). The dataset takes the same
underlying corpus as the character level task, but with tokens representing
words instead of characters. This results in a smaller dataset with a larger
vocabulary size, with 888 thousand, 70 thousand, and 79 thousand words as
training, validation, and testing sets, and a vocabulary of 10 thousand words.

_Implementation Details:_ We trained NC-GRU with three layers with dimensions
(400, 1150, 400). We have a learning rate set to $5\times 10^{-4}$ for both
the $A$ matrix and the rest of the model, optimized using Adam (Kingma and
Ba,, 2014). The dropout after the NC-GRU cell has a coefficient of 0.4, the
embedding layer dropout has a coefficient of 0.4, and the output dropout has a
coefficient of 0.25. The number of negative ones for matrix $D$ for each layer
was 50.

_Results:_ Results were evaluated using the perplexity (PPL) metric and are
shown in table 6. We show improved performance over both baseline GRU and LSTM
models.

Table 6: Word Level PTB Results: Evaluated perplexity (PPL) for every model (the smaller the better). ∗ \- result obtained from our experiments; † \- result quoted from (Bai et al.,, 2018) |  PPL $\downarrow$
---|---
LSTM |  $78.93^{\dagger}$
GRU |  $92.48^{\dagger}$ ($80.73^{*}$)
scoRNN |  $123.90^{*}$
NC-GRU($U_{c}$) (ours) |  $77.00$

##  5 Ablation Studies

In this section, we consider several ablation studies that help us justify the
use of Neumann series, orthogonal matrices, as well as scaled-Cayley
transforms.

###  5.1 Neumann Series Method vs Inverse

In this experiment we study the sharpness of approximating
$\left(I+A^{(k-1)}\right)^{-1}$ with Neumann series in the NC-GRU($U_{c}$)
model on the Parenthesis task, see 4.1 for implementation details and
description of NC-GRU($U_{c}$) model. We consider using Neumann series
approximation of order 1, 2, and 3 and compare them to the method of Least
Squares for taking a matrix inverse which is one of the widely used methods
from Deep Learning libraries such as TensorFlow, PyTorch, and NumPy.

Figure 5: Inverse vs Neumann series: In this graphs, Inverse represents NC-
GRU($U_{c}$) model with $\left(I+A^{(k-1)}\right)^{-1}$ computed using the
Least Squares method and Neumann i represents NC-GRU($U_{c}$) model with
$\left(I+A^{(k-1)}\right)^{-1}$ computed using the Neumann series method of
order $i$.

Our experiments showed that the Neumann series approximation method achieves
better results in comparison to the classical Least Squares method for finding
matrix inverse. Particularly, Figure 5 shows that the Neumann series method of
order 2 performs marginally better than order 1 and order 3 Neumann series
methods and significantly better than the Least Squares method.

Figure 6: Computational time comparisons: The amount of time (in seconds) it
takes to train one epoch of NC-GRU($U_{c}$) model on Character Level PTB
dataset using a single NVIDIA® Tesla® V100 GPU.

Additionally, we have compared the time it takes to train our model using
mentioned methods. Figure 6 shows the time, it takes to train one epoch of NC-
GRU($U_{c}$) model on a Character Level PTB dataset on a single NVIDIA® Tesla®
V100 GPU using Inverse (Least Square), Neumann 1, Neumann 2, and Neumann 3
methods.

The observed behavior appears to be quite general and we have conducted all of
the experiments in section 4 using the second order Neumann series method.

###  5.2 Orthogonality in the Hidden Weights of GRU

For our second ablation study, we have studied the effect of Neumann-Cayley
transformation orthogonal weights applied to the hidden units inside the GRU
cell (1). We have considered three models. The first model only had $U_{c}$
weight replaced with orthogonal weight preserved by the Neumann-Cayley method,
we previously called such a model NC-GRU($U_{c}$). The second model, NC-
GRU($U_{r},U_{c}$) had two weights $U_{c}$ and $U_{r}$ replaced with Neumann-
Cayley transformation orthogonal weights, and finally, the third model, NC-
GRU($U_{r},U_{u},U_{c}$) had all three weights $U_{r}$, $U_{u}$, and $U_{c}$
replaced.

The results are shown in Figure 7. We see that implementing one or two
orthogonal weight models, NC-GRU($U_{c}$) or NC-GRU($U_{r},U_{c}$), would have
the most benefits, while the three orthogonal weight model NC-
GRU($U_{r},U_{u},U_{c}$) does not perform as well. This can also be seen in
experiments in section 4 where both one and two orthogonal weight models are
used.

Figure 7: Orthogonal matrices for different combinations of hidden weights:
NC-GRU one hidden parameter orthogonal vs three hidden parameters orthogonal

###  5.3 Necessary Condition for the Neumann Series Method

As it was mentioned in 3.2, the assumption that allowed us to use the Neumann
series is $\left\lvert\left\lvert\left(I+A^{(k-1)}\right)^{-1}\delta
A^{(k)}\right\rvert\right\rvert<1$ where
$\left\lvert\left\lvert\cdot\right\rvert\right\rvert$ satisfies
$\left\lvert\left\lvert AB\right\rvert\right\rvert\leq\left\lvert\left\lvert
A\right\rvert\right\rvert\left\lvert\left\lvert B\right\rvert\right\rvert$ for
some matrices $A$ and $B$ of appropriate dimensions.

Figure 8: Norm condition for the Neumann approximation: Values of
$\left\lvert\left\lvert\left(I+A^{(k-1)}\right)^{-1}\delta
A^{(k)}\right\rvert\right\rvert_{2}$ using NC-GRU($U_{c}$) model with second-
order Neumann series approximation method.

For this experiment, we chose the Spectral Norm, i.e. $L_{2}$-norm, and
recorded values $\left\lvert\left\lvert\left(I+A^{(k-1)}\right)^{-1}\delta
A^{(k)}\right\rvert\right\rvert_{2}$ during training the NC-GRU($U_{c}$) model
on the Parenthesis task with the second-order Neumann series approximation. We
observed that values of the norm were changing during the training however
they are well below 1 and satisfy the necessary condition as depicted in
Figure 8.

##  6 Conclusion

In this paper, we have presented a thorough analysis of the Gated Recurrent
Unit (GRU) model’s gradients. Based on this analysis, we introduced the
Neumann-Cayley Gated Recurrent Unit model, which we called NC-GRU. Our model
incorporates orthogonal weights in the hidden states of the GRU model which
are trained using a newly proposed method of Neumann-Cayley transformation for
maintaining the desired orthogonality in those weights. We have conducted a
series of experiments that demonstrate the superiority of our proposed method
outperforming GRU, LSTM, scoRNN, and GORU moduls on different synthetic and
real-world tasks. Moreover, we conducted several ablation studies that
empirically confirmed our theoretical results.

##  7 Acknowledgment

We would like to thank the University of Kentucky Center for Computational
Sciences and Information Technology Services Research Computing for their
support and use of the Lipscomb Compute Cluster and associated research
computing resources. This research was supported in part by NSF under grants
DMS-1821144, DMS-2053284 and DMS-2151802, and University of Kentucky Start-up
fund.

## References

  * Arjovsky et al., (2016) Arjovsky, M., Shah, A., and Bengio, Y. (2016).  Unitary evolution recurrent neural networks.  In Proceedings of the 33rd International Conference on International Conference on Machine Learning - Volume 48, ICML’16, page 1120–1128. JMLR.org.
  * Ba et al., (2016) Ba, J. L., Kiros, J. R., and Hinton, G. E. (2016).  Layer normalization.
  * Bai et al., (2018) Bai, S., Kolter, J. Z., and Koltun, V. (2018).  An empirical evaluation of generic convolutional and recurrent networks for sequence modeling.
  * Bengio et al., (1993) Bengio, Y., Frasconi, P., and Simard, P. (1993).  The problem of learning long-term dependencies in recurrent networks.  In Proceedings of 1993 IEEE International Conference on Neural Networks (ICNN ’93), pages 1183–1195, San Francisco, CA. IEEE Press.
  * Cho et al., (2014) Cho, K., van Merrienboer, B., Gulcehre, C., Bahdanau, D., Bougares, F., Schwenk, H., and Bengio, Y. (2014).  Learning phrase representations using rnn encoder-decoder for statistical machine translation.
  * Cooijmans et al., (2017) Cooijmans, T., Ballas, N., Laurent, C., Gülçehre, Ç., and Courville, A. C. (2017).  Recurrent batch normalization.  In 5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings. OpenReview.net.
  * Demmel, (1997) Demmel, J. W. (1997).  Applied Numerical Linear Algebra.  SIAM.
  * Dorado-Rojas et al., (2020) Dorado-Rojas, S. A., Vinzamuri, B., and Vanfretti, L. (2020).  Orthogonal laguerre recurrent neural networks.  34th Conference on Neural Information Processing Systems.
  * Foerster et al., (2016) Foerster, J. N., Gilmer, J., Chorowski, J., Sohl-Dickstein, J., and Sussillo, D. (2016).  Input switched affine networks: An rnn architecture designed for interpretability.
  * Helfrich et al., (2018) Helfrich, K., Willmott, D., and Ye, Q. (2018).  Orthogonal Recurrent Neural Networks with Scaled Cayley Transform.  In Proceedings of ICML 2018, Stockholm, Sweden, PMLR 80.
  * Helfrich and Ye, (2019) Helfrich, K. and Ye, Q. (2019).  Eigenvalue normalized recurrent neural networks for short term memory.
  * Henaff et al., (2017) Henaff, M., Szlam, A., and LeCun, Y. (2017).  Recurrent orthogonal networks and long-memory tasks.  In Proceedings of the 33rd International Conference on Machine Learning (ICML 2017), volume 48.
  * (13) Hochreiter, S. and Schmidhuber, J. (1997a).  Long short-term memory.  Neural Comput., 9(8):1735–1780.
  * (14) Hochreiter, S. and Schmidhuber, J. (1997b).  Long short-term memory.  Neural computation, 9(8):1735–1780.
  * Hopfield, (1982) Hopfield, J. J. (1982).  Neural networks and physical systems with emergent collective computational abilities.  Proceedings of the National Academy of Sciences, 79(8):2554–2558.
  * Hyland and Gunnar, (2017) Hyland, S. L. and Gunnar, R. (2017).  Learning unitary operators with help from u(n).  In Proceedings of the 31st AAAI Conference on Artificial Intelligence (AAI 2017), pages 2050–2058, San Francisco, CA.
  * Ioffe and Szegedy, (2015) Ioffe, S. and Szegedy, C. (2015).  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  arXiv preprint arXiv:1502.03167.
  * Jing et al., (2019) Jing, L., Gulcehre, C., Peurifoy, J., Shen, Y., Tegmark, M., Soljacic, M., and Bengio, Y. (2019).  Gated orthogonal recurrent units: On learning to forget.  Neural Computation, 31(4):765–783.
  * Jing et al., (2017) Jing, L., Shen, Y., Dubcek, T., Peurifoy, J., Skirlo, S., LeCun, Y., Tegmark, M., and Soljačić, M. (2017).  Tunable efficient unitary neural networks (eunn) and their application to rnns.  In Proceedings of the 34th International Conference on Machine Learning - Volume 70, ICML’17, page 1733–1741. JMLR.org.
  * Kahan, (2006) Kahan, W. (2006).  Is there a small skew cayley transform with zero diagonal?  Linear Algebra and its Applications, 417(2):335–341.  Special Issue in honor of Friedrich Ludwig Bauer.
  * Kingma and Ba, (2014) Kingma, D. and Ba, J. (2014).  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980.
  * Lezcano-Casado and Martínez-Rubio, (2019) Lezcano-Casado, M. and Martínez-Rubio, D. (2019).  Cheap orthogonal constraints in neural networks: A simple parametrization of the orthogonal and unitary group.  In Chaudhuri, K. and Salakhutdinov, R., editors, Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pages 3794–3803. PMLR.
  * Maduranga et al., (2018) Maduranga, K. D. G., Helfrich, K. E., and Ye, Q. (2018).  Complex unitary recurrent neural networks using scaled cayley transform.
  * (24) Marcus, M. P., Marcinkiewicz, M. A., and Santorini, B. (1993a).  Building a large annotated corpus of english: The penn treebank.  Comput. Linguist., 19(2):313–330.
  * (25) Marcus, M. P., Santorini, B., and Marcinkiewicz, M. A. (1993b).  Building a large annotated corpus of English: The Penn Treebank.  Computational Linguistics, 19(2):313–330.
  * Merity et al., (2018) Merity, S., Keskar, N. S., and Socher, R. (2018).  Regularizing and optimizing LSTM language models.  In International Conference on Learning Representations.
  * Mhammedi et al., (2017) Mhammedi, Z., Hellicar, A. D., Rahman, A., and Bailey, J. (2017).  Efficient orthogonal parametrisation of recurrent neural networks using householder reflections.  In Proceedings of ICML 2017, volume 70, pages 2401–2409. PMLR.
  * Rumelhart et al., (1986) Rumelhart, D. E., Hinton, G. E., and Williams, R. J. (1986).  Learning Internal Representations by Error Propagation, page 318–362.  MIT Press.
  * Salimans and Kingma, (2016) Salimans, T. and Kingma, D. P. (2016).  Weight normalization: A simple reparameterization to accelerate training of deep neural networks.  In Lee, D., Sugiyama, M., Luxburg, U., Guyon, I., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 29. Curran Associates, Inc.
  * Saxe et al., (2014) Saxe, A. M., McClelland, J. L., and Ganguli, S. (2014).  Exact solutions to nonlinear dynamics of learning in deep linear neural networks.
  * Tagare, (2011) Tagare, H. D. (2011).  Notes on optimization on stiefel manifolds.  Yale.  ”Technical report, Yale University”.
  * Tieleman and Hinton, (2012) Tieleman, T. and Hinton, G. (2012).  Lecture 6.5-rmsprop: Divide the gradient by a running average of its recent magnitude.  COURSERA: Neural networks for machine learning, 4(2).
  * Ulyanov et al., (2017) Ulyanov, D., Vedaldi, A., and Lempitsky, V. (2017).  Instance normalization: The missing ingredient for fast stylization.
  * (34) Vorontsov, E., Trabelsi, C., Kadoury, S., and Pal, C. (2017a).  On orthogonality and learning recurrent networks with long term dependencies.
  * (35) Vorontsov, E., Trabelsi, C., Kadoury, S., and Pal, C. (2017b).  On orthogonality and learning recurrent networks with long term dependencies.  arXiv preprint arXiv:1702.00071.
  * Wisdom et al., (2016) Wisdom, S., Powers, T., Hershey, J., Le Roux, J., and Atlas, L. (2016).  Full-capacity unitary recurrent neural networks.  In Advances in Neural Information Processing Systems 29, pages 4880–4888.
  * Wu and He, (2018) Wu, Y. and He, K. (2018).  Group normalization.  In Proceedings of the European Conference on Computer Vision (ECCV).
  * Xu et al., (2019) Xu, J., Sun, X., Zhang, Z., Zhao, G., and Lin, J. (2019).  Understanding and improving layer normalization.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc.

##  8 Supplementary Materials

In this section, we present the proofs for the theorem and corollaries stated
in section 3.

######  Theorem 1.

Let $h_{t-1}$ and $h_{t}$ be two consecutive hidden states from the GRU model
stated in (1). Then

$\left\lvert\left\lvert\dfrac{\partial h_{t}}{\partial
h_{t-1}}\right\rvert\right\rvert_{2}\leq\alpha+\beta\left\lvert\left\lvert
U_{c}\right\rvert\right\rvert_{2}$ (23)

where

$\displaystyle\alpha$
$\displaystyle=\delta_{u}\left(\max_{i}\left\\{[h_{t-1}]_{i}\right\\}+\max_{i}\left\\{[c_{t}]_{i}\right\\}\right)\left\lvert\left\lvert
U_{u}\right\rvert\right\rvert_{2}+\max_{i}\left\\{(1-[u_{t}]_{i})\right\\}$
(24)

and

$\displaystyle\beta$
$\displaystyle=\max_{i}\left\\{[u_{t}]_{i}\right\\}\left(\delta_{r}\left\lvert\left\lvert
U_{r}\right\rvert\right\rvert_{2}\max_{i}\left\\{[h_{t-1}]_{i}\right\\}+\max_{i}\left\\{[r_{t}]_{i}\right\\}\right),$
(25)

with constants $\delta_{u}$ and $\delta_{r}$ be as follows:

$\delta_{u}=\max_{i}\left\\{\left[u_{t}\right]_{i}\left(1-\left[u_{t}\right]_{i}\right)\right\\}$
(26)

and

$\delta_{r}=\max_{i}\left\\{\left[r_{t}\right]_{i}\left(1-\left[r_{t}\right]_{i}\right)\right\\}.$
(27)

###### Proof.

Since $h_{t}$ depends on $u_{t}$ and $c_{t}$ (which is also depends on
$r_{t}$), we will start by finding $\dfrac{\partial r_{t}}{\partial h_{t-1}}$,
and $\dfrac{\partial u_{t}}{\partial h_{t-1}}$, and $\dfrac{\partial
c_{t}}{\partial h_{t-1}}$ as well as corresponding bounds of these gradients.
Recall that

$r_{t}=\sigma\left(W_{r}x_{t}+U_{r}h_{t-1}+b_{r}\right)$ (28)

then

$\dfrac{\partial r_{t}}{\partial
h_{t-1}}=diag\left(\sigma^{\prime}\left(W_{r}x_{t}+U_{r}h_{t-1}+b_{r}\right)\right)U_{r}=diag\left(r_{t}\odot(1-r_{t})\right)U_{r}$
(29)

since $\sigma^{\prime}(x)=\sigma(x)\left(1-\sigma(x)\right)$. Moreover, we can
bound $\dfrac{\partial r_{t}}{\partial h_{t-1}}$ in the following way

$\left\lvert\left\lvert\dfrac{\partial r_{t}}{\partial
h_{t-1}}\right\rvert\right\rvert_{2}\leq\delta_{r}\left\lvert\left\lvert
U_{r}\right\rvert\right\rvert_{2}$ (30)

with

$\delta_{r}:=\max_{i}\left\\{\left[r_{t}\right]_{i}\left(1-\left[r_{t}\right]_{i}\right)\right\\}.$
(31)

The constant $\delta_{r}$ is defined to be the largest entry of the vector
$r_{t}\odot(1-r_{t})$ and it is bounded by $\dfrac{1}{4}$. Similarly,

$\dfrac{\partial u_{t}}{\partial
h_{t-1}}=diag\left(\sigma^{\prime}\left(W_{u}x_{t}+U_{u}h_{t-1}+b_{u}\right)\right)U_{u}=diag\left(u_{t}\odot(1-u_{t})\right)U_{u},$
(32)

and

$\left\lvert\left\lvert\dfrac{\partial u_{t}}{\partial
h_{t-1}}\right\rvert\right\rvert_{2}\leq\delta_{u}\left\lvert\left\lvert
U_{u}\right\rvert\right\rvert_{2}$ (33)

with

$\delta_{u}:=\max_{i}\left\\{\left[u_{t}\right]_{i}\left(1-\left[u_{t}\right]_{i}\right)\right\\}.$
(34)

The definition of the vector $c_{t}$ is

$c_{t}=\Phi\left(W_{c}x_{t}+U_{c}\left(r_{t}\odot h_{t-1}\right)+b_{c}\right)$
(35)

and

$\dfrac{\partial c_{t}}{\partial
h_{t-1}}=diag\left(\Phi^{\prime}\left(W_{c}x_{t}+U_{c}(r_{t}\odot
h_{t-1})+b_{c}\right)\right)U_{c}\left(diag\left(h_{t-1}\right)\dfrac{\partial
r_{t}}{\partial h_{t-1}}+diag\left(r_{t}\right)\right),$ (36)

with $\Phi^{\prime}$ applied entrywise, and

$\displaystyle\left\lvert\left\lvert\dfrac{\partial c_{t}}{\partial
h_{t-1}}\right\rvert\right\rvert_{2}$
$\displaystyle\leq\max_{i}\left\\{\left[\Phi^{\prime}\left(W_{c}x_{t}+U_{c}\left(r_{t}\odot
h_{t-1}\right)+b_{c}\right)\right]_{i}\right\\}\left(\left\lvert\left\lvert\dfrac{\partial
r_{t}}{\partial
h_{t-1}}\right\rvert\right\rvert_{2}\max_{i}\left\\{[h_{t-1}]_{i}\right\\}+\max_{i}\left\\{[r_{t}]_{i}\right\\}\right)\left\lvert\left\lvert
U_{c}\right\rvert\right\rvert_{2}$ (37)
$\displaystyle\leq\left(\delta_{r}\left\lvert\left\lvert
U_{r}\right\rvert\right\rvert_{2}\max_{i}\left\\{[h_{t-1}]_{i}\right\\}+\max_{i}\left\\{[r_{t}]_{i}\right\\}\right)\left\lvert\left\lvert
U_{c}\right\rvert\right\rvert_{2}$ (38)

since $\Phi^{\prime}$ is bounded by 1.

Finally,

$h_{t}=\left(1-u_{t}\right)\odot h_{t-1}+u_{t}\odot c_{t}$ (39)

with

$\dfrac{\partial h_{t}}{\partial
h_{t-1}}=-diag\left(h_{t-1}\right)\dfrac{\partial u_{t}}{\partial
h_{t-1}}+diag\left(1-u_{t}\right)+diag\left(c_{t}\right)\dfrac{\partial
u_{t}}{\partial h_{t-1}}+diag\left(u_{t}\right)\dfrac{\partial c_{t}}{\partial
h_{t-1}}$ (40)

and

$\left\lvert\left\lvert\dfrac{\partial h_{t}}{\partial
h_{t-1}}\right\rvert\right\rvert_{2}\leq\max_{i}\left\\{[h_{t-1}]_{i}\right\\}\left\lvert\left\lvert\dfrac{\partial
u_{t}}{\partial
h_{t-1}}\right\rvert\right\rvert_{2}+\max_{i}\left\\{(1-[u_{t}]_{i})\right\\}+\max_{i}\left\\{[c_{t}]_{i}\right\\}\left\lvert\left\lvert\dfrac{\partial
u_{t}}{\partial
h_{t-1}}\right\rvert\right\rvert_{2}+\max_{i}\left\\{[u_{t}]_{i}\right\\}\left\lvert\left\lvert\dfrac{\partial
c_{t}}{\partial h_{t-1}}\right\rvert\right\rvert_{2}.$ (41)

Furthermore, using (30), (33), and (38), we get

$\displaystyle\left\lvert\left\lvert\dfrac{\partial h_{t}}{\partial
h_{t-1}}\right\rvert\right\rvert_{2}$
$\displaystyle\leq\delta_{u}\left(\max_{i}\left\\{[h_{t-1}]_{i}\right\\}+\max_{i}\left\\{[c_{t}]_{i}\right\\}\right)\left\lvert\left\lvert
U_{u}\right\rvert\right\rvert_{2}+\max_{i}\left\\{(1-[u_{t}]_{i})\right\\}$
$\displaystyle\quad\quad+\max_{i}\left\\{[u_{t}]_{i}\right\\}\left(\delta_{r}\left\lvert\left\lvert
U_{r}\right\rvert\right\rvert_{2}\max_{i}\left\\{[h_{t-1}]_{i}\right\\}+\max_{i}\left\\{[r_{t}]_{i}\right\\}\right)\left\lvert\left\lvert
U_{c}\right\rvert\right\rvert_{2}$ (42)
$\displaystyle=:\alpha+\beta\left\lvert\left\lvert
U_{c}\right\rvert\right\rvert_{2}$ (43)

where

$\alpha=\delta_{u}\left(\max_{i}\left\\{[h_{t-1}]_{i}\right\\}+\max_{i}\left\\{[c_{t}]_{i}\right\\}\right)\left\lvert\left\lvert
U_{u}\right\rvert\right\rvert_{2}+\max_{i}\left\\{(1-[u_{t}]_{i})\right\\}$
(44)

and

$\beta=\max_{i}\left\\{[u_{t}]_{i}\right\\}\left(\delta_{r}\left\lvert\left\lvert
U_{r}\right\rvert\right\rvert_{2}\max_{i}\left\\{[h_{t-1}]_{i}\right\\}+\max_{i}\left\\{[r_{t}]_{i}\right\\}\right).$
(45)

∎

######  Corollary 1.

For the hyperbolic tangent activation function in (1) (i.e. $\Phi=\,$`tanh`),
we have $\delta_{u},\delta_{r}\leq\frac{1}{4}$, $[h_{t}]_{i}\leq 1$ for any
$i$ and $t$ as well as

$\alpha\leq\dfrac{1}{2}\left\lvert\left\lvert
U_{u}\right\rvert\right\rvert_{2}+1\quad\text{and}\quad\beta\leq\dfrac{1}{4}\left\lvert\left\lvert
U_{r}\right\rvert\right\rvert_{2}+1.$ (46)

###### Proof.

The function $\sigma^{\prime}(x)=\sigma(x)\left(1-\sigma(x)\right)$ is bounded
above by $\frac{1}{4}$, thus both $\delta_{u},\delta_{r}\leq\frac{1}{4}$.

Now, to show that $[h_{t}]_{i}\leq 1$ for any $i$ and $t$, we first need to
note that $h_{0}$ is initialized to zero (i.e. $[h_{0}]_{i}=0$ for all $i$)
and $0\leq[u_{t}]_{i},[c_{t}]_{i}\leq 1$ for all $i$ and $t$ from the
definition of GRU cell in (1). Then for any fixed $i$

$\displaystyle[h_{0}]_{i}$ $\displaystyle=0$ (47) $\displaystyle[h_{1}]_{i}$
$\displaystyle=[1-u_{1}]_{i}\cdot[h_{0}]_{i}+[u_{1}]_{i}\cdot[c_{1}]_{i}$ (48)
$\displaystyle=[u_{1}]_{i}\cdot[c_{1}]_{i}\leq 1$ (49)

Furthermore, if we assume that $[h_{\tau}]_{i}\leq 1$ for some $\tau\geq 1$,
then

$\displaystyle[h_{\tau+1}]_{i}$
$\displaystyle=[1-u_{\tau+1}]_{i}\cdot[h_{\tau}]_{i}+[u_{\tau+1}]_{i}\cdot[c_{\tau+1}]_{i}$
(50) $\displaystyle\leq[1-u_{\tau+1}]_{i}\cdot 1+[u_{\tau+1}]_{i}\cdot 1$ (51)
$\displaystyle=[1-u_{\tau+1}]_{i}+[u_{\tau+1}]_{i}=1.$ (52)

Thus, by induction we can conclude that $[h_{t}]_{i}\leq 1$ for any $i$ and
$t$. Finally, using these obtained bounds, we can bound constants $\alpha$ and
$\beta$ as follows

$\displaystyle\alpha$
$\displaystyle=\delta_{u}\left(\max_{i}\left\\{[h_{t-1}]_{i}\right\\}+\max_{i}\left\\{[c_{t}]_{i}\right\\}\right)\left\lvert\left\lvert
U_{u}\right\rvert\right\rvert_{2}+\max_{i}\left\\{(1-[u_{t}]_{i})\right\\}$
(53)
$\displaystyle\leq\dfrac{1}{4}\cdot\left(1+1\right)\cdot\left\lvert\left\lvert
U_{u}\right\rvert\right\rvert_{2}+1$ (54)
$\displaystyle=\dfrac{1}{2}\left\lvert\left\lvert
U_{u}\right\rvert\right\rvert_{2}+1$ (55)

and

$\displaystyle\beta$
$\displaystyle=\max_{i}\left\\{[u_{t}]_{i}\right\\}\left(\delta_{r}\left\lvert\left\lvert
U_{r}\right\rvert\right\rvert_{2}\max_{i}\left\\{[h_{t-1}]_{i}\right\\}+\max_{i}\left\\{[r_{t}]_{i}\right\\}\right)$
(56) $\displaystyle\leq 1\cdot\left(\dfrac{1}{4}\cdot\left\lvert\left\lvert
U_{r}\right\rvert\right\rvert_{2}\cdot 1+1\right)$ (57)
$\displaystyle=\dfrac{1}{4}\left\lvert\left\lvert
U_{r}\right\rvert\right\rvert_{2}+1.$ (58)

Note, that all of these bounds are independent of $i$ and $t$.

∎

The following two corollaries use a notation $x\approx y$ and $x\lesssim y$ to
denote that $x$ is approximately equal to $y$ and $x$ is bounded by a quantity
that is approximately equal to $y$, respectively.

######  Corollary 2.

When the elements of GRU gates $u_{t}$ and $r_{t}$ are nearly either 0 or 1,
then constants $\alpha$ and $\beta$ from Theorem 1 satisfy the following
inequality:

$\alpha+\beta\lesssim 2.$ (59)

Moreover if $u_{t}$ and $r_{t}$ are nearly either the zero vector or the
vector of all ones, then

$\alpha+\beta\lesssim 1.$ (60)

###### Proof.

Recall the definition of $\alpha$ and $\beta$

$\alpha=\delta_{u}\left(\max_{i}\left\\{[h_{t-1}]_{i}\right\\}+\max_{i}\left\\{[c_{t}]_{i}\right\\}\right)\left\lvert\left\lvert
U_{u}\right\rvert\right\rvert_{2}+\max_{i}\left\\{(1-[u_{t}]_{i})\right\\}$
(61)

and

$\beta=\max_{i}\left\\{[u_{t}]_{i}\right\\}\left(\delta_{r}\left\lvert\left\lvert
U_{r}\right\rvert\right\rvert_{2}\max_{i}\left\\{[h_{t-1}]_{i}\right\\}+\max_{i}\left\\{[r_{t}]_{i}\right\\}\right).$
(62)

If we assume that $[u_{t}]_{i}\approx 0$ and $[u_{t}]_{j}\approx 1$ for some
$i\neq j$, then $\delta_{u}\approx 0$ and $\alpha\lesssim 1$. Additionally, if
we assume that the elements of $r_{t}$ are nearly either 0 or 1, then
$\delta_{r}\approx 0$ and $\beta\lesssim 1$. Putting these two inequalities
together yields

$\alpha+\beta\lesssim 2.$ (63)

Now, if we assume that $u_{t}$ is nearly a zero vector (i.e.
$\max_{i}\left\\{[u_{t}]_{i}\right\\}\approx 0$), then constant
$\delta_{u}\approx 0$, $\alpha\approx 1$, and $\beta\approx 0$. On the other
hand, if we assume that $u_{t}$ is nearly the vector of all ones (i.e.
$\min_{i}\left\\{[u_{t}]_{i}\right\\}\approx 1$) then $\delta_{u}\approx 0$
and $\alpha\approx 0$. Moreover, if we also assume that $r_{t}$ is nearly a
zero vector (i.e. $\max_{i}\left\\{[r_{t}]_{i}\right\\}\approx 0$) then
$\delta_{r}\approx 0$ and $\beta\approx 0$. However, if we assume that $r_{t}$
approaches a vector of all ones (i.e.
$\min_{i}\left\\{[r_{t}]_{i}\right\\}\approx 1$), then $\delta_{r}\approx 0$
but $\beta\approx 1$ for this case. Putting all of these cases together, we
obtain

$\alpha+\beta\lesssim 1.$ (64)

∎

######  Corollary 3.

Let $h_{t-1}$ and $h_{t}$ be two consecutive hidden states from the NC-GRU
model defined in (19). Then $\left\lvert\left\lvert
U_{r}\right\rvert\right\rvert_{2}=\left\lvert\left\lvert
U_{c}\right\rvert\right\rvert_{2}=1$ and if element of the gates $u_{t}$ and
$r_{t}$ are nearly either 0 or 1, then the following inequality is satisfied:

$\left\lvert\left\lvert\dfrac{\partial h_{t}}{\partial
h_{t-1}}\right\rvert\right\rvert_{2}\lesssim 2.$ (65)

Furthermore, if $u_{t}$ and $r_{t}$ are nearly either the zero vector or the
vector of all ones,

$\left\lvert\left\lvert\dfrac{\partial h_{t}}{\partial
h_{t-1}}\right\rvert\right\rvert_{2}\lesssim 1.$ (66)

###### Proof.

By the definition of NC-GRU model, weights $U_{r}$ and $U_{c}$ are orthogonal
which implies $\left\lvert\left\lvert
U_{r}\right\rvert\right\rvert_{2}=\left\lvert\left\lvert
U_{c}\right\rvert\right\rvert_{2}=1$, and by Theorem 1 and Corollary 2, we
conclude

$\left\lvert\left\lvert\dfrac{\partial h_{t}}{\partial
h_{t-1}}\right\rvert\right\rvert_{2}\lesssim 2$ (67)

when element of the gates $u_{t}$ and $r_{t}$ approach either 0 or 1; and

$\left\lvert\left\lvert\dfrac{\partial h_{t}}{\partial
h_{t-1}}\right\rvert\right\rvert_{2}\lesssim 1$ (68)

if $u_{t}$ and $r_{t}$ approach either zero vector or vector of all ones.

∎


Figure 3 (d). Thus, we demonstrated that MuSe-GNN can learn gene embeddings
from multiple species and identify genes with similar functions across various
species.

###  Appendix K Lung cancer analysis

In this section, we applied MuSe-GNN to analyze differentially co-expressed
genes in two types of datasets, comprising lung cancer samples and healthy
samples. We converted the scRNA-seq results into distinct graphs according to
our graph construction strategy and trained MuSe-GNN to generate gene
embeddings. We then employed the Leiden algorithm to identify common
functional gene groups. After collecting essential genes with cancer-related
special functions, we used GOEA and IPA to uncover more biological
information.

As shown in Figure 19 (a)-(c), we observed that for lung tissue, the co-
expression relationships of genes in lung tumor cells and normal cells differ
significantly. Figure 19 (d) displays the top pathways identified by GOEA for
genes recognized as specifically co-expressed in cancer cells, with most
pathways related to tissue development and cell migration. These pathways
align with the fact that cancer cells often undergo multiple divisions and
frequently relocate. Furthermore, using IPA, we analyzed causal networks and
disease functions. Figure 19 (e) illustrates an example of a causal network
discovered in a specific gene group, primarily regulated by ethyl
protocatechuate, indicating a chemical-dominated regulatory network. Thus,
with the genes identified by MuSe-GNN, we can uncover various specific
regulatory networks present in cancer and investigate strategies to inhibit
cancer cell growth and metastasis. Lastly, Figure K (f) reveals that functions
with low FDR rates include tissue development and cancer disease.

Our analysis demonstrates that MuSe-GNN can be utilized to discover functional
gene clusters exhibiting specificity in cancer, assisting researchers in
gaining a deeper understanding of gene expression and regulatory relationships
in cancer.

Figure 19: Gene representation learning results for samples with cancer and
healthy samples. (a) represents the UMAPs of gene embeddings colored by
functional groups. (b) represents the UMAPs of gene embeddings colored by
datasets. (c) represents the gene embeddings colored by the conditions, and
the red circle reflects the differential co-expression genes. (d) shows the
top6 pathways related to the genes in the special cluster discovered by GOEA.
(e) represents the causal network existing in the special cluster discovered
by IPA. (f) represents the top diseases & biological functions discovered by
IPA.

###  Appendix L Gene Function Analysis

Here we continued discussing the application of gene embeddings in the gene
function prediction task. We further intended to predict the characters of
genes in the regulation process (as transcript factor or not) [37]. We used
MuSe-GNN to generate gene embeddings for different datasets based on an
unsupervised learning framework and utilized the gene embeddings as training
dataset to predict the function of genes based on k-NN classifier.

In this task, we evaluated the performance of MuSe-GNN based on scRNA-seq
datasets from different tissues, comparing it to the prediction results based
on raw data or PCA. On average, the gene embeddings from MuSe-GNN are the most
powerful embeddings, yielding highest accuracy in the validation set.
Moreover, the accuracy based on gene embeddings from MuSe-GNN also surpass the
score based on raw data or output of PCA in different tissues: Heart, Lung,
Liver, Thymus, PBMC, Cerebellum and Pancreas. Such result is shown in Figure
20.

Figure 20: Accuracy for Gene-TF prediction across different tissues.

###  Appendix M Datasets information

Figure 21: Information about multimodal biological datasets from Human. Figure
22: Information about multimodal biological datasets from Lemur and Mouse.

The statistics of different graphs can be found in Supplementary File 1.