Multiway Multislice PHATE: Visualizing Hidden Dynamics of RNNs through Training

## Abstract

Abstract Recurrent neural networks (RNNs) are a widely used tool for sequential data analysis, however, they are still often seen as black boxes of computation. Understanding the functional principles of these networks is critical to developing ideal model architectures and optimization strategies. Previous studies typically only emphasize the network representation post-training, overlooking their evolution process throughout training. Here, we present Multiway Multislice PHATE (MM-PHATE), a novel method for visualizing the evolution of RNNs’ hidden states. MM-PHATE is a graph-based embedding using structured kernels across the multiple dimensions spanned by RNNs: time, training epoch, and units. We demonstrate on various datasets that MM-PHATE uniquely preserves hidden representation community structure among units and identifies information processing and compression phases during training. The embedding allows users to look under the hood of RNNs across training and provides an intuitive and comprehensive strategy to understanding the network’s internal dynamics and draw conclusions, e.g., on why and how one model outperforms another or how a specific architecture might impact an RNN’s learning ability.

## 1 Introduction

Recurrent neural networks (RNNs) are artificial neural networks (ANNs) designed for sequential data. Unlike feedforward neural networks (FNNs), which treat each input as independent, RNNs model sequences of inputs to produce one or a sequence of outputs. RNNs achieve this by processing each sequence element in order while retaining a memory of past inputs through memory units (e.g., LSTMs or GRUs) or through recurrent feedback connections [, ]. This memory effectively passes information forward in time, updating the internal state of the RNN with each input to reflect the sequence’s context. RNNs are thus particularly suited for tasks where the order and relationship between elements are crucial for understanding the whole sequence, e.g., neural time-series, postural action data, etc. [].
These intrinsic properties, as well as the developments of different architectures [], have led to RNNs showcasing exceptional performance across numerous domains, such as natural language processing (NLP), neuroscience and biomedical signal processing [, , ].

Despite their extensive use, RNN dynamics remain poorly understood. The opaque nature of their learned representations complicates the interpretation of their performance and robustness, hinders the selection of appropriate architectures and training parameters, and slows the development of more effective models. Significant efforts have been made to address the "black-box” nature of FNNs [, , ]. However, similar advances have not been as prevalent for RNNs, which have primarily been analyzed from the standpoint of dynamical systems analyses after training (e.g., via fixed points) [], or by comparing the performance of different RNN architectures at the network level and investigating the role of the various components within them [, ].
Thus, new methods that facilitate the interpretation of RNNs’ latent representations and their evolution during training are clearly needed.

In the field of explainable deep learning, dimensionality reduction methods are one of the more popular techniques [, , , , , ]. These approaches are useful as they allow human interpreters to readily visualize the structure of high-dimensional data, allowing them to gain intuition on structures and interactions within the data quickly. Nevertheless, they may inadvertently obscure meaningful data relationships by either emphasizing local (e.g., t-SNE []) or global (e.g., PCA [] or Isomap []) structures, or possess sensitivity to noise and outliers []. These challenges are compounded in cases such as RNNs where structures along multiple dimensions (time, epoch, unit) must all be preserved to produce summaries that capture the data’s complexity. As a result, traditional dimensionality reduction techniques fall short in visualizing learning in RNNs.

To study learning in FNNs, Gigante et al. [] introduced Multislice PHATE (M-PHATE) to visualize the network’s hidden state’s representations during training.
M-PHATE combines a multislice graph representation to model the community structure and temporal relationship between hidden units over epochs, with Potential of Heat Diffusion for Affinity-based Transition Embedding (PHATE) [] for dimensionality reduction.
Each slice in the multi-slice graph captures the network’s state at a specific epoch during training. This collection of graphs represents the dynamical system that governs the evolution of the network’s hidden states, where each hidden unit is connected to itself across epochs.

While M-PHATE requires only the activations of hidden units during training, the generated visualization captures key properties of network performance, such as test error accuracy, without any held-out validation data. This property is particularly advantageous in small-data scenarios where omitting data for testing can impact model performance.

Despite its numerous advantages, M-PHATE is designed to visualize the evolution of FNNs across epochs but fails to account for the recurrent and sequential nature of RNNs. Research has shown that hidden states from all time-steps are crucial for RNNs’ representations []. Therefore, to fully understand the evolution of these representations, it is essential to consider network dynamics across both time-steps and epochs. To address this need, we propose Multiway Multislice PHATE (MM-PHATE). Our method captures the latent dynamics of RNNs by visualizing hidden states across both time-steps and training epochs, providing deeper insights into the complex learning processes of RNNs. Our findings indicate that MM-PHATE retains considerably more dynamic details necessary for understanding RNNs’ performance compared to PCA, t-SNE, and Isomap.

Our main contributions are as follows:

- •
We present MM-PHATE, a novel framework for visualizing the hidden dynamics of RNNs across both time-steps and epochs simultaneously, providing a new perspective to RNN’s leaning trajectory, learned representation, and model performance.
- •
We demonstrate that MM-PHATE uniquely preserves the hidden representation community structure of hidden units throughout training by tracking each unit’s learning trajectory and the correlations among their activations.
- •
Applying MM-PHATE to RNN dynamics identifies phases of information processing and compression during learning, an observation that aligns with information bottleneck theory.

## 2 Related Work

We group existing methods to interpret RNNs into two categories: 1) performance-oriented and 2) application-oriented post-training analysis.

Performance-oriented analyses include investigating the role of components within various RNN architectures, as well as comparing different architectures and training parameters based on their network-level performance post-training. For example, Chung et al.[] conducted a performance evaluation comparing gated RNNs, particularly GRUs and LSTMs. Similarly, Greff et al.[] carried out a comprehensive analysis of LSTM components. While these studies interpret the performance of different RNNs at the network level, they did not delve deeply into analyzing their hidden states, thus providing limited insight into the nature of the representations learned by these networks.

In contrast, application-oriented analyses of trained networks often involve visualizing and interpreting activation maps that depict the representations learned by hidden units after training, often in the context of a specific task or network configuration. Numerous studies have applied this approach within the field of NLP. For instance, Karpathy et al. overlaid activation heat maps over texts, demonstrating that certain units developed interpretable representations, such as tracking text length, new line beginnings, and quote initiations []. Li et al. used saliency heat maps to identify words critical to the learned representations []. Strobelt et al. and Ming et al. created interactive visualizations to correlate hidden state patterns with phrases in texts [, ]. Beyond NLP applications, some studies explore other domains, including speech recognition [], earth sciences []), and medical states []. While these studies provide intuitive insights into the representations learned by the networks post-training, they are tailored for a particular task and may not generalize to other tasks.
Other studies tailored to understand general RNN properties include applying Proper Orthogonal Decomposition (POD) to the internal states of encoder and decoder units in Seq2Seq RNNs[] or employing PCA to visualize the activation of recurrent units and its link to generalization[].

