Gated Recurrent Neural Networks: Exploration of LSTMs, GRUs, and MGUs

Project Report for STAT 362 Spring 2026

Abstract

Gated recurrent neural networks are a specialized type of recurrent neural network that have become popular in the field of ML, and are among the most popular types of neural networks among convolutional neural networks and attention transformers. Among gRNNs include the long short-term memory network (LSTM), gated recurrent units (GRU), and minimal gated units (MGU). Each of these networks controls information flow in and out of the network using separate “gates”, which are learned functions separate from the main information flow. The resulting computational power is a mitigation of the notorious vanishing gradients problem, which is the issue in which normal recurrent neural networks are unable to make major weight changes via gradient descent algorithms due to so many time points. By controlling the learning through gates, networks can more efficiently learn throughout a long time series. In this report, we will explore the theory behind the structure and training of gRNNs. Additionally, we will describe how these models work in practice. Our results highlight the stability of steady training and improvements seen in gRNN, while also highlighting the subtle differences and advantages of certain gRNNs over others.

1) Introduction

Recurrent neural networks (RNNs) are neural nets specialized for time-series and sequence-structured data. These models work by incorporating inputs through each time step while also having the internal hidden state “recur” onto itself, acting as if each time step was a new layer in a feed-forward neural network. RNNs are very powerful, but suffer from major drawbacks, notably the vanishing gradients problem. Gradients become very small during training due to having to propogate through large time sequences, leading to almost no weight changes. The result is a neural network that cannot effectively learn long series. To mitigate this, several modificiations to the original RNN model to incorporate “gates”, which control the flow of information through the network. We refer to these models as gated recurrent neural networks (GRNNs).

gRNNs were first formulated by [3] through the long short-term memory network (LSTM). In an LSTM, the network has two states to control general information (long-term memory) and information local in time (short-term memory) with information regulated through various gate mechanisms. The resulting model quickly became popular in its ability to mitigate the vanishing gradients problem. However, the resulting model was very complicated, so various modifications were proposed to help simplify the amount of parameters. The main one is the gated recurrent unit (GRU), which removed the two types of memory and put greater emphasis on the gating mechanics [1]. A subclass of this approach is the minimal gated unit (MGU), which is a simplification of the GRU model to include just one gate [8]. Together, these models make up some of the most popular approaches to formulating a gRNN.

First, we will look at the theory of the gRNN models. We will understand how the LSTM is formulated based on its attempt to solve the vanishing gradients problem using a special unit that ensures error signals are not diminished during training, and then see how the gated mechanisms help facilitate the resulting model. Then, we take a look at how the GRU and MGU build off of this structure. Finally, we will understand how the models are trained using backpropogation through time [6] by first exploring how this algorithm works on a vanilla RNN and then explaining the modifications necessary for gRNNs [2].

Next, we will explore the performance of the GRNN models over the vanilla RNN. This will be done using a complex EEG dataset for participants imagining motor movements [4]. We will train each of our models on this data, and compare their performance against the vanilla RNN and against each other.

2) Theory and Methods

Formulating the gRNN Model

We will first understand why gRNNs are structured the way they are. To do, we will first describe how LSTMs were first theorized, which will give us an understanding of the gating structure and why is can help mitigate vanishing gradients.

Long Short-Term Memory (LSTM) Units

Let us focus on a single recurrent cell, which takes in an input \(x^{<t>}\) and input from the previous timepoint \(h^{<t-1>}\) and outputs an output \(h^{<t>}\) according to the equation: \[ \begin{align*} h^{<t>} &= f(z^{<t>}), & z^{<t>} &= W^{hh}h^{<t-1>} + W^{xh}x^{<t>}+\vec{b} \\ \end{align*} \]

In essence, the hidden states \(h^{<t>}\) already encode what can be described as short-term memory, meaning they are only able to “remember” information for small amounts of time.

Introducing the Cell State

To incorporate aspects of long-term memory into our model, we define a new state \(c^{<t>}\) called the cell state, the goal of which is to retain information for long periods of time. \[ \begin{align*} c^{<t>} &= g(z^{c, <t>}), & z^{c, <t>} &= W^{cc}c^{<t-1>}+ W^{hc}h^{<t-1>} + W^{xc}x^{<t>}+\vec{b}^c \\ h^{<t>} &= f(z^{h, <t>}), & z^{h, <t>} &=W^{ch}c^{<t>} + W^{hh}h^{<t-1>} + W^{xh}x^{<t>}+\vec{b}^h \end{align*} \]

Formulating the CEC

Now, to have the cell state retain information, we want the weight updates for \(W^{cc}\) to not vanish as we train it. According to the formal for backpropogation derived in class, weights of \(W^{cc}\) are changed according to an error propogated through the backward phase according to:

\[\delta^{<t>} = (W^{cc})^T\delta^{<t+1>} \odot g'(z^{c, <t>})\]

For now, we assume there is one observation, so \(\delta^{<t>}, \delta^{<t+1>}, g'(z^{c, <t>})\) are vectors. Imagine we wanted the error to fully propogate back. That is, we would like \(\delta^{<t>} = \delta^{<t+1>}\). For simplicity, imagine we only have self connections, meaning \(W^{cc}\) is only non-zero on its diagonal. If we denote the \(j^{th}\) component by the subscript \(j\), then \[\delta_{j}^{<t>} = g_j'(z_j^{c, <t>})w^{cc}_{jj}\delta_{j}^{<t+1>}\] Hence, we need \(g_{j}'(z_j^{c, <t>})w^{cc}_{jj}=1\). We can then put this in a form so that we can integrate: \[ \begin{align*} \left.\frac{\partial g_{j}}{\partial z_j^{c}}\right\vert_{z_j^{c,<t>}} &= \frac{1}{w^{cc}_{jj}}\\ g_{j}(z_j^{c,<t>})&=\frac{z_j^{<t>}}{w^{cc}_{jj}} && \text{integrate both sides}\\ c^{<t>}&=c^{<t-1>} + \vec{i} \odot (W^{hc}h^{<t-1>} + W^{xc}x^{<t>}+\vec{b^c}) \end{align*} \]

Where \(\vec{i}_j = 1 / w^{cc}_{jj}\). This defines the basics for the equation of the central error carrousel (CEC), which gives recurrent self connections of weight 1 on the cell state and is the original “naive” derivation of the LSTM model given in [3]. The equation is split into two parts: one in which the previous cell state is retained, and another in which information from the incoming inputs and previous hidden states are added. By this formulation, errors within the carrousel cannot change through BPTT, ensuring information is kept. We may additionally allow for information added to the cell state to be non-linear by passing it through an activation function. This gives the set of equations

\[ \begin{align*} \tilde{c}^{<t>} &= g(z^{\tilde{c}, <t>}), & z^{\tilde{c}, <t>} &= W^{h\tilde{c}}h^{<t-1>} + W^{x\tilde{c}}x^{<t>}+\vec{b}^{\tilde{c}}\\ c^{<t>} &= c^{<t-1>}+ \vec{i} \odot \tilde{c}^{<t>}\\ h^{<t>} &= f(z^{h, <t>}), & z^{h, <t>} &=W^{ch}c^{<t>} + W^{hh}h^{<t-1>} + W^{xh}x^{<t>}+\vec{b}^h \end{align*} \]

Note \(\tilde{c}\) and \(c\) have the same shape.

Adding Output, Input, and Forget Gates

However, we still have a small problem. Since the hidden state depends on the previous cell state, the cell state also depends on previous cell states according to the weight matrix \(W^{ch}\), which is not necessarily diagonal with unit weights. Hence our argument of constant propagating error does not hold. We could use the same reformulation as before, but that would just kick the problem down the line, as the cell state would still be dependent on previous cell states according to \(W^{h\tilde{c}}\). Instead, we allow for the hidden state to use some information, but not all, from the cell state. This is done having \(h\) be a multiplicative function of the cell state, in which the multiplicative weights are learned by some output \(o\):

\[ \begin{align*} \tilde{c}^{<t>} &= g(z^{\tilde{c}, <t>}), & z^{\tilde{c}, <t>} &= W^{h\tilde{c}}h^{<t-1>} + W^{x\tilde{c}}x^{<t>}+\vec{b}^{\tilde{c}}\\ o^{<t>} &= \sigma (z^{o, <t>}), & z^{o, <t>} &=W^{ho}h^{<t-1>} + W^{xo}x^{<t>}+\vec{b}^o\\ c^{<t>} &= c^{<t-1>}+ \vec{i} \odot \tilde{c}^{<t>}\\ \tilde{h}^{<t>} &= f(c^{<t>})\\ h^{<t>} &= o^{<t>} \odot \tilde{h}^{<t>} \end{align*} \]

\(o\) is called the output gate, and determines what information from the “output” \(\tilde{h}\) is actually outputed. Here, we explicity use the sigmoid activation function \(\sigma\). This ensures that the values for \(\vec{i}\) are kept between 0 and 1: 1 meaning the information is completley kept, and 0 meaning the information is completely ignored. Now, note again how the cell state is made up of two different parts: one that remembers information, and one that learns information. For the latter, one could improve the model by allowing the model to learn what it should learn. This is done by letting \(\vec{i}\) change in time by being a function of the inputs and previous hidden state: \[ \begin{align*} \tilde{c}^{<t>} &= g(z^{\tilde{c}, <t>}), & z^{\tilde{c}, <t>} &= W^{h\tilde{c}}h^{<t-1>} + W^{x\tilde{c}}x^{<t>}+\vec{b}^{\tilde{c}}\\ o^{<t>} &= \sigma (z^{o, <t>}), & z^{o, <t>} &=W^{ho}h^{<t-1>} + W^{xo}x^{<t>}+\vec{b}^o\\ \vec{i}^{<t>} &= \sigma(z^{i, <t>}), & z^{i, <t>} &= W^{hi}h^{<t-1>} + W^{xi}x^{<t>}+\vec{b}^{i}\\ c^{<t>} &= c^{<t-1>} + \vec{i}^{<t>} \odot \tilde{c}^{<t>}\\ \tilde{h}^{<t>} &= f(c^{<t>})\\ h^{<t>} &= o^{<t>} \odot \tilde{h}^{<t>} \end{align*} \]

\(\vec{i}\) is referred to as the input gate, and controls what information from the “input” \(\tilde{c}\) is actually inputted. Here, we again explicity use the sigmoid activation function \(\sigma\) to keep the weights between 0 and 1 to cover the spectrum of not learning and fully learning. Similarly, we can also let the model learn what it should remember, or in other words, learn what it should forget. Currently, the cell state has constant recurrent self-connections of weight 1. If the weights were between 0 and 1, it would allow the model to forget some information to prevent it from remembering uncessary information. This motivates a time-varying parameter \(\vec{k}\) similar to that of \(\vec{i}\):

\[ \begin{align} \tilde{c}^{<t>} &= g(z^{\tilde{c}, <t>}), & z^{\tilde{c}, <t>} &= W^{h\tilde{c}}h^{<t-1>} + W^{x\tilde{c}}x^{<t>}+\vec{b}^{\tilde{c}}\\ o^{<t>} &= \sigma (z^{o, <t>}), & z^{o, <t>} &=W^{ho}h^{<t-1>} + W^{xo}x^{<t>}+\vec{b}^o\\ \vec{i}^{<t>} &= \sigma(z^{i, <t>}), & z^{i, <t>} &= W^{hi}h^{<t-1>} + W^{xi}x^{<t>}+\vec{b}^{i}\\ \vec{k}^{<t>} &= \sigma(z^{k, <t>}), & z^{k, <t>} &= W^{hk}h^{<t-1>} + W^{xk}x^{<t>}+\vec{b}^{k}\\ c^{<t>} &= \vec{k}^{<t>} \odot c^{<t-1>} + \vec{i}^{<t>} \odot \tilde{c}^{<t>}\\ \tilde{h}^{<t>} &= f(c^{<t>})\\ h^{<t>} &= o^{<t>} \odot \tilde{h}^{<t>} \end{align} \]

\(\vec{k}\) is then referred to as the forget gate. However if the forget gate does not initially start as all ones (an instead as a value between 0 and 1), it would cause gradients to vanish as per are argument above, definiting the purpose of the cell state. In practice, it is recommended to initialize \(\vec{b}^{k}\) to be all ones (or some other large value), which forces \(\vec{k}\) to initialize close to all ones. Hence, by lowering the values of \(\vec{k}\), the network ``learns to forget”.

In practice, \(f\) and \(g\) are usually taken as \(\tanh\) to prevent values of the cell state and hidden state from “blowing up” by keeping the values constrained between -1 and 1. The final set of equations above then makeup what is known as the forward phase of an LSTM unit.

This unit has several variations, including peephole connections, which allow for the cell state to influence each of the 3 gates, and coupled LSTMs, which combine the forget and input gates under the assumption that \(\vec{i} = 1 - \vec{k}\). However, we do not formulate those variations here.

Gated Recurrent Units (GRUs)

As you can see, LSTMs are made up of a considerable number trainable parameters. To mitigate this, efforts have been made to reduce the number of moving parts in the model. The most popular way has been the Gated Recurrent Unit (GRU), which attempts to simplify an LSTM unit by combining the cell state and hidden state into one state, given by \(h\). The idea is that we hope we can encode both long term and short term information within one state. Looking at the LSTM equations, this change already simplifies our model by eliminating the need for an explicit output gate:

\[ \begin{align*} \tilde{h}^{<t>} &= g(z^{\tilde{h}, <t>}), & z^{\tilde{h}, <t>} &= W^{h\tilde{h}}h^{<t-1>} + W^{x\tilde{h}}x^{<t>}+\vec{b}^{\tilde{h}}\\ \vec{i}^{<t>} &= \sigma(z^{i, <t>}), & z^{i, <t>} &= W^{hi}h^{<t-1>} + W^{xi}x^{<t>}+\vec{b}^{i}\\ \vec{k}^{<t>} &= \sigma(z^{k, <t>}), & z^{k, <t>} &= W^{hk}h^{<t-1>} + W^{xk}x^{<t>}+\vec{b}^{k}\\ h^{<t>} &= \vec{k}^{<t>} \odot h^{<t-1>} + \vec{i}^{<t>} \odot \tilde{h}^{<t>} \end{align*} \]

In essence, this is the biggest change from an LSTM. However, the GRU also makes some other minor changes. First, the input gate is applied to \(h\) before it is transformed into \(\tilde{h}\), allowing it to better control what information is ignored or “reset”, and hence renamed the reset gate and denoted by \(\vec{r}\) (Although, strictly speaking, it is \(1-\vec{r}\) which controls what is “reset”)

\[ \begin{align*} \vec{r}^{<t>} &= \sigma(z^{r, <t>}), & z^{r, <t>} &= W^{hr}h^{<t-1>} + W^{xr}x^{<t>}+\vec{b}^{r}\\ \tilde{h}^{<t>} &= g(z^{\tilde{h}, <t>}), & z^{\tilde{h}, <t>} &= W^{h\tilde{h}}(\vec{r}^{<t>} \odot h^{<t-1>}) + W^{x\tilde{h}}x^{<t>}+\vec{b}^{\tilde{h}}\\ \vec{k}^{<t>} &= \sigma(z^{k, <t>}), & z^{k, <t>} &= W^{hk}h^{<t-1>} + W^{xk}x^{<t>}+\vec{b}^{k}\\ h^{<t>} &= \vec{k}^{<t>} \odot h^{<t-1>} + \tilde{h}^{<t>} \end{align*} \]

Additionally, the forget gate is modified to instead control what is remembered or “updated” from the new information, and hence renamed the updated gate which we denote by \(\vec{u}\). This gate also controls what old information is forgotten (like the original forget gate) by using \(1-\vec{u}\). The resulting equations gives the forward pass for the GRU:

\[ \begin{align*} \vec{r}^{<t>} &= \sigma(z^{r, <t>}), & z^{r, <t>} &= W^{hr}h^{<t-1>} + W^{xr}x^{<t>}+\vec{b}^{r}\\ \tilde{h}^{<t>} &= g(z^{\tilde{h}, <t>}), & z^{\tilde{h}, <t>} &= W^{h\tilde{h}}(\vec{r}^{<t>} \odot h^{<t-1>}) + W^{x\tilde{h}}x^{<t>}+\vec{b}^{\tilde{h}}\\ \vec{u}^{<t>} &= \sigma(z^{u, <t>}), & z^{u, <t>} &= W^{hu}h^{<t-1>} + W^{xu}x^{<t>}+\vec{b}^{u}\\ h^{<t>} &= (1 - \vec{u}^{<t>}) \odot h^{<t-1>} \\ &+ \vec{u}^{<t>} \odot\tilde{h}^{<t>} \end{align*} \]

Minimal Gated Units (MGUs)

The gated mechanism has been further simplified into the MGU. This unit does away with the reset gate, and instead assumes that the information that is not reset is also the information updated. That is, \(\vec{u}^{<t>} = \vec{r}^{<t>}\), which gives the forward pass equations

\[ \begin{align} \vec{u}^{<t>} &= \sigma(z^{u, <t>}), & z^{u, <t>} &= W^{hu}h^{<t-1>} + W^{xu}x^{<t>}+\vec{b}^{u}\\ \tilde{h}^{<t>} &= g(z^{\tilde{h}, <t>}), & z^{\tilde{h}, <t>} &= W^{h\tilde{h}}(\vec{u}^{<t>} \odot h^{<t-1>}) + W^{x\tilde{h}}x^{<t>}+\vec{b}^{\tilde{h}}\\ h^{<t>} &= (1 - \vec{u}^{<t>}) \odot h^{<t-1>} \\ &+ \vec{u}^{<t>} \odot\tilde{h}^{<t>} \end{align} \]

Most literature referes to \(\vec{u}\) as the forget gate in the context of MGUs since it most closely mimicks the role of the forget gate \(\vec{k}\) from the LSTM.

Backpropogation Through Time (BPTT)

The algorithm to train gRNNs is backpropogation through time (BPPT) [7]. Backpropogation through time uses the same algorithm as normal backpropogation for feedforward neural networks, except each time point is “unfolded” through time and treated as a separate layer. That is, a recurrent neural network of 10 time points can be thought of as a feedforward neural network with 10 layers constrained to have the same weight matrices. Each iteration of backpropogation then leads to an weight update for the same weight matrix.

We start with a vanilla recurrent neural network with time points \(0 \leq t \leq T\) such that for \(t \neq 0\) \[ \begin{align*} h^{<t>} &= f(z^{<t>}), & z^{<t>} &= W^{hh}h^{<t-1>} + W^{xh}x^{<t>}+\vec{b} \\ \end{align*} \]

Consider some loss function \(L\). We wish to determine the gradients to update \(W^{hh}, W^{xh},\) and \(\vec{b}\). That is, we wish to calculate \(\partial L / \partial W^{hh}, \partial L / \partial W^{xh},\) and \(\partial L / \partial \vec{b}\). Consider first \(\partial L / \partial W^{hh}\), which we expand via the chainrule: \[ \begin{align*} \frac{\partial L}{\partial W^{hh}} &= \frac{\partial L}{\partial h^{<T>}} \cdot \frac{\partial h^{<T>}}{\partial z^{<T>}} \cdot \frac{\partial z^{<T>}}{\partial W^{hh}} = \left[\frac{\partial L}{\partial h^{<T>}} \odot f'(z^{<T>})\right](h^{<T-1>})^T &&\text{[7], eq. (69)} \end{align*} \] This derivation is remnant of that of normal backpropogation. We define an error signal \(\delta^{<T>} = \frac{\partial L}{\partial h^{<T>}} \odot f'(z^{<T>})\). Similary to backpropogation, we do not need to terminate the chain rule at \(z^{T}\). Instead, we can expand to the previous layer to obtain a different expression. Unlike normal backpropogation, however, this results in a gradient for the same weight matrix: \[ \begin{align*} \frac{\partial L}{\partial W^{hh}} &= \frac{\partial L}{\partial h^{<T>}} \cdot \frac{\partial h^{<T>}}{\partial z^{<T>}} \cdot \frac{\partial z^{<T>}}{\partial h^{<T-1>}} \cdot \frac{\partial h^{<T-1>}}{\partial z^{<T-1>}} \cdot \frac{\partial z^{<T-1>}}{\partial W^{hh}}\\ &= \delta^{<T>} \cdot \frac{\partial z^{<T>}}{\partial h^{<T-1}}\cdot \frac{\partial h^{<T-1>}}{\partial z^{<T-1}}\cdot \frac{\partial z^{<T-1>}}{\partial W^{hh}}\\ &= (W^{hh})^T\delta^{<T>} \cdot \frac{\partial h^{<T-1>}}{\partial z^{<T-1}}\cdot \frac{\partial z^{<T-1>}}{\partial W^{hh}} &&\text{[7], eq. (69)}\\ &= \left[(W^{hh})^T\delta^{<T>} \odot f'(z^{<T-1>}) \right ](h^{<T-2>})^T &&\text{[7], eq. (69)}\\ \end{align*} \]

Again, like normal backpropogation, we define the next error signal to be \(\delta^{<T-1>} = (W^{hh})^T\delta^{<T>} \odot f'(z^{<T-1>})\). It is clear then that this pattern will repeat as normal. However, we again emphasize that each time step results in a gradient for the same weight matrix and hence, the same weight matrix is updated each time and used in the next weight update. Let \(W^{hh, <t>}\) be the updated weight matrix \(W^{hh}\) after using the gradient obtained from time point \(t\). At the end of the epoch, we take \(W^{hh} = W^{hh, <0>}\). We then arrive the folllowing updates for \(W^{hh}\), for some learning rate \(\eta\):

\[ \begin{align*} \delta^{<t>} &= (W^{hh, <t+1>})^T\delta^{<t+1>} \odot f'(z^{<t>}) \\ W^{hh, <t>} &= W^{hh, <t+1>} - \eta \delta^{<t>}(h^{<t-1>})^T \end{align*} \]

It is also not hard to see, then, that \(x\) and \(\vec{b}\) follow the same mathematical derivation, with the special note that \(\frac{\partial z^{<t>}}{\partial \vec{b}} = \vec{1}\) for all time points. Following similar notation from before, this gives

\[ \begin{align*} \delta^{<t>} &= (W^{hh, <t+1>})^T\delta^{<t+1>} \odot f'(z^{<t>}) \\ W^{hh, <t>} &= W^{hh, <t+1>} - \eta \delta^{<t>}(h^{<t-1>})^T\\ W^{xh, <t>} &= W^{xh, <t+1>} - \eta \delta^{<t>}(x^{<t-1>})^T\\ \vec{b}^{<t>} &= \vec{b}^{<t+1>} - \eta \delta^{<t>}\vec{1} \end{align*} \]

BPTT For Gated RNNs

The derivation for BPTT for gRNNs is very similar. Let’s consider the LSTM. Note that the structure of the block inputs and the three gates are very similar to that of the vanilla RNN:

\[ \begin{align} \tilde{c}^{<t>} &= g(z^{\tilde{c}, <t>}), & z^{\tilde{c}, <t>} &= W^{h\tilde{c}}h^{<t-1>} + W^{x\tilde{c}}x^{<t>}+\vec{b}^{\tilde{c}}\\ o^{<t>} &= \sigma (z^{o, <t>}), & z^{o, <t>} &=W^{ho}h^{<t-1>} + W^{xo}x^{<t>}+\vec{b}^o\\ \vec{i}^{<t>} &= \sigma(z^{i, <t>}), & z^{i, <t>} &= W^{hi}h^{<t-1>} + W^{xi}x^{<t>}+\vec{b}^{i}\\ \vec{k}^{<t>} &= \sigma(z^{k, <t>}), & z^{k, <t>} &= W^{hk}h^{<t-1>} + W^{xk}x^{<t>}+\vec{b}^{k}\\ \end{align} \]

Hence, the equations to update each of the weight matrices is similar. The difference, however, is that each error update must additionally pass through the output equations

\[ \begin{align*} c^{<t>} &= \vec{k}^{<t>} \odot c^{<t-1>} + \vec{i}^{<t>} \odot \tilde{c}^{<t>}\\ \tilde{h}^{<t>} &= f(c^{<t>})\\ h^{<t>} &= o^{<t>} \odot \tilde{h}^{<t>} \end{align*} \]

which means there are additional multiplications at every step. Consider, for example, updating \(W^{hi}\). If we derive the backpropogation equations, we get:

\[ \begin{align*} \frac{\partial L}{\partial W^{hi}} &= \frac{\partial L}{\partial h^{<T>}} \cdot \frac{\partial h^{<T>}}{\partial \tilde{h}^{<T>}} \cdot \frac{\partial \tilde{h}^{<T>}}{\partial c^{<T>}} \cdot \frac{\partial c^{<T>}}{\partial \vec{i}^{<T>}} \cdot \frac{\partial \vec{i}^{<T>}}{\partial z^{i, <T>}} \cdot \frac{\partial z^{i, <T>}}{\partial W^{hi}}\\ &= \frac{\partial L}{\partial h^{<T>}} \odot o^{<T>} \odot f'(c^{<T>}) \odot \tilde{c}^{<T>} \cdot \frac{\partial \vec{i}^{<T>}}{\partial z^{i, <T>}} \cdot \frac{\partial z^{i, <T>}}{\partial W^{hi}}\\ &= o^{<T>} \odot f'(c^{<T>}) \odot \tilde{c}^{<T>} \odot \left[\frac{\partial L}{\partial h^{<T>}} \odot \sigma'(z^{i, <T>})\right] \cdot \frac{\partial z^{i, <T>}}{\partial W^{hi}}\\ \end{align*} \]

The part in brackets is the familar error signal \(\delta^{i,<T>}\), while the parts to the left are new multiplicative factors due to the extra output layers which we call \(\lambda^{<i,T>}\). In the next time step, its easy to see that these multiplicative factors arise again. In general, our BPTT equation for \(W^{hi}\) stays the same, with the inclusion of the added factors:

\[ \begin{align*} \delta^{i,<t>} &= \lambda^{<i,t+1>}(W^{hi, <t+1>})^T\delta^{i, <t+1>} \odot \sigma'(z^{<t>}) \\ W^{hi, <t>} &= W^{hi, <t+1>} - \eta \delta^{i,<t>}(h^{<t-1>})^T \end{align*} \]

The weight updates for other parameters, as well as for the other models GRU and MGU, follow similar modifications using the same logic.

3) Data

import os
import time
import numpy as np
import pandas as pd

import mne
from mne.preprocessing import ICA
from sklearn.preprocessing import LabelEncoder
from keras.utils import to_categorical
from sklearn.model_selection import train_test_split

import tensorflow as tf
from keras.models import Sequential
from keras.layers import Layer, RNN, SimpleRNN, LSTM, GRU, InputLayer, Dense, Dropout
from keras.callbacks import EarlyStopping

from sklearn.metrics import accuracy_score, confusion_matrix

from scipy.signal import spectrogram
import matplotlib.pyplot as plt

mne.set_log_level('ERROR')

def eeg_data_preprocess(file_path):
    raw = mne.io.read_raw_gdf(file_path,preload=True)
    raw.set_channel_types({
        raw.ch_names[-3]: 'eog',
        raw.ch_names[-2]: 'eog',
        raw.ch_names[-1]: 'eog'
    })
    raw.rename_channels(lambda x: x.replace('EEG-', ''))
    montage = mne.channels.make_standard_montage('standard_1020')
    raw.set_montage(
        montage,
        match_case=False,
        on_missing='ignore'
    )

    raw.filter(1., 40.)
    picks_eeg = mne.pick_types(
        raw.info,
        eeg=True,
        eog=False
    )
    ica = ICA(
        n_components=20,
        method='fastica',
        random_state=42
    )

    ica.fit(raw, picks=picks_eeg)
    eog_indices, eog_scores = ica.find_bads_eog(raw)
    ica.exclude = eog_indices
    raw_clean = ica.apply(raw.copy())
    raw_clean.pick_types(eeg=True, eog=False)
    raw_clean.filter(8., 30.)

    events, event_id = mne.events_from_annotations(raw_clean)
    event_dict = {
        'left': event_id['769'],
        'right': event_id['770'],
        'foot': event_id['771'],
        'tongue': event_id['772']
    }
    epochs = mne.Epochs(
        raw_clean,
        events,
        event_id=event_dict,
        tmin=0.5,
        tmax=4.5,
        baseline=None,
        preload=True
    )

    X = epochs.get_data()
    y_ = epochs.events[:, -1]

    le = LabelEncoder()
    y = le.fit_transform(y_)

    X = (X - X.mean(axis=2, keepdims=True)) / (X.std(axis=2, keepdims=True))

    X = X.transpose(0, 2, 1)
    return X, y, ica, epochs

X, y, ica, epochs = eeg_data_preprocess('BCICIV_2a_gdf/A01T.gdf')

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.3, stratify = y, random_state = 1)
y_train_dummy = to_categorical(y_train, 4)

The data comes from a motor imagery experiment in which test subjects imaged making a specific motor movement when cued on a screen [5]. Depending on the cue, participants wer asked to imagine moving either their left hand, right hand, foot, or tongue. Electroencephalogram (EEG) data was collected from 22 electrodes spatially distributed acorss their scalp, along with Electrooculogram (EOG) data to measure the movements of the participants eyes. The EOG data was not used in classification and only used to help clean the EEG data from signals due to eye movements. These signals were measured for 4 seconds after each cue at a sampling frequency of 250Hz, resulting in 1001 time points.

This data was selected for due to its large complexity overtime. EEG signals are an important measurement for computational psychology and neuroscience, and undertsanding these signals are crucial to understanding how computations are spatially organized in the brain. Each EEG signal is likely to be specific to each of the four movements due to the activity in the sensiromotor cortex being topographically organized. RNNs therefore can be hugely beneficial for helping learn the large time-series data. However, learning that specificity across 22 channels 1001 time points is challenging. Hence, we would expect the vanishing gradient problem to play a role and hinder learning of a vanilla RNN. This makes the data suitable for comparing the benefits of gRNNs.

The data was cleaned from its original form. Firstly only 1 out of 9 sessions from the original experiment was used in order to make training faster. Each of the four movements were isolated from this session, resulting in 288 total observations. Signals were then filtered through a standard bandpass filter to remove irrelevant frequencies. As previously stated, each observation was also cleaned to remove irrelevant signals due to eye movements using Independent Component Analysis (ICA) [4]. The resulting data was then normalized in order to make signals comparable across trials. This resulted in 288 observations of (1001, 22) dimensional data. The data was then split using a 70%-30% train test split, and the observations were one-hot-encoded.

fs = 250

signal = X[0, :, 0]
time = np.arange(len(signal)) / fs

plt.figure(figsize=(12, 4))

plt.plot(time, signal, linewidth=1)

plt.title('Example EEG Signal\nObservation 1, Channel 1', fontsize=16)
plt.xlabel('Time (s)', fontsize=12)
plt.ylabel('Normalized Amplitude', fontsize=12)

plt.xlim(time[0], time[-1])

plt.tight_layout()
plt.show()

Here are the average spectrograms for each resulting channel grouped by imagined movement. Channels are spatially organized in the same manner that the experiment describes across the participants scalp.

fs = 250

label_specs = {}

labels = ['Left Hand', 'Right Hand', 'Foot', 'Tongue']

for label in range(4):
    X_label = X[y == label]
    avg_specs = []

    for ch in range(22):
        specs = []

        for obs in X_label:
            f, t, Sxx = spectrogram(
                obs[:, ch],
                fs=fs,
                nperseg=128,
                noverlap=64
            )

            freq_mask = (f >= 7.5) & (f <= 40)

            specs.append(Sxx[freq_mask])

        avg_specs.append(np.mean(specs, axis=0))

    label_specs[label] = avg_specs

layout = [
    [0],
    [1, 2, 3, 4, 5],
    [6, 7, 8, 9, 10, 11, 12],
    [13, 14, 15, 16, 17],
    [18, 19, 20],
    [21]
]

fig, axes = plt.subplots(2, 2, figsize=(12, 12))
axes = axes.flatten()

fig.suptitle(
    'Average Spectrogram Over Space For Each Motor Imagery Movement',
    fontsize=22,
    y=0.98
)

for label, label_name in enumerate(labels):
    ax = axes[label]

    for row_idx, row in enumerate(layout):
        n = len(row)

        x_positions = np.linspace(
            -((n - 1) * 0.9) / 2,
            ((n - 1) * 0.9) / 2,
            n
        )

        y_base = (len(layout) - row_idx - 1) * 42

        for xpos, ch in zip(x_positions, row):
            Sxx = label_specs[label][ch]

            Sxx = np.log1p(Sxx)
            Sxx = (Sxx - Sxx.min()) / (Sxx.max() - Sxx.min())

            extent = [
                xpos - 0.38,
                xpos + 0.38,
                y_base + 7.5,
                y_base + 40
            ]

            ax.imshow(
                Sxx,
                aspect='auto',
                origin='lower',
                extent=extent
            )

    ax.set_title(label_name, fontsize=16)
    ax.set_xlim(-3.2, 3.2)
    ax.set_ylim(0, 255)
    ax.set_xticks([])
    ax.set_yticks([])

    for spine in ax.spines.values():
        spine.set_visible(False)

plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()

Please note that the EEG signals are still exceptionally noisy. It is likely that much more preprocessing is needed, and more complex network architectures should be utilized (for example, a CNN to capture spatial information). However, we are only interested in comparing how these models perform during training to understand how they combat the vanishing gradients problem. We should therefore expect accuracies to be fairly low, but this should not hinder our main results.

4) Results and Discussion

To exemplify the power of gRNNs, we will construct 4 similar neural networks to help classify the type of motor imagery during each trial. The first neural network is uses a vanilla RNN (vRNN) as the control, while the other use an LSTM, GRU, and MGU respectively. The EEG data is likely very complicated, and a robust model would take into account many more components to achieve better prediction. However, for the sake of comparison, we will be constructing fairly simple RNN networks to fit to the data. We do not want them too simple, though, as to create meaningless prediction. Hence, for each of the four networks, we will follow this same general layout:

1. An input layer taking in (10001, 25) dimensional data.
2. 2 recurrent blocks (either vRNN, LSTM, GRU, or MGU), which take in 500 and 100 dimensional inputs, respectively.
3. A dense layer of 20 inputs to help compress the outputs from (2).
4. Dropout of 0.15 on layers from (2), (3) to help prevent large overfitting.
5. An output layer of size 4 (for left, right, foot, tongue).

Additionally, since the original GRNN models use the \(\tanh\) function, we will be keeping that default activation in (2) for the sake of comparison. Note as well the while LSTMs and GRUs are included in keras, MGUs are not. Hence, we have coded a simple MGUcell function to help compare against the other models.

class MGUCell(Layer):
    def __init__(self, units, **kwargs):
        super().__init__(**kwargs)
        self.units = units
        self.state_size = units

    def build(self, input_shape):
        d = input_shape[-1]

        # W^{xu}, W^{hu}
        self.W_xu = self.add_weight(shape=(d, self.units))
        self.W_hu = self.add_weight(shape=(self.units, self.units))
        self.b_u  = self.add_weight(shape=(self.units,), initializer="ones")

        # W^{x\tilde{h}}, W^{h\tilde{h}}
        self.W_xh = self.add_weight(shape=(d, self.units))
        self.W_hh = self.add_weight(shape=(self.units, self.units))
        self.b_h  = self.add_weight(shape=(self.units,), initializer="zeros")

    def call(self, inputs, states):
        h_prev = states[0]

        # forget gate
        u = tf.sigmoid(
            tf.matmul(h_prev, self.W_hu) +
            tf.matmul(inputs, self.W_xu) +
            self.b_u
        )

        uh = u * h_prev

        # recurrent layer
        z_h = (
            tf.matmul(uh, self.W_hh) +
            tf.matmul(inputs, self.W_xh) +
            self.b_h
        )

        # output
        h_tilde = tf.tanh(z_h)
        h = (1 - u) * h_prev + u * h_tilde

        return h, [h]

model_labels = ['vRNN', 'LSTM', 'GRU', 'MGU']

rnn = Sequential()
rnn.add(InputLayer(shape = (10001,22)))
rnn.add(SimpleRNN(500, return_sequences=True))
rnn.add(Dropout(0.15))
rnn.add(SimpleRNN(100))
rnn.add(Dropout(0.15))
rnn.add(Dense(20, activation = 'relu'))
rnn.add(Dropout(0.15))
rnn.add(Dense(4, activation = 'softmax'))

lstm = Sequential()
lstm.add(InputLayer(shape = (10001,22)))
lstm.add(LSTM(500, return_sequences=True))
lstm.add(Dropout(0.15))
lstm.add(LSTM(100))
lstm.add(Dropout(0.15))
lstm.add(Dense(20, activation = 'relu'))
lstm.add(Dropout(0.15))
lstm.add(Dense(4, activation = 'softmax'))

gru = Sequential()
gru.add(InputLayer(shape = (10001,22)))
gru.add(GRU(500, return_sequences=True))
gru.add(Dropout(0.15))
gru.add(GRU(100))
gru.add(Dropout(0.15))
gru.add(Dense(20, activation = 'relu'))
gru.add(Dropout(0.15))
gru.add(Dense(4, activation = 'softmax'))

mgu = Sequential()
mgu.add(InputLayer(shape = (10001,22)))
mgu.add(RNN(MGUCell(500), return_sequences=True))
mgu.add(Dropout(0.15))
mgu.add(RNN(MGUCell(100)))
mgu.add(Dropout(0.15))
mgu.add(Dense(20, activation = 'relu'))
mgu.add(Dropout(0.15))
mgu.add(Dense(4, activation = 'softmax'))

models = [rnn, lstm, gru, mgu]
params = [rnn.count_params(), lstm.count_params(), gru.count_params(), mgu.count_params()]

We will also conduct early stopping using a epoch window of length 10 and validation set size of 10% to help prevent overfitting. Hence, we will run fo 50 epochs, but most models stop at epoch 10.

os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'

model_history = []
runtime = []
early = EarlyStopping('val_loss', patience = 10)
for model in models:
    model.compile(loss = 'categorical_crossentropy', optimizer = 'adam', metrics = ['accuracy'],)

    start = time.time()
    history = model.fit(X_train, y_train_dummy, 
                        epochs = 50, 
                        validation_split = 0.1,
                        callbacks = early,
                        verbose=0)
    end = time.time()

    model_history.append(history)
    runtime.append(np.round(end - start, 2))

The resulting sizes and runtimes of our models are given below:

model_comparisons = pd.DataFrame({
    'Model': model_labels,
    'Parameters': params,
    'Runtime (s)': runtime
})

model_comparisons['Parameters'] = model_comparisons['Parameters'].map('{:,}'.format)
model_comparisons['Runtime (s)'] = model_comparisons['Runtime (s)'].map('{:.3f}'.format)

import pandas as pd

model_comparisons = pd.DataFrame({
    'Model': model_labels,
    'Parameters': params,
    'Runtime (s)': runtime
})

model_comparisons['Parameters'] = model_comparisons['Parameters'].map('{:,}'.format)
model_comparisons['Runtime (s)'] = model_comparisons['Runtime (s)'].map('{:.3f}'.format)

styled_table = (
    model_comparisons.style
    .hide(axis='index')
    .set_caption("Model Comparison Summary")
    .set_table_styles([
        {
            'selector': 'table',
            'props': [
                ('width', '600px'),
                ('margin', '0 auto'),
                ('border-collapse', 'collapse')
            ]
        },
        {
            'selector': 'caption',
            'props': [
                ('font-size', '16px'),
                ('font-weight', 'bold'),
                ('text-align', 'center'),
                ('padding', '10px'),
                ('white-space', 'nowrap')
            ]
        },
        {
            'selector': 'th',
            'props': [
                ('font-size', '13px'),
                ('text-align', 'center'),
                ('font-weight', 'bold'),
                ('padding', '8px')
            ]
        },
        {
            'selector': 'td',
            'props': [
                ('text-align', 'center'),
                ('padding', '8px')
            ]
        }
    ])
)

styled_table

Table 1: Model Comparison Summary

Model	Parameters	Runtime (s)
vRNN	323,704	216.100
LSTM	1,288,504	1086.510
GRU	968,704	1190.850
MGU	645,304	516.020

Note that both the parameter size and runtime of the LSTM is exceptionally larger than a vRNN. However, the GRU and MGU are much smaller than the vRNN, highlighting their role as a “simplified” LSTM. Interestingly, while the MGU runtime is comparable to that of a vRNN, the GRU runtime is even longer than an LSTM. This may highlight that the GRUs structure is less well suited for training than an LSTM. It also highlights that the GRU is not necessarily better than an LSTM: it is just different. The relative strengths and weakeness of gRNNs are likely heavily dependent on the data and task involved.

Next, we look at how the models perform over BPTT:

colors = ['#9c179e','#0d0887','#ed7953','#f0f921']

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

for i, history in enumerate(model_history):
    axes[0].plot(history.epoch,
                 history.history['loss'],
                 label=model_labels[i],
                 color=colors[i])
    axes[1].plot(history.epoch,
                 history.history['accuracy'],
                 label=model_labels[i],
                 color=colors[i])

axes[0].set_title('Training Loss Over Epochs')
axes[0].set_xlabel('Epoch')
axes[0].set_ylabel('Loss')
axes[0].legend()
    

axes[1].set_title('Training Accuracy Over Epochs')
axes[1].set_xlabel('Epoch')
axes[1].set_ylabel('Accuracy')
axes[1].legend()

plt.tight_layout()
plt.show()

The vRNN struggles to derease in error or increase in accuracy throughout the duration of 16 epochs. However, each gRNN is able to steadily decrease in error and relatively steadily increase in accuracy and terminates in 10 or under epochs.

This showcases two things. One, gRNNs are able to steadily decrease the loss function by effectively changing their weights according to BPTT, while vRNN struggles to make any notable improvements. This is likely due to the vanishing gradients problem: with so many time steps, vRNN’s gradients get so small that it struggles to make notable parameter improvements from epoch to epoch. Hence, gRNNs architecture is able to mitigate the effects of the vanishing gradients.

Secondly, by mitigating these effects, gRNNs likely also uitilize more information by being able to learn throughout all timepoints. This extra information allows it to learn quicker and terminate sooner than vRNNs, while also more effectively increase in accuracy.

Finally, we take a look at how each model performs on the test data:

accuracies = []
conf_matrices = []

for label, model in zip(model_labels, models):
    y_preds_proba = model.predict(X_test, verbose=0)
    y_preds = np.argmax(y_preds_proba, axis=1)

    acc = accuracy_score(y_test, y_preds)
    cm = confusion_matrix(y_test, y_preds)

    accuracies.append(acc)
    conf_matrices.append(cm)

fig = plt.figure(figsize=(16, 8))

left_grid = fig.add_gridspec(
    2, 2,
    left=0.1,
    right=0.5,
    top=0.92,
    bottom=0.10,
    wspace=0.1,
    hspace=0.15
)

for idx, (cm, label) in enumerate(zip(conf_matrices, model_labels)):
    ax = fig.add_subplot(left_grid[idx // 2, idx % 2])

    im = ax.imshow(cm, cmap='Greys')

    ax.set_title(label)
    ax.set_xlabel("Predicted")
    ax.set_ylabel("True")

    ax.set_xticks([])
    ax.set_yticks([])

ax_bar = fig.add_axes([0.55, 0.20, 0.25, 0.60])

bars = ax_bar.bar(model_labels, accuracies)

ax_bar.set_ylabel("Accuracy")
ax_bar.set_ylim(0, 1)
ax_bar.set_title("Model Test Accuracies")

ax_bar.axhline(
    y=0.25,
    color='black',
    linestyle='--',
    linewidth=1
)

ax_bar.text(
    len(model_labels) - 0.5,
    0.26,
    "Chance",
    ha='right',
    va='bottom'
)

for i, acc in enumerate(accuracies):
    ax_bar.text(i, acc + 0.01, f"{acc:.3f}", ha='center')

plt.show()

Again, we note that this data is complicated, so it is not surprising that the test accuracy is much lower. Each of the models are likely acceptable to overfitting from so many weight updates through time, and since the data size is so small. However, of note is that since there are 4 classes, an accuracy of 0.25 signifies classification via chance. The vRNN just misses this threshold, while each gRNN is able to surpass it. Thus, while its numerical improvement on the test data is much smaller than the training data, its effective improvement is extremely powerful compared to the vRNN.

It is also worth noting the comparable accuracy and efficiency of each gRNN. One model does not seem to vastly outperform the other in any metric, except notably time, in which the MGU surpasses that of the LSTM and GRU. Therefore, while the performance is likely dependent on the scenario, it is worth noting that this shows the MGU may be more efficient at using the gating architecture to achieve better performance. It also highlights the original claim from the authors of the MGU [8], which is that the forget gate, or the gate that controls what information is added to each hidden state, is likely the most important one.

5) Conclusion

Gated RNNs are powerful architectures. By uitilizing gate mechanisms to control learning of information throughout time, they are computationally more efficient than vanilla RNN architectures for complex sequence problems over many time points. The LSTM architecture, the most common gRNN, divides information between the short-term hidden state and long-tern cell state, and controls how information is entered and updated in the cell state via the input and forget gates, respectively, along with what information flows into the hidden state via the output gate. The resulting performance of LSTM is superior to that of the vRNN during the complex EEG motor imagery classification task, and is able to learn steadily over fewer epochs. However, its large parameters size impedes runtime compared to vRNNs. GRUs attempt to simplify the LSTM by putting greater emphases on the gates and combines short-term and long-term information. However, in practice GRUs and LSTMs perform very similarly, and which is better is more determined on a case-by-case bases. A second improvement, known as the MGU, uses only a single forget gate. It’s performance matches that of the LSTM and GRU while being significantly quicker, highlighting not only its computational efficiency but also the importance of the forget gate among the others.

References

[1] K. Cho et al., “Learning Phrase Representations using RNN Encoder–Decoder for Statistical Machine Translation,” in Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP), A. Moschitti, B. Pang, and W. Daelemans, Eds., Doha, Qatar: Association for Computational Linguistics, Oct. 2014, pp. 1724–1734. doi: 10.3115/v1/D14-1179.

[2] K. Greff, R. K. Srivastava, J. Koutník, B. R. Steunebrink, and J. Schmidhuber, “LSTM: A Search Space Odyssey,” IEEE Transactions on Neural Networks and Learning Systems, vol. 28, no. 10, pp. 2222–2232, Oct. 2017, doi: 10.1109/TNNLS.2016.2582924.

[3] S. Hochreiter and J. Schmidhuber, “Long Short-Term Memory,” Neural Comput., vol. 9, no. 8, pp. 1735–1780, Nov. 1997, doi: 10.1162/neco.1997.9.8.1735.

[4] José C. García Alanis, “Introduction to EEG-preprocessing.” Accessed: Jun. 10, 2026. [Online]. Available: https://g0rella.github.io/gorella_mwn/preprocessing_eeg.html#now-that-the-raw-file-has-been-created-we-can-plot-the-data

[5] R. Leeb, F. Lee, C. Keinrath, R. Scherer, H. Bischof, and G. Pfurtscheller, “Brain–Computer Communication: Motivation, Aim, and Impact of Exploring a Virtual Apartment,” IEEE Transactions on Neural Systems and Rehabilitation Engineering, vol. 15, no. 4, pp. 473–482, Dec. 2007, doi: 10.1109/TNSRE.2007.906956.

[6] Michael C. Mozer, “A Focused Backpropagation Algorithm for Temporal Pattern Recognition | Request PDF.” Accessed: Jun. 10, 2026. [Online]. Available: https://www.researchgate.net/publication/243781476_A_Focused_Backpropagation_Algorithm_for_Temporal_Pattern_Recognition

[7] K. B. Petersen and M. S. Pedersen, [ http://matrixcookbook.com ], Version 20121115. Technical University of Denmark, 2012.

[8] G.-B. Zhou, J. Wu, C.-L. Zhang, and Z.-H. Zhou, “Minimal gated unit for recurrent neural networks,” Int. J. Autom. Comput., vol. 13, no. 3, pp. 226–234, Jun. 2016, doi: 10.1007/s11633-016-1006-2.