While these application-oriented studies have identified critical aspects of RNNs’ hidden representation, they largely overlook the network’s learning trajectory over training epochs, a crucial process that must be understood to improve RNN-based machine learning systems. To the best of our knowledge, no existing visualizations of RNNs’ hidden states interpret the hidden dynamics simultaneously across time-steps and epochs.

## 3 Background

##### PHATE:

The PHATE algorithm[] optimizes the diffusion kernel [] for the visualization of high dimensional data. Let $\bm{x}_{i}$ be a point in a high-dimensional dataset. PHATE begins by computing the Euclidean distance matrix $\bm{D}$ between all data points, where $\bm{D}_{ij}=\|\bm{x}_{i}-\bm{x}_{j}\|_{2}$. These distances are then transformed into affinities using an adaptive $\alpha$-decay kernel
$\bm{K}_{k,\alpha}(\bm{x}_{i},\bm{x}_{j})=\frac{1}{2}\exp\left(-\left(\frac{\bm
{D}_{ij}}{\epsilon_{k}(\bm{x}_{i})}\right)^{\alpha}\right)+\frac{1}{2}\exp
\left(-\left(\frac{\bm{D}_{ij}}{\epsilon_{k}(\bm{x}_{j})}\right)^{\alpha}\right)$,
which adapts to the data density around each point and captures local information. The parameters $\epsilon_{k}(x_{i})$ and $\epsilon_{k}(x_{j})$ are the $k$-nearest-neighbor distance of $x_{i}$ and $x_{j}$, and $\alpha$ controls the decay rate. The affinities are then row-normalized to obtain the diffusion operator $\bm{P}=\bm{D}^{-1}\bm{K}_{k,\alpha}$ that represents the single-step transition probabilities between data points, where $\bm{D}$ is a diagonal matrix of row sums of $\bm{K}_{k,\alpha}$.
PHATE calculates the information distance between points based on their transition probabilities $\textrm{dist}_{ij}=\sqrt{\|\log{\bm{P}_{i}^{t}}-\log{\bm{P}_{j}^{t}}\|^{2}}$, where $\bm{P}^{t}$ captures the transition probabilities of a diffusion process on the data over $t$ steps and $i$ and $j$ are rows in the matrix.
These distances are embedded into low dimensions using Multidimensional Scaling (MDS) [] for visualization. For further details, see Moon et al. [].

##### M-PHATE:

Gigante et al. [] model the evolution of the hidden units and their community structure using a multislice graph. Let $F$ be an FNN with a total of $m$ hidden units, and let $F^{(\tau)}$ be the representation of the network after being trained for $\tau\in\{1,...,n\}$ epochs on the training data $\bm{X}$. The algorithm first calculates a shared feature space using the normalized activations of all hidden units $(1,...,m)$ on the input data, as a 3-tensor:

$$ $\bm{T}(\tau,i,k)=\frac{\bm{F}_{i}^{(\tau)}(\bm{Y}_{k})-\frac{1}{p}\sum_{\ell} \bm{F}_{i}^{(\tau)}(\bm{Y}_{\ell})}{\sqrt{\mathrm{Var}_{\ell}[\bm{F}_{i}^{( \tau)}(\bm{Y}_{\ell})]}},$ $$

where $\bm{F}_{i}^{(\tau)}(\bm{Y}_{k}):\mathbb{R}^{d}\rightarrow\mathbb{R}$ denote the activation of the $i$-th hidden unit of $\bm{F}$ for the $k$-th sample from input data $\bm{Y}$. Here, $\bm{Y}$ is a subset of $p$ samples from the $d$-dimensional training data $\bm{X}$, with an equal number of samples from each input class. The training data $\bm{X}$ is drawn from a larger dataset $\bm{\Pi}$, and $p\ll|\bm{X}|$. This activation tensor $\bm{T}$ is then used to calculate intraslice affinities between pairs of hidden units within an epoch $\tau$ during the training, as well as the interslice affinities between a hidden unit $i$ and itself at different epochs:

$$ $\bm{K}_{\text{intraslice}}^{(\tau)}(i,j)=\exp\left(\frac{-\|\bm{T}(\tau,i)-\bm {T}(\tau,j)\|_{2}^{\alpha}}{\sigma^{\alpha}_{(\tau,i)}}\right),\bm{K}_{\text{ interslice}}^{(i)}(\tau,\upsilon)=\exp\left(\frac{-\|\bm{T}(\tau,i)-\bm{T}( \upsilon,i)\|_{2}^{2}}{\epsilon^{2}}\right)$ $$

where $\alpha$ is the $\alpha$-decay parameter, $\sigma_{(\tau,i)}$ is the intraslice bandwidth for unit $i$ in epoch $\tau$, and $\epsilon_{n}$ is the fixed interslice bandwidth.
These matrices are then combined to form an $nm\times nm$ multislice kernel matrix $\bm{K}$, which is then symmetrized, row-normalized, and visualized using PHATE in 2D or 3D.

## 4 Multiway Multislice PHATE

M-PHATE was shown to be a powerful tool for visualizing FNNs. However, to effectively visualize the evolution of RNNs’ hidden representations, we need to consider hidden state dynamics across *time-steps* within the sequence and training epochs concurrently. In RNNs, the output from previous *time-steps* is fed as an input to current *time-steps*. This is useful in the treatment of sequences and building a memory of the previous inputs into the network. The network iteratively updates a hidden state $h$. At each *time-step* $t$, the next hidden state $h_{t+1}$ is computed using the input $x_{t}$ and the current hidden state $h_{t}$. Importantly, the network uses the same weights $W$ and biases $b$ for each *time-step*. Thus the output $y_{t}$ at *time-step* $t$ is $y_{t}=f(W\cdot h_{t}+b)$, where $f$ is some activation function.

Let $\bm{R}^{(\tau)}$ be the representation of an $m$-unit RNN after being trained for $\tau\in\{1,\ldots,n\}$ epochs on the training data $\bm{X}$ split from dataset $\bm{\Pi}$. We denote $R_{i,w}^{(\tau)}(\bm{Y}_{k}):\mathbb{R}^{d}\rightarrow\mathbb{R}$ the activation of the $i$-th hidden unit of $R$ at time-step $w\in\{1,\ldots,s\}$ in epoch $\tau$ for the $k$-th sample of $\bm{Y}$, where $\bm{Y}$ consists of $p$ samples from the training data $\bm{X}$.
We construct the 4-way tensor $\bm{T}$ using the hidden unit activations as a shared feature space, which we use to calculate unit affinities across all epochs and time-steps.
The tensor $\bm{T}$ is an $n\times s\times m\times p$ tensor containing the activations at each epoch $\tau\in\{1\dots n\}$ and time-step $w\in\{1\dots s\}$ of each hidden unit $\bm{R}_{i}$ ($i\in\{1\dots m$}) with respect to each sample $\bm{Y}_{k}\subset\bm{X}$. To eliminate the variability in $\bm{T}$ due to the bias term $b$, we $z$-score the activation of each hidden unit at time-step $w$ and epoch $\tau$:

$$ $\bm{T}(\tau,\omega,i,k)=\frac{\bm{R}_{i,\omega}^{(\tau)}(\bm{Y}_{k})-\frac{1}{ p}\sum_{\ell}\bm{R}_{i,\omega}^{(\tau)}(\bm{Y}_{\ell})}{\sqrt{Var_{\ell}[\bm{R }_{i,\omega}^{(\tau)}(\bm{Y}_{\ell})]}}.$ (1) $$

We construct a kernel over $\bm{T}$ utilizing our prior knowledge of the temporal aspect of $\bm{T}$ to capture its dynamics over epochs and time-steps. This constructed kernel, denoted $\bm{K}$, represents the weighted edges in the multislice graph of the hidden units. In this representation, each unit has two types of connections: edges between the unit to itself across epochs and time steps and, within a fixed epoch and time-step, edges between a unit and its community—the other units which have the most similar representation. The edges are weighted by the similarity in activation pattern.
We define $\bm{K}$ as a $nsm\times nsm$ kernel matrix between all $m$ hidden units at all $s$ time-steps in all $n$ training epochs. The $((\tau-1)sm+(\omega-1)m+j)_{th}$ row or column of $\bm{K}$ refers to the $j_{th}$ unit at time-step $w$ in epoch $\tau$. We henceforth refer to the row as $\bm{K}((\tau,\omega,j),:)$ and the column as $\bm{K}(:,(\tau,\omega,j))$.
In order to capture the evolution of hidden units of $R$ across time-steps and epochs, while preserving the unit’s community structure, we construct a multiway multislice kernel matrix reflecting two types of connections simultaneously:

- •
Intrastep affinities between hidden units $i$ and $j$ at time-step $\omega$ in epoch $\tau$:
$\bm{K}_{\text{intrastep}}^{(\tau,\omega)}(i,j)=\exp(-\parallel\bm{T}(\tau,
\omega,i)-\bm{T}(\tau,\omega,j)\parallel^{\alpha}_{2}/\sigma^{\alpha}_{(\tau,
\omega,i)})$
- •
Interstep affinities between a hidden unit $i$ and itself at different time-steps and epochs:
$\bm{K}_{\text{interstep}}^{(i)}((\tau,\omega),(\eta,\nu))=\exp(-\parallel\bm{T
}(\tau,\omega,i)-\bm{T}(\eta,\nu,i)\parallel^{2}_{2}/\epsilon^{2})$

where $\alpha$ is the $\alpha$-decay parameter, $\sigma_{(\tau,\omega,i)}$ is the intrastep bandwidth for unit $i$ at time-step $w$ and epoch $\tau$, and $\epsilon$ is the fixed interstep bandwidth.

The bandwidth $\sigma_{(\tau,\omega,i)}$ of the $\alpha$-decay kernel is set to be the distance of unit $i$ at time-step $w$ from epoch $n$ to its $k$-th nearest neighbor across units at that time-step and epoch: $\sigma_{(\tau,\omega,i)}=d_{k}(\bm{T}(\tau,\omega,i),\bm{T}(\tau,\omega,:))$, where $d_{k}(z,Z)$ denotes the $\ell_{2}$ distance from $z$ to its $k$-th nearest neighbor in $Z$. We used $k=5$ in all the results presented. The use of this adaptive bandwidth means the kernel is not symmetric and thus requires symmetrization. In the interstep affinities $\bm{K}^{(i)}_{\textrm{interstep}}$, we use a fixed-bandwidth Gaussian kernel $\epsilon=\frac{1}{nsm}\sum_{\tau=1}^{n}\sum_{\omega=1}^{s}\sum_{i=1}^{m}d_{k}(
\bm{T}(\tau,\omega,i),\bm{T}(:,:,i))$, the average across all time-steps in all epochs and all units of the distance of unit $i$ at time-step $t$ to its $k_{th}$ nearest neighbor among the set consisting of the same unit $i$ at all steps.

The combined kernel matrix of these two matrices contains one row and column for each unit at each time-step in each epoch, such that the intrastep affinities form a block diagonal matrix and the interstep affinities form off-diagonal blocks composed of diagonal matrices (Fig. 1a, 1b).

$$ $\bm{K}((\tau,\omega,i),(\eta,\nu,j))=\begin{cases}\bm{K}^{(\tau,\omega)}_{ \text{intrastep}}(i,j)&\text{if }(\tau,\omega)=(\eta,\nu)\\ \bm{K}^{(i)}_{\text{interstep}}((\tau,\omega),(\eta,\nu))&\text{if }i=j\\ 0&\text{otherwise}\end{cases}$ $$

We symmetrize this final kernel as $\bm{K}^{\prime}=\frac{1}{2}(\bm{K}+\bm{K}^{T})$, and row normalize it to obtain $\bm{P}=\bm{D}^{-1}\bm{K}^{\prime}$, which represents a random walk over all units at all time-steps in all epochs, where propagating from one state to another is conditional on the transition probabilities between time-step $\omega$ in epoch $\tau$ and time-step $\nu$ in epoch $\eta$. PHATE is applied to $\bm{P}$ to visualize the tensor $\bm{T}$ in two or three dimensions. This resulting visualization thus simultaneously captures information regarding the evolution of the units across both time-steps and epochs.

Figure: (a)
Refer to caption: extracted/5641925/figures/multiway_multislice_graph.png

## 5 Results

We demonstrate the ability of MM-PHATE to capture useful properties of RNN learning on two datasets: 1) The Area2Bump dataset [] consisting of neural activity recorded from Brodmann’s area 2 of the somatosensory cortex while macaque monkeys performed a slightly modified version of a standard center-out reaching task and 2) The Human Activity Recognition (HAR) Using Smartphones dataset [], kinematic recordings of 30 subjects performing daily living activities with a smartphone embedded with an inertial measurement unit.

### 5.1 Neural activity

We begin with the Area2Bump dataset, which consists of spiking activity data from macaques []. We trained an LSTM network comprised of a single layer of 20 units to classify the direction of arm-reaching movements, and applied MM-PHATE to visualize the learning of the LSTM.
In an example session where our network achieved a validation accuracy of 74%, we applied the MM-PHATE visualization to the tensor $\bm{T}$ that consisted of the network activations of all 20 units over 600 time-steps for each of 200 training epochs. In practice, we sampled time-steps and epochs to reduce memory load. Each point in the visualization represents a hidden unit at a given time-step in a given epoch (Fig. 2). Through this visualization, we observed a smooth transition of the hidden states across both time-steps and epochs, reflecting the dynamic changes in the network’s internal representations throughout training.
We compare MM-PHATE to three other dimensionality reduction techniques: PCA, t-SNE, and Isomap, using the same T tensor. We first flattened $\bm{T}$ along the epoch, time-step, and unit axis, and embedded it with each of the dimensionality reduction methods.
Notably, the MM-PHATE visualization reveals a distinct split in representations during the later time-steps and epochs, highlighting unique learning patterns as the model converges—patterns that are not discernible with other methods. PCA and Isomap, while showing a seemingly smooth transition, do not show distinct differences between the early and late epochs, failing to capture how the representation transforms during learning. On the other hand, t-SNE resulted in a visualization that lacked smooth transitions across time-steps and epochs, possibly due to its sensitivity to noise.

Figure: Figure 2: Area2Bump: Visualization of a 20-unit LSTM network trained for 200 epochs. The visualizations are generated using MM-PHATE, PCA, t-SNE, and Isomap, from left to right, respectively. Points are colored based on epoch (top row) or time-step (bottom-row).
Refer to caption: extracted/5641925/figures/combined_plots.png

#### 5.1.1 Intra-step entropy

To evaluate the effectiveness of our embedding in capturing meaningful structures within the network’s hidden dynamics, we next analyzed the flow of information during training (Fig. 3). Specifically, we aimed to analyze the spread of points in the MM-PHATE space, which corresponds to the diversity of the internal representations of the network. To quantify this property, we computed the entropy between hidden units in the embedded space at each time-step $\omega$ and compared these intra-step entropies with the accuracy and loss metrics recorded at the end of each training epoch. Conceptually, we model a general RNN trained on dataset $\bm{X}$ with label $\bm{L}$ as a Markov chain $(\bm{L}\rightarrow\bm{X}\rightarrow\bm{R})$. This allows us to compute the mutual information between $\bm{X}$ and $\bm{R}$ as $I(\bm{X},\bm{R})=H(\bm{R})-H(\bm{R}|\bm{X})$, where $H(\bm{R})$ and $H(\bm{R}|\bm{X})$ are the marginal and conditional entropies. Given the deterministic nature of RNNs, $H(\bm{R}|\bm{X})$ equates to zero, indicating that $H(\bm{R})$ at each time-step $\omega$ reflects the input information the network retains during training [, ].

Figure: Figure 3: Area2Bump: Intra-step entropy of all hidden units in embedding space at each time-step in each epoch, compared to training and validation accuracy (top) and losses (bottom), comparing embeddings of MM-PHATE, PCA, t-SNE, and Isomap.
Refer to caption: extracted/5641925/figures/combined_plots_intrastep_entropy.png

Our analysis of the MM-PHATE embedding revealed a general increase in the entropy (Fig. 3) that aligns with the training epoch where the network begins to overfit (approximately epoch 100). We hypothesize that this increase may be due to memorization of noise or other nuisance variance in the training data. The observed changes in entropy at specific time-steps and epochs coincided with shifts in model performance throughout the training. This suggests that our algorithm successfully captures and retains dynamical information critical for the model’s learning process. Notably, transitions in validation accuracy are not always reflected in changes in the training loss curve, e.g., fluctuation of validation accuracy curve before epoch 30 cannot be observed in the training loss (Fig. 3). The MM-PHATE embedding, however, successfully identified relevant transitions coinciding with shifts in both performance measures, thus offering insights beyond those provided by accuracy and loss functions alone. Moreover, the MM-PHATE embedding reveals details on information processing through time-steps. This is not the case for the loss and accuracy metrics, which are measured once per epoch. With MM-PHATE, we obtain intra-step entropy plots for each time-step, and can thus compared dynamics across them. Notably, we observe diverging patterns between the entropy of earlier and later time steps. In this network, earlier time steps’ entropy increases and plateaus around epoch 100, while later time steps mimic a spiking activity around that epoch, eventually returning back to the entropy value they had in earlier epochs. While the effects of these diverging patterns remains uncertain, MM-PHATE clearly offers additional insights into the model’s dynamics.

We performed the same intra-step entropy analysis on other methods (Fig. 3). PCA failed to capture the critical dynamics associated with early subtle changes in model performance, likely due to its linear nature and emphasis on the data’s global structures. Similarly, Isomap failed to capture significant dynamics related to performance changes before the onset of overfitting. Although both PCA and Isomap exhibited a rise in entropy, this increase occurred after the loss had plateaued. In contrast, MM-PHATE demonstrated an earlier rise in entropy, aligning with the transition into the plateau and thus more accurately capturing the underlying dynamics. On the other hand, t-SNE struggled to identify the increase in retained information corresponding to model overfitting and was generally inadequate in capturing the global structural dynamics evident in other techniques.

#### 5.1.2 Inter-step entropy

To further assess the quality of our embedding, we quantified the inter-step entropy of hidden unit activations across various time-steps $\omega$ (Fig. 4a). This metric measures the entropy between a unit’s activations at different time-steps, providing insights into the temporal dynamics of each hidden unit. Our analysis identified distinct patterns among the hidden units. Notably, certain units exhibited significantly higher inter-step entropy, indicative of higher sensitivity to input changes over time and their potential role in capturing complex dependencies. These units peaked around epoch 80, aligning with the peak in intra-step entropy at later time-steps. This suggests that these units are largely responsible for the increase in mutual information at these time-steps.

Figure: Figure 4: Area2Bump: a) Inter-step entropy of each hidden unit across epochs, plotted alongside model accuracy (left) and loss (right). b) Clusters of hidden units from the model’s inter-step entropy trajectories across epochs. Left: Trajectories of the units in cluster 0 (12 units) and cluster 1 (8 units). Right: cluster center trajectories across epochs. c) Confusion matrices of clusters on training and d) validation data.
Refer to caption: extracted/5641925/figures/interstep_clustering_combined.png

To validate the insights from our entropy analyses and confirm the significance of the units’ distinct learning behaviors with respect to the learned representation and model performance, we clustered the hidden units into 2 groups, expecting one group to showcase the higher inter-step entropy. We then analyzed each group’s specific properties. For clustering, we used k-means clustering with the Dynamic Time Warping (DTW) metric to group the hidden units into clusters based on their inter-step entropy trajectories (Fig. 4b). DTW determines an optimal match between two time series, making it more suitable as a distance metric for our sequence data. Clustering revealed distinct trends in inter-step entropy trajectories across epochs. Cluster 1 (8 units) corresponds to the units previously identified as exhibiting higher inter-step entropy, and units in cluster 0 (12 units) showed a sharp decrease in inter-step entropy around epoch 100.
Following the clustering, we built new networks comprised of the post-training weighted units contained in one particular cluster, forming two sub-networks of the original network architecture to analyze their learning capabilities. By computing confusion matrices for each sub-network and comparing their predictive performance on both training and validation data, we assessed the quality of their learned representation (Fig. 4c-d). Interestingly, despite having fewer units, cluster 1 performed significantly better with both training and validation data. This outcome aligns with our previous prediction, where cluster 1 was identified as learning more complex dependencies from the input. This differential learning capability between clusters demonstrates the utility of our visualization and clustering approach in revealing critical differences in how information is processed and represented within the network.

Previous studies have indicated that the clustering property of hidden unit activation is crucial for understanding the quality of the learned representation [, , ]. Our findings affirm the importance of capturing each unit’s temporal dynamics, their community structure, as well as the evolution of these structures across time-steps and epochs.

Comparisons with other dimensionality reduction techniques highlighted their limitations (Fig. S1). PCA and Isomap failed to capture subtle dynamic variations, particularly in early epochs before the onset of overfitting. t-SNE displayed a general trend corresponding to each performance transition, however, the trajectories are noisy. More importantly, all these methods failed to differentiate individual units’ learning behavior or to capture their community structure effectively. This emphasizes the superior performance of MM-PHATE in maintaining the integrity of the hidden dynamics.

### 5.2 Analysis with Human Activity Recognition Model

We next considered an action recognition dataset.
We trained 30-unit LSTM network designed for kinematics-based Human Activity Recognition (HAR) []. The model was trained for 1000 epochs, achieving a final validation accuracy of 84% (Fig. 5).
We applied MM-PHATE to the tensor of hidden unit activations and repeated our analysis as in the Area2Bump dataset.
Similar to the patterns observed in the Area2Bump LSTM network, the HAR network displayed an increase in intra-step entropy with the onset of overfitting. Particularly noteworthy was the gradual increase in entropy across time-steps before epoch 300. However, after this epoch, entropy at later time-steps significantly dropped, while entropy at early time-steps remained stable throughout the training epochs. This indicates a reduction in mutual information between the input and the hidden states as the model processed more input over time. This behavior aligns with findings by Farrell et al.[], who reported similar dynamics of information expansion and compression across time-steps in trained RNNs. Their results further elucidated how gradient-based learning mechanisms contribute to the development of robust representations by balancing these processes of expansion and compression.

Further insights were gained from the inter-step entropy analysis, which indicated that the entropy of many hidden units began to rise significantly around epoch 300, coinciding with a decrease in intra-step entropy and an improvement in model performance. As previously discussed, high inter-step entropy indicates that the units are more sensitive to changes in input across time and may be crucial for the model to learn complex dependencies. This observation also implies that despite the increase in time-dependent information retained by many hidden units, the network successfully compresses this information in its overall learned representation across epochs. Such compression seems beneficial to model performance and aligns with the principles in information bottleneck theory. This theory states that deep network learning involves a fitting phase followed by a compression phase, during which useful information is distilled from the input to enhance generalization [, , ].

Building on the discussion of information compression in the HAR model, it is pertinent to question why the Area2Bump model does not exhibit a similar compression phase. Cheng et al. [] suggest that models with insufficient generalizability often fail to demonstrate a compression phase in practice, especially when applied to complex datasets with relatively simple network architectures. Given that our models have comparable structures, we investigated the complexity of the two datasets by analyzing the number of principal components required to explain 95% of the variance. We discovered that the Area2Bump data is significantly more complex than the HAR data, necessitating 35 principal components compared to only 6 for the HAR data (Fig. 6). This difference in complexity likely accounts for the absence of a noticeable compression phase in the Area2Bump model.

Figure: Figure 5: HAR: a) MM-PHATE visualization of a 30-unit LSTM network trained on HAR data, colored by epoch (left) and time-step (right). b) Intra-step entropy, c) and inter-step entropy.
Refer to caption: extracted/5641925/figures/kinematics/har.png

Our analysis demonstrates that information compression occurs across both time-steps and epochs in RNNs, closely aligning with performance improvements. These results affirm the practical utility of the information bottleneck theory in RNNs and confirm that MM-PHATE effectively reveals detailed insights into the model’s hidden representation and its evolution. Further research should investigate the implications and interactions of various compression dynamics in RNNs, which could lead to more robust and generalizable network architectures.

Figure: Figure 6: PCA on the 2 datasets. The Area2Bump data requires 35 PCs (left) to cover 95% of the variance. The HAR data (right) requires 6 PCs to cover 95% of variance.
Refer to caption: extracted/5641925/figures/PCA_2data.png

## 6 Limitations

While we demonstrate here MM-PHATE’s ability to provide a detailed qualitative descriptions an RNN’s internal representation during training, our approach rests on a number of assumptions. First we assume continuity in time and over epochs (as encoded in the structured graph kernel). Using large learning rates, having multiple significant restarts, or using discontinuous activation functions, while not common in the literature, can violate this assumption and provide less informative visualizations.
Moreover there are many domains in which RNNs are applied, such as LLMs etc. A number of computational stages, such as the graph creation, can be slow for extremely large systems, requiring additional efficiency advances to analyze and compare models with larger architectures or training data []. Finally, MM-PHATE does not use any infomation of the internal structure of the RNN, such as attention connections in LLMs, which are another direction of future research.

## 7 Conclusion

Here we introduced MM-PHATE, a novel dimensionality reduction technique designed to visualize the hidden dynamics of RNNs throughout training. MM-PHATE uniquely captures the evolution of these dynamics over both time-steps and epochs, surpassing the capabilities of existing visualization tools.
Our intra- and inter-step entropy analysis demonstrates that MM-PHATE can unveil learning behaviors obscured by traditional metrics such as accuracy and loss curves. The ability to directly observe learned representations is particularly valuable when data is limited, as it does not rely on external validation data. Furthermore, MM-PHATE can distinguish between unique learning behaviors among subgroups of hidden units, elucidating their roles in modulating information flow within the network. This analysis aligns closely with the principles of the information bottleneck theory, offering a robust framework to dissect the network’s training into distinct phases of information expansion and compression.

While this study marks an important step in validating the effectiveness of MM-PHATE, future studies can work on reducing memory consumption and delve more deeply into various RNN architectures and further analyze the unique learning behaviors of individual hidden units. Such explorations could greatly enrich our understanding of RNN dynamics and pave the way for more robust and effective network designs. Finally, our work could be extended to transformers, which rival RNNs in sequential data treatment, and whose hidden dynamics of learned representation have yet to be elucidated.

## Supplemental Materials

### 7.1 Codes

Codes will be shared upon publication.

Figure: Figure S1: Area2Bump: Inter-step entropy of all hidden units in embedding space of the Area2Bump model at each time-step in each epoch, compared to training and accuracies (top) and losses (bottom). From left to right, the dimensionality reduction metrics used are MM-PHATE, PCA, t-SNE, and Isomap.
Refer to caption: extracted/5641925/figures/combined_plots_interstep_entropy.png

### 7.2 Datasets

We employed two datasets.

- 1.
The Area2Bump dataset [] consists of neural activity recorded from Brodmann’s area 2 of the somatosensory cortex while macaque monkeys performed a slightly modified version of a standard center-out reaching task. Dataset license: CC0 1.0. According to the authors, data was collected consistently "with the guide for the care and use of laboratory animals and approved by the institutional animal care and use committee of Northwestern University under protocol #IS00000367". This dataset, collected from monkeys thus does not contain personally identifiable information. Since it is neural data, we do not consider it to be offensive content.
- 2.
The Human Activity Recognition (HAR) Using Smartphones dataset [], kinematic recordings of 30 subjects performing daily living activities with a smartphone embedded with an inertial measurement unit. Dataset license: CC BY 4.0. Information relating to participant consent was not found relating to this dataset. This dataset does not reveal participant name or identifiable information. The kinematics contained in the dataset are not considered offensive content.

### 7.3 Model Training

We used TensorFlow’s Keras API for all models training and validation.

The network in Section 5.1 was trained as follows. The Area2Bump dataset was randomly split into training and validation subsets containing an equal number of samples for each input class with an 8 to 2 ratio. Additional samples that would make the training data uneven were added back to the validation subset to make use of all samples. The network consists of an LSTM layer with 20 units. This is followed by a Flatten layer that converts the LSTM’s output into a one-dimensional vector. Finally, a Dense layer with 8 units and a softmax activation function produces the output for the 8-class classification tasks. The network was trained with a batch size of 64. We used an Adam optimizer with a learning rate of $1e^{-4}.$ During the training process, we recorded the activations from the LSTM layer into the activation tensor $\bm{T}$ for visualization.

The network in Section 5.2 was trained as follows. The HAR dataset was preprocessed and split by the authors into training and validation subsets according to the subjects with a 7 to 3 ratio. The network consists of an LSTM layer with 30 units and a Dense layer with 6 units and a softmax activation function produces the output for the 6-class classification tasks. The network was trained with a batch size of 32. We used an Adam optimizer with a learning rate of $2e^{-5}.$ During the training process, we recorded the activations from the LSTM layer into the activation tensor $\bm{T}$ for visualization.

### 7.4 MM-PHATE

Due to memory constraints, we only used a subset of the tensor $\bm{T}$ for MM-PHATE computation. Specifically, for the Area2Bump model in section 5.1, epochs were sampled using an array combining the first 29 epochs with every 5th epoch thereafter to cover the initial rapid learning phase, and intrinsic steps were sampled using a linear space from 0 to the end (600), resulting in 100 evenly spaced steps. For the HAR model in section 5.2, we sampled the epochs using an array combining the first 29 epochs with every 10th epoch thereafter. We utilized the M-PHATE package to construct our multiway multislice graphs and for the application of PHATE.

### 7.5 PCA

PCA was performed using the "PCA" class from the "sklearn.decomposition" package to reduce the dimensionality of the time trace tensor $\bm{T}$—recorded during training of the Area2Bump model—to three principal components.

### 7.6 t-SNE

In this analysis, we first conducted an initial dimensionality reduction on the same time trace tensor $\bm{T}$ with PCA to 15 principle components, which explains 99.93% of the variance in the activations. Subsequently, t-SNE with the Barnes-Hut approximation was performed using the "sklearn.manifold.TSNE" class.

### 7.7 Isomap

Due to memory constraints, we only used a subset of the tensor $\bm{T}$ for Isomap computation. Specifically, epochs were sampled using an array combining the first 29 epochs with every 10th epoch thereafter, and intrinsic steps were sampled using a linear space from 0 to the end (600), resulting in 50 evenly spaced steps. We first conducted an initial dimensionality reduction on the sampled tensor with PCA to 15 principle components using "sklearn.decomposition.PCA" package. Then, we applied Isomap using "sklearn.manifold.Isomap" class to reduce the dimensionality to 3.

### 7.8 Area2Bump with GRU

Here is the same analysis as section 5.1 using GRU (Fig. S2, S3, S4). Other parameters were kept the same.

From these figures, it is evident that PCA and t-SNE present similar visualizations of the hidden dynamics across different network architectures, while MM-PHATE distinctly captures the unique learning behaviors of each model. Consistent with Section 5.1, PCA displays an increasing intra-step entropy even after model accuracy has plateaued, and t-SNE produces a noisy visualization. In contrast, MM-PHATE uniquely aligns its transitions well with the learning curve. Notably, the GRU model’s representation appears more compact and organized compared to the LSTM model, potentially reflecting its superior performance and reduced overfitting.

Figure: Figure S2: Area2Bump GRU: Visualization of a 20-unit GRU network trained for 200 epochs. Each point represents a hidden unit at a specific time-step in a given epoch throughout the entire training process. The visualizations are generated using MM-PHATE, PCA, and t-SNE, from left to right, respectively. Points are colored based on epoch (top row) or time-step (bottom-row).
Refer to caption: extracted/5641925/figures/gru/combined_plots.png

Figure: Figure S3: Area2Bump GRU: Intra-step entropy of all hidden units in embedding space at each time-step in each epoch, compared to training and validation accuracy (top) and losses (bottom), comparing embeddings of MM-PHATE, PCA, and t-SNE.
Refer to caption: extracted/5641925/figures/gru/combined_plots_intrastep_entropy.png

Figure: Figure S4: Area2Bump GRU: Inter-step entropy of all hidden units in embedding space of the Area2Bump model at each time-step in each epoch, compared to training and accuracies (top) and losses (bottom). From left to right, the dimensionality reduction metrics used are MM-PHATE, PCA, and t-SNE.
Refer to caption: extracted/5641925/figures/gru/combined_plots_interstep_entropy.png

### 7.9 Area2Bump with Vanilla RNN

Here is the same analysis as section 5.1 using vanilla RNN (Fig. S5, S6, S7). Other parameters were kept the same.

From these figures, it is evident that regardless of the network architectures, PCA exhibits a revolving pattern with overly smooth transitions across epochs, while t-SNE produces a noisy visualization. In contrast, the MM-PHATE visualization reveals that the Vanilla RNN displays a more chaotic pattern compared to the LSTM and GRU models, which is likely associated with its reduced performance and increased overfitting. Furthermore, the intra-step entropies of MM-PHATE show reduced variation across time-steps, indicating that the model struggles to process the input data effectively to generate meaningful representations.

Figure: Figure S5: Area2Bump Vanilla: Visualization of a 20-unit Vanilla RNN trained for 200 epochs. Each point represents a hidden unit at a specific time-step in a given epoch throughout the entire training process. The visualizations are generated using MM-PHATE, PCA, and t-SNE, from left to right, respectively. Points are colored based on epoch (top row) or time-step (bottom-row).
Refer to caption: extracted/5641925/figures/SimpleRNN/combined_plots.png

Figure: Figure S6: Area2Bump Vanilla: Intra-step entropy of all hidden units in embedding space at each time-step in each epoch, compared to training and validation accuracy (top) and losses (bottom), comparing embeddings of MM-PHATE, PCA, and t-SNE.
Refer to caption: extracted/5641925/figures/SimpleRNN/combined_plots_intrastep_entropy.png

Figure: Figure S7: Area2Bump Vanilla: Inter-step entropy of all hidden units in embedding space of the Area2Bump model at each time-step in each epoch, compared to training and accuracies (top) and losses (bottom). From left to right, the dimensionality reduction metrics used are MM-PHATE, PCA, and t-SNE.
Refer to caption: extracted/5641925/figures/SimpleRNN/combined_plots_interstep_entropy.png

### 7.10 Computing Infrastructure

All but the t-SNE computation was done on a 14-core laptop running Windows 11 Home with a NVIDIA GeForce RTX 3070 Ti Laptop graphics card and 40GB of RAM. The t-SNE was done on a single 95-core internal cluster running Ubuntu 18.04.6 LTS with 10 Quadro RTX 5000 graphics cards and 755GB of RAM.

## Peer Reviews

confidence: 1
contribution: 3
initial_rating: 6
presentation: 3
questions: see weakness
soundness: 3
strengths: ## Strength
1. The visualization results look cool.

2. The motivation is clear and the understanding of RNN is meaningful.

3. The analyses are comprehensive.
summary: ## Summary
This paper aims to visualize the internal dynamics of RNNs. The authors claim that the previous papers merely achieve it via post-training, overlooking their evolution process. To this end, they proposed MM-PHATE, which is a graph-based embedding using structured kernels across the multiple dimensions spanned by RNNs: time, training epoch, and units. Experiments show that the proposed method helps users to understand RNNs.
weaknesses: ## Weakness
1. The codes and datasets are missing, limiting the reproducibility. 

2. Recommend generalizing the proposed method to more recent sequential models like transformers. 

3. The statistical information of the datasets are missing.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: Please see weaknesses.

- Can the authors please expand upon the utility of the visualizations and entropy based plots? I found them visually interesting but I could not see myself utilizing them when training or comparing amongst models? 

- Can the authors clarify a bit how different length sequences are represented in T. My understanding is that T is n×s×m×p, so is there padding for the different length s, or is this a non-issue for the formulation of K? 

- The two datasets evaluated on appear to be classification datasets. Is the author's formulation useful for other tasks? Please describe what insights could be derived from the MMPhate visualizations there. 

- I noticed that the authors did not compare to [1], which has been used in the past to analyze the dynamics of neural networks. I bring it up because it may be a stronger baseline to the TSNE/PCA. 

[1] SVCCA: Singular Vector Canonical Correlation Analysis for Deep Learning Dynamics and Interpretability
soundness: 2
strengths: - *Interesting Problem.* Understanding the learning dynamics of RNN representations is an interesting problem. It will enable additional insights into improving their architectures. Though transformers are certainly more popular these days, I think that RNNs are still useful in some areas, so this is a useful task.

- *Some Nice Insights.* The author's demonstration of the relationship between the entropy of the representation and loss/over-fitting and related experiments are interesting observations, especially considering existing hypotheses for what happens when learning RNN representations.
summary: This paper proposes a technique to visualize the hidden states of RNNs while training and across units. They extend PHATE and M-PHATE to account for within sequence time-steps, in addition to epoch-level time steps. They compute kernels over a 4-way tensor that creates two types of edges: (1) edges that look at the same unit evolving over the course of a sequence and epochs, and (2) edges that look at how different units are related to each other for a fixed timestep/epoch. The intrastep edges form block-diagonal matrices and the interstep edges create off diagonal diagonal matrices. The utility of the proposed representation is evaluated on two datasets (Area2Bump, Human Activity Recognition dataset). They demonstrate that looking at the inter and intra-step entropy provides insights into the loss and when models may begin to overfit. They also show that visualizations obtained are better than TSNE/PCA/etc.

UPDATE: I thank the authors for their feedback. I have read the other reviews and responses as well. While some of my concerns were addressed, I do not feel this paper is quite read for acceptance and thus will maintain my score. I encourage the authors to take into account some of the suggestions made by reviewers for expanding their analyses.
weaknesses: -  *Novelty.*  This paper builds off of PHATE and MMPhate to introduce a visualization model specifically for RNNs, which accounts for the sequential nature of the data. This involves some changes to the tensor and kernel to include additional edges amongst the sequence steps, in addition to the epoch steps. I do not see this as significant novelty.

- *Weak baselines/Limited Datasets.* The authors limit their analysis to one tasks (classification), on two datasets. The visualizations compared are to tsne/pca/isomap, which are known to be insufficient for neural network representations. I would have liked to seen more recently comparisons as well. I feel that the authors are missing some citations and recent work. 

- *Utility/Practicality*. While the visualizations are interesting, I am not persuaded that they are particularly useful or would generalize across different datasets, sizes, or variations in architecture. Nor am I convinced the additional computational time/memory to compute these visualizations would be useful if I did.

---

confidence: 2
contribution: 2
initial_rating: 3
presentation: 2
questions: What is the relationship between the visualization of MM-PHATE and the visualization of M-PHATE at each time step?
soundness: 2
strengths: - This paper is clearly written, well organized. 
- Experiments can effectively reflect the intended objectives of the model.
summary: This work introduces Multiway Multislice PHATE (MM-PHATE), an innovative visualization method designed to explore the internal dynamics of recurrent neural networks (RNNs). This method employs a graph-based embedding with structured kernels to capture the multiple dimensions of RNNs, including time, training epochs, and units. Through experiments on various datasets, the authors demonstrate that MM-PHATE effectively preserves the community structure of hidden representations and highlights distinct phases of information processing and compression during training.
weaknesses: - The technical contribution is incremental. The authors merely add the time step dimension to the M-PHATE and do not make any additional adaptation for RNN. 
- MM-PHATE is closely related to M-PHATE, but the experiment lacks analysis of M-PHATE.

---

confidence: 3
contribution: 1
initial_rating: 5
presentation: 2
questions: 1. The authors claim that MM-PHATE identifies information processing and compression phases. However, this appears to be based only on entropy change observations. Could you provide mathematical/theoretical proof of how this relates to information bottleneck theory?

2. How do you differentiate between useful information compression and mere information loss in your entropy analysis?

3. Regarding the claimed correlation between MM-PHATE visualization and model performance, what specific visualization characteristics quantitatively correlate with which performance metrics?

4. Has the correlation between visualization patterns and model performance been validated across different architectures and datasets? If so, how consistent are these correlations?

5. The paper claims to preserve community structure well. How do you define and measure "well-preserved community structure"? What metrics validate this preservation?

6. The authors emphasize the analysis of RNN temporal dynamics as a key objective yet use sampling (e.g., reducing 600-time steps to 100 in the Area2Bump dataset) due to memory constraints. How do you moderate this substantial temporal sampling with the goal of preserving temporal dynamics, particularly given RNNs' inherent purpose of capturing sequential dependencies? Could you analyze information loss at different sampling rates, empirical validation that your chosen sampling strategy is optimal, and quantitative metrics demonstrating that temporal patterns are preserved?

7. Concerning the visualizations across LSTM, GRU, and Vanilla RNN, how can you verify that visualization differences genuinely reflect internal model behavior differences rather than artifacts of the visualization method?

8. Specifically in Appendix D.3, how do you distinguish whether the observed 'chaotic pattern' in Vanilla RNN represents actual learning instability rather than limitations of the visualization approach?

9. The paper claims that PCA and t-SNE "failed to capture" certain aspects. How do you verify that the uncaptured information is indeed significant? What evidence supports that MM-PHATE's additional patterns reflect meaningful model dynamics?

10. How does the diffusion kernel accurately capture RNN temporal dependencies? What theoretical guarantees can you provide for this?

11. What is the theoretical justification for using different methods (α-decay vs. Gaussian) for intrastep and interstep affinity calculations in the multiway multislice kernel design?

12. What assumptions are made about the geometric/topological properties of the manifold formed by RNN hidden state dynamics?

13. The authors claim that "changes in entropy coincided with shifts in model performance" (line 323, p.6). However, this claim requires rigorous validation. Could you provide correlation coefficients between entropy changes and performance metrics at each time step, along with statistical tests validating their significance?  Furthermore, can you explain why these entropy changes should reflect performance shifts? 

14. Can you derive statistical bounds for the sampling effects? What convergence guarantees can you provide in the finite sample regime?

15. Given the $O(n^2 s m^2)$ memory complexity, what theoretical approaches could improve this computational burden while maintaining the method's effectiveness?
soundness: 2
strengths: - The paper presents a comprehensive approach to visualizing RNN dynamics that simultaneously captures both temporal and training evolution, addressing a significant gap in existing visualization methods. 
- The authors provide thorough empirical validation through entropy analysis that reveals meaningful patterns in model behavior, particularly in identifying different learning phases and community structures among hidden units. 
- The method demonstrates practical utility by successfully differentiating between different RNN architectures (LSTM, GRU, Vanilla RNN) and providing insights into their learning behaviors without requiring additional validation data.
summary: - This paper introduces MM-PHATE, a novel visualization method designed to analyze hidden dynamics of RNNs during training by capturing both temporal and epoch-wise changes. 
- The method extends PHATE by incorporating a multiway multislice graph representation that preserves community structure among hidden units while tracking their evolution across time steps and training epochs. 
- Through extensive experiments on neurological and human activity recognition datasets, the authors demonstrate that MM-PHATE can identify distinct phases of information processing and compression.
weaknesses: - The authors fail to provide rigorous mathematical foundations, particularly in connecting entropy changes to information bottleneck theory and distinguishing between meaningful compression and information loss. 
- The method suffers from significant scalability issues due to its $O(n^2sm^2)$ complexity, and the necessary sampling of temporal data may compromise the analysis of RNN dynamics. 
- The validation of their method lacks quantitative metrics for evaluating preserved structures and performance correlations, while experimental evidence is limited to small-scale networks and only two datasets. 
- Claims of superiority over baseline methods (PCA, t-SNE) are not sufficiently supported by objective metrics or theoretical guarantees.

---

confidence: 3
contribution: 2
initial_rating: 6
presentation: 3
questions: 1. Line 461-469: This paragraph attempts to address the mismatch in the insight mentioned the previous paragraph (information compression). It presents a hypothesis that the Area2Bump dataset is too complex for the size of LSTMs considered. I believe, it is important to test this hypothesis by training wider/deeper LSTMs. This insight is only valid if you're able to see a similar information compression stage when training with Area2Bump dataset with a sufficiently large model.
2. I find the ability of the presented methods to find "high-information" neural units that are "important" for the task, as presented in figure 4, to be very good for gaining confidence in this method's ability. Repeating this analysis and showing the same for the HAR dataset would be helpful in estabilishing that this phenomena is valid across datasets.
soundness: 3
strengths: - Well detailed presentation of the proposed method as well as prior work.
- Application to two datasets and presentation of gained insights.
- The ability of this method to detect "important" neural units in RNNs (Fig 4c,d).
summary: This paper extends the M-PHATE framework for visualization of the inner time-dependent dynamics of RNNs across training epochs. Example visualizations are provided for two datasets, while also providing insights gained through these visualizations.
weaknesses: The method itself is lacking in novelty as it simply extends the kernel-matrix slicing framework set up in M-PHATE by adding another slicing dimension in the kernel matrix. One could argue that once the slicing framework is established, users could choose to slice across any number of dimensions - epochs, RNN time index, network depth, etc. 

In light of the above, I evaluate the scientific contribution of this paper only in terms of the insights gained through visualizations of RNN dynamics. While insights were presented from the generated visualizations on two datasets, many times these were post-hoc explainations. Some phenomena was only observed in one of the two datasets. Rigorous testing for the generality of these insights is missing.

Meta-review: This work presents MM-PHATE, a visualization technique designed to analyze the internal dynamics of Recurrent Neural Networks (RNNs) during the training process. By extending PHATE with a multiway multislice graph representation, MM-PHATE captures both temporal and epoch-wise changes in RNN hidden states, offering a comprehensive insight into RNN behavior.

The methodology proposed leverages the strengths of PHATE to provide a nuanced understanding of community structure among hidden units and distinct phases of information processing. Experimental evaluations across various datasets, including those from neurological and human activity recognition tasks, demonstrate the effectiveness and broader applicability of MM-PHATE.

However, upon reviewing the feedback from multiple reviewers, it becomes evident that while MM-PHATE shows promise in visualizing RNN dynamics during training, there are several areas that require improvement for the paper to be published. Some key concerns include:
- Theoretical Foundations: Reviewers have noted a lack of rigorous theoretical underpinnings, particularly in connecting the observed entropy changes with information bottleneck theory. Strengthening these foundations would significantly enhance the paper's contribution.
- Scalability and Efficiency: The computational complexity of MM-PHATE, currently at $O(n^2 s m^2)$, poses significant scalability issues. Exploring optimizations or integrating faster versions of PHATE could mitigate these concerns. While the authors have described the potential of new versions of PHATE to reduce this complexity in the rebuttal, the work should make clear how to specifically incorporate those in MM-PHATE.

Overall, the paper introduces an interesting new visualization method for RNN training. Unfortunately, the current version of the manuscript lacks a theoretical grounding and a clear application of the visualization method to improve RNNs. Reviewers acknowledge the superiority of MM-PHATE's visualizations over traditional methods like t-SNE and PCA but emphasize the need for addressing these limitations.

Given these considerations, my overall recommendation leans towards a borderline reject.

Decision: Reject