Graph Convolutional Network with Connectivity Uncertainty for EEG-based Emotion Recognition

We represent an undirected graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ with node set $\mathcal{V}=\{v_{1},v_{2},\ldots,v_{N}\}$ and edge set $\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$. The edges in the graph may be alternately represented using an adjacency matrix ${\bf A}\in\mathbb{R}^{N\times N}$, where ${\bf A}_{ij}=0$ indicates that node $i$ and $j$ are not connected and $N=|\mathcal{V}|$ is the number of nodes. Let the node attribute matrix be ${\bf X}\in\mathbb{R}^{N\times f}$, where $f$ denotes the input feature dimension and each node may be endowed with a $f$-dimensional node attribute vector ${\bf X}_{i:}$. Let ${\bf D}=diag(d_{1},\ldots,d_{N})$ be the diagonal degree matrix of ${\bf A}$ with $d_{i}=\sum_{j}{\bf A}_{ij}$.

Spectral Graph Convolution [35] - Spectral graph theory [48] study the graph topology property by means of the eigenvalues and eigenvectors of the Laplace matrix of graphs. Laplacian matrix is defined as $\bf L=D-A$ and its normalization matrix is $\bf\widetilde{L}=I-D^{-1/2}AD^{-1/2}$, which is a semi-definite matrix with eigenvalues as $\bf\Lambda$ and eigenvectors as $\bf U$. Thus, $\bf\widetilde{L}=U\Lambda U^{T}$. Similar to the role of CNNs in image processing, spectral convolution on graphs can be understood as a filtering operation in the spatial domain. It is achieved by performing a multiplication between the raw graph signal $\mathbf{X}$ and the filter kernel $\mathbf{G}$ in the spectral domain, utilizing techniques such as Fourier transform [35] or wavelet transform [49, 50]. This spectral convolution operation allows us to extract meaningful features from graph-structured data, analogous to the way CNNs extract features from images.

where $\odot$ denotes the element-wise multiplication and ${\bf\hat{G}}=\text{diag}(\hat{g}_{1},\ldots,\hat{g}_{N})$ denotes a diagonal matrix of spectral filter coefficients. To avoid the high computation overload of matrix decomposition of $\bf\widetilde{L}$, [35] approximated $\bf\hat{G}$ with $K$-order Chebyshev polynomials $T_{k}(\bf\hat{\Lambda})$:

and the corresponding filter coefficients $\hat{g}(\lambda_{j})=\sum_{k}\theta_{k}\lambda_{j}^{k}$, where ${\bf\hat{\Lambda}}=\frac{1}{2\lambda_{max}}{\bf\Lambda}-{\bf I}$ is a scaled Laplacian normalization, $\lambda_{max}$ denotes the largest eigenvalue of $\bf L$ and $\theta\in\mathbb{R}^{K}$ is a learnable vector of Chebyshev coefficients, $k\in\{0,\dots,K\}$ means the $k$-th term of a $K$-order polynomials.

Graph convolution neural network (GCN) [51] further approximate the Eq. (2) with $K=1,\theta_{0}=2$, and $\theta_{1}=-1$, so that the convolution operation in Eq. (1) translate into $\bf X\ast G=(I+D^{-1/2}AD^{-1/2})X$. Each GCN layer transformation function $\bf H$ is defined as:

where ${\bf\tilde{A}}_{sym}={\bf\tilde{D}}^{-\frac{1}{2}}{\bf\tilde{A}}{\bf\tilde{D}}^{-\frac{1}{2}}$ is a normalized adjacency matrix, $\bf\tilde{D}=I+D,\tilde{A}=I+A$ is a renormalization trick introduced in GCN, ${\bf H}^{(0)}={\bf X}$, ${\bf W}^{(l)}$ is the convolutional parameter of each layer to realize the feature transformation, $\sigma(\cdot)$ denotes a nonlinear activation function and $k$ denotes the layer number. The normalized Laplacian ${\bf\tilde{L}}_{sym}={\bf I-\tilde{D}}^{-\frac{1}{2}}{\bf\tilde{A}\tilde{D}}^{-\frac{1}{2}}={\bf\tilde{U}\tilde{\Lambda}\tilde{U}^{T}}$ and its corresponding graph filter can be described as $\hat{g}(\hat{\lambda}_{i})={(1-\hat{\lambda}_{i})}^{L}$, where $\hat{\lambda}_{i}$ are eigenvalues of ${\bf\tilde{A}}_{sym}$.

For brevity and without loss of generality, we denote $\bf A$ as the adjacency matrix with a self-loop and $\bf D$ as the degree matrix with a self-loop.

The aggregation progress of neighborhood information can be seen as a graph filtering step over node features with a filter expressed in Eq. (2). Obviously, the preserved spectral component in $\bf\hat{G}_{\Lambda}$ is dominated by the eigenvalues of $T_{k}(\bf\hat{\Lambda})$, i.e. ${\bf\hat{\Lambda}}=\{\hat{\lambda}_{1},\ldots,\hat{\lambda}_{K}\}$ and $\hat{\lambda}_{1}>,\ldots,>\hat{\lambda}_{K}$, which in turn represent the low to high-frequency components. Currently, practical models are usually shallow (number of layers $K$ is 2 - 4), as such a choice offers the best empirical performance which is later proved attributed to the low-pass characteristic. Besides, the transition matrix will finally collapse to an independent matrix of input without a linear interposition in Eq. (3), $\lim\limits_{k\to\infty}{\bf H}^{(k)}={\bf H}^{(\infty)}{\bf\leftarrow}{\bf\tilde{A}}_{sym}^{\infty}{\bf H}^{(0)}$. And the final decision derived from such an undiscriminating, ${\bf Y}_{pred}={\rm softmax}(\sigma({\bf H}^{(L)}{\bf W}))$ is merely determined by the $L$-hop node attributes with $L$ denotes the total layers, provided that the graph is irreducible and aperiodic. More specifically, the learned representations will lose discriminative information provided by nodes with different topological and feature characteristics as the model goes deeper.

In practical applications, each EEG channel is regarded as a node, with node vectors encapsulating the computed features across various frequency bands, specifically, the delta, theta, alpha, beta, and gamma bands. Our central proposition rests on leveraging the information derived from both the inherent topology of the graph in the spatial domain and the node features in the time-frequency domain to overcome prevalent challenges such as the identification of topological connections and distinct effective frequency band responses. This two-pronged approach allows us to pinpoint crucial information to be assimilated into the propagation of knowledge during the design phase of the GCN.

A crucial step in this procedure is the introduction of an edge predictor (EP), which computes the connection probabilities between any two nodes, paving the way for node feature learning. The edge predictor can be formally expressed as $f_{ep}:\boldsymbol{\rm A,X}\rightarrow\boldsymbol{\rm Z}$, where $\boldsymbol{\rm Z}$ denotes an edge existence mask. We leverage the mask matrix ${\bf Z}$ to establish deterministic links to high-scoring non-edges and remove poorly scored linked edges. The propagation procedure is then reformulated as the following metric:

It is easy to observe that implementing node-wise and feature-wise uncertainty quantization leads to different probabilistic interpretations. To provide a comprehensive measure of model uncertainty, we further extend ${\bf Z}^{(l)}={\bf Z}_{uv}^{(l)}\cdot{\bf Z}_{f}^{(l)}$ as a 3-dimensional matrix, where ${\bf Z}_{uv}^{(l)}\in\mathbb{R}^{N\times N\times 1}$ denotes the connection of any two nodes $u,v$, and ${\bf Z}_{f}^{(l)}\in\mathbb{R}^{N\times f_{l}}$ denotes the mask vector of layer $l$ on features of each node. We can rewrite and reorganize Eq. 4 in a feature-wise view as

where

that means we allocate each frequency component with a different topological connection strategy for a better group of useful information. ${\bf Z}^{(l)}[:,i]$ is a binary mask sampled from a Bernoulli distribution with a success rate as $p_{l}$ for each layer. The success rate $p_{l}$ is adaptive for different representations, guaranteeing the effectiveness of topology construction.

To generate an uncertain binary mask ${\bf Z}^{(l)}$, we endeavor to leverage the power of Bayesian inference to train our predictor, using the predictive posterior derived from training data. Intriguingly, we illustrate in the context of this study that the dynamic process of adaptive connection learning can be elegantly transposed from the realm of output feature space to the more abstract parameter space, thereby casting it as a credible Bayesian process encapsulated within the Graph Convolutional Network (GCN) framework.

To further clarify our approach, we revisit Eq. 4 and reinterpret it from a node-wise perspective of a GCN layer, which can be described as:

where

In this context, we introduce the notation ${\bf W}^{(l)}_{uv}=\text{diag}({\bf Z}_{uv}^{(l)}){\bf W}^{(l)}$, thereby facilitating the transference of the learning process for the uncertainty-based edge predictor to the acquisition of a weighted adjacency matrix along the edges $e=(u,v)\in\mathcal{E}$ for each respective layer.

The model parameters that require optimization are denoted as $\boldsymbol{\omega}=\{\boldsymbol{\omega}_{l}\}_{l=1}^{L}$, with $\boldsymbol{\omega}_{l}=\{{\bf W}_{e}^{(l)}\}_{e=1}^{|\mathcal{E}|}$ serving as the weight set for the $e$-th edge. Here, $|\mathcal{E}|$ represents the number of edges. Initially, each row of ${\bf W}_{e}^{(l)}$ is assigned a distribution according to $p(\boldsymbol{\omega})$. Moreover, we postulate a vector ${\bf m}_{l}$ of dimensions $f_{l}\times f_{l-1}$ for each layer. Consequently, the predictive probability of the deep GNN model can be expressed as:

However, the posterior distribution $p(\boldsymbol{\omega}|{\bf A,X})$ proves intractable. In light of this, we adopt a variational interpretation, where a binary mask (each row of $\boldsymbol{\omega}$) can be viewed as an approximate distribution $q_{\theta}(\boldsymbol{\omega})$ that approximates the posterior distribution $p(\boldsymbol{\omega}|{\bf A,X})$. We define $q_{\theta}(\boldsymbol{\omega})$ as follows:

Here, $\theta=\{p_{l},{\bf m}_{l}\}_{l=1}^{L}$ represents the variational parameters, and ${\bf m}_{l}$ denotes the mean weight matrices. The variable ${\bf z}_{l,e}=\{0,1\}^{f_{l}}$ corresponds to the removal of features, leading to the omission of nodes in layer $l-1$ as connected neighbors within the graph ${\bf A}^{(l)}$, solely when the entire set of features is eliminated.

To optimize the aforementioned model, we introduce the Kullback-Leibler (KL) divergence ${\rm KL}(q_{\theta}(\boldsymbol{\omega})||(p(\boldsymbol{\omega}))$ in order to minimize the discrepancy between the approximate posterior $q_{\theta}(\boldsymbol{\omega})$ and the prior distribution $p(\boldsymbol{\omega})$. The discrete quantized Gaussian prior distribution is commonly employed to analytically evaluate this intractable KL divergence. Considering the distribution of the edges, we have $q_{\theta}(\boldsymbol{\omega})=\prod_{l=1}^{L}\prod_{e=1}^{|\mathcal{E}|}q_{\theta_{l}}({\bf W}_{e}^{(l)})$. By approximating each edge distribution as $q_{\theta_{l}}({\bf W}_{e}^{(l)})=p_{l}\delta({\bf W}_{e}^{(l)}-0)+(1-p_{l})\delta({\bf W}_{e}^{(l)}-{\bf m}_{l})$, we can express the general KL term as:

Consequently, the edge predictor loss $\mathcal{L}_{ep}$ can be approximated as:

where $\mathcal{H}$ denotes the entropy of a Bernoulli random variable with success rate $p_{l}$, defined as:

Given a fixed edge connection probability $p_{l}$, the entropy term remains unaffected by model weights and can be viewed as a constant regularization term during optimization, thereby rendering it negligible. However, the connectivity sampling probability inherently relies on the specific functional structure of each frequency component. The minimization of the KL divergence term equates to maximizing the entropy term with a probability of $1-p_{l}$, which is optimally achieved when $p_{l}$ approximates 0.5. It becomes apparent that constant proportion sampling fails to meet the requirements for intricate functional connections and may introduce extraneous information. In order to train $p_{l}$ during the model training process, we must compute the derivative of the entropy objective, which presents a formidable challenge.

Several commonly used estimators, such as Reinforce [52] or pathwise derivative re-parametrization methods [53, 44, 54], are inapplicable due to the discrete nature of the masks. Rather than sampling the random variable from a discrete Bernoulli distribution, we employ an approximation by substituting the discrete distribution with its continuous relaxation. This can be seen as a relaxation of the ”max” function in the Gumbel-max trick to a ”softmax” function. The transition from a discrete Bernoulli random variable $\bf z$ to a continuous variable $\bf{\tilde{z}}$ is reparametrized as $\bf{\tilde{z}}=g(\theta,\epsilon)$, where $\theta$ represents the model parameters and $\epsilon$ is a random variable independent of $\theta$. Assuming that $\epsilon$ follows a uniform distribution, i.e., $\epsilon\sim\text{Unif}(0,1)$, we utilize the Sigmoid function to constrain the continuous value within the interval [0, 1]:

where $t$ is a temperature parameter. This distribution assigns most of its mass to the boundaries of the interval 0 and 1. With this concrete relaxation of the adjacency masks, we can now optimize the connection probability.

Our investigation enlisted two established EEG datasets, the SJTU emotion EEG dataset (SEED) [4] and its successor SEEDIV dataset [5], in the execution of emotion recognition tasks. These datasets encompass EEG recordings from 15 youthful participants, evenly distributed by gender. The data, collected via a 62-channel ESI NeuroScan System as the subjects viewed emotionally charged film clips, were gathered over three sessions held on different days under uniform conditions to ensure data integrity. The raw EEG signals were then preprocessed, down-sampled from 1000 Hz to 200 Hz, and cleansed of artifacts stemming from eye movements (EOG) and muscle activity (EMG).

The SEED dataset comprises 15 film clips of around 4 minutes each, carefully crafted by psychologists to evoke specific emotional states - Positive, Neutral, or Negative. Hence, each participant in the SEED dataset produced 45 distinct EEG signal trials. The SEEDIV dataset, following a similar data acquisition procedure, increased the number of trials by featuring 24 film clips that lasted roughly 2 minutes each and elicited four types of emotions: Neutral, Sad, Fear, and Happy. This resulted in each participant contributing 72 EEG recording trials in the SEEDIV dataset.

In the pursuit of equitable comparisons, we conducted both subject-dependent and subject-independent classifications on the SEED and SEED-IV datasets, in adherence with protocols laid out in previous investigations [65, 4, 10, 11, 20, 8]. In the subject-dependent trials on the SEED dataset, we designated the initial 3 trials from each emotional category (making up a total of 9 trials) for training purposes, while the residual 6 trials were assigned to the testing set, thereby ensuring an even distribution across categories. Notably, these experiments drew upon two sessions from each subject. We reported the final outcome as the mean accuracy across all 15 subjects.

Turning to the subject-dependent trials on the SEED-IV dataset, we allocated the first 4 trials from each emotional category (a total of 16 trials) for training, whereas the remaining 8 trials served as the testing set. All three sessions were incorporated into the evaluation of model performance. As with the SEED dataset, we reported accuracy over all 15 subjects. Regarding the subject-independent trials, we espoused the leave-one-out methodology. Specifically, cross-validation was performed by excluding one session from each subject in the case of the SEED dataset, and all sessions from each subject in the case of the SEED-IV dataset. We evaluated model performance across all 15 trials.

We conducted our experiments using the PyTorch Geometric toolbox and trained the models on a single NVIDIA Tesla A100 GPU, utilizing the Adam optimizer. To determine the optimal hyperparameters, we performed a search on the validation set, exploring the following ranges: (a) initial learning rate within the range [1e-3, 3e-2]; (b) the number of convolutional layers within range {2, 4, 6, 8} and hidden feature dimension within range {32, 64, 128}. We set the temperature parameter $t$ in the concrete distribution as 0.67, the Beta distribution weight $\psi$ as 2, and the batch size as 32. The hyper-parameters were hand-tuned with the best performance over the validation set. A cosine annealing learning rate scheduler was adopted for all experiments. We adopt an early stopping strategy when the validation loss did not decrease for ten consecutive epochs.

Following the procedure of Algorithm 1, we compare the proposed model with a quantity of popular used models: (1) traditional machine learning classifying models, e.g. SVM [5], DBN [4]; (2) non-topological deep learning methods, e.g. CNN alone or combined with RNN, such as BiDANN [61], BiHDM [63] and R2G-STNN [62]; (3) recently proposed GNN models: vanilla GCN (vGCN) [51], DGCNN [10], RGNN [11], GCB-net with BLS [64] and V-IAG [20]. To further explore the relation of spectral information and emotion state, we conduct these experiments from the separate and total five frequency bands, known as $\delta$ band (1-4 Hz), $\theta$ band (4-8 Hz), $\alpha$ band (8-14 Hz), $\beta$ band (14-30 Hz) and $\gamma$ band (30-50 Hz). For a fair comparison, we adopt the pre-extracted DE, PSD, DASM and RASM features after linear dynamic systems (LDS) as our input source data.

Table 1 shows the performance of our model and the SOTA methods. In view of the DE feature is more effective in emotion recognition, we only take experiments with DE feature on SEEDIV, with results shown in Table 2. Our model performs on par or better than other methods both on separate and total frequency bands on all features. We analyze the results from three aspects: (1) input source - DE feature are experimentally proved as the most effective metric; (2) classifier - deep learning methods among which CNN is used to acquire spatial attention and RNN is used to learn temporal relation, outperformed traditional machine learning methods for a lot. While GCN-based methods utilize the EEG topological structure and make a more reasonable channel grouping, finally present the best performance; (3) spectral effectiveness: leverage features on $\beta$ and $\gamma$ band alone are more effective for emotion activation but greatly under-performed than using features on whole bands, somehow indicating that high-frequency information is more active with emotional activities. Our model obtained the best performance on both datasets. We attribute such a performance gain to the leverage of uncertainty-guided instance-wise adjacency matrix and the improvement of model ability on the spectral richness and alleviation to over-smoothing.

Physiological signals, such as EEG, present severe individual variation and greatly hinder the model generalization. Thus, we conduct the subject-independent emotion classification to further verify our model stability. We use DE feature as the exclusive input for its remarkable discriminating property. Table 3 shows the performance of three categorical classifiers on SEED and SEEDIV. Note that for brief comparison, we conduct the test using features from all five bands on SEEDIV following previous literature settings. Not surprisingly, despite numerical value presenting a similar tendency as subject-dependent experiments in the view of all three analysis aspects, all accuracy suffered a great decrease. Our method achieves the best performance in SEEDIV and comparable results in the SEED datasets, which verifies the generalization of CU-GCN with subject-independent EEG emotion recognition. R2G-STNN leveraged a bidirectional long short term memory (BiLSTM) to learn the region to global spatial-temporal information. BiHDM proposed a four-directed RNN based on two spatial orientations to obtain the discrepancy information between two separate brain regions and used a domain discriminator to generate the domain-invariant feature. RGNN eliminates the individual invariant by introducing a similar domain adversarial training pattern. V-IAG was believed to promote cross-subject performance by designing a more reasonable adjacency matrix that contains an instance adaptive branch and a variational branch to learn the deterministic graphs and an uncertain graph. All these methods solve the subject-independent problem by designing a data-driven node connection graph from both regional and global space. It indicates that the key point of emotion recognition should lie in the precise combination of EEG channels.

GCN-based methods from Table 1, 2, 3, including DGCNN, RGNN, and V-IAG, have demonstrated superior performance in EEG-based emotion recognition tasks compared to other approaches. DGCNN pioneered the utilization of two-layer vanilla GCN to process EEG signals by leveraging graph structure. Building upon this, RGNN further emphasized the importance of node aggregation in graph structure learning and employed a two-layer Spectral Graph Convolution (SGC) as the backbone. However, both vanilla GCN and SGC suffer from the over-smoothing problem, which arises as the network depth increases, leading to the loss of local information. V-IAG employed an eight-layer GCN architecture that preserved features from local to global regions by aggregating the outputs of each layer for the final decision.

In light of our previous theoretical analysis, our enhanced GCN (eGCN) is designed to address the over-smoothing problem while enhancing spectral information. In this subsection, we perform ablation studies by replacing eGCN with SGC and GCN to validate our hypothesis. Table 4 presents the results, where ”eGCN $\rightarrow$ GCN(2)” indicates the replacement of eGCN with a two-layer GCN, and similarly for SGC. As shown in Table 4, it is evident that eGCN achieves significantly better results than the common benchmark methods.

Furthermore, we conduct additional experiments by varying the number of layers ($L$ corresponding to the number of aggregated filters in CU-GCN). The results in Table 5 consistently demonstrate that our method outperforms others, with the best performance achieved when $L=6$, equivalent to six layers. These findings indicate that eGCN effectively mitigates the over-smoothing problem as the network depth increases. Notably, the aggregation of information from six-hops neighborhoods performs optimally, potentially due to ensuring the global connection between two hemisphere brain regions, as the relative node distance in a 2D squared matrix is nine. Overall, these results highlight the superior ability of CU-GCN to aggregate information over larger neighborhoods compared to other methods.

To compare the effectiveness of our uncertainty-guided edge predictor (EP) with commonly used graph construction approaches, we trained CU-GCN for subject-independent emotion recognition on the SEED and SEEDIV datasets. The second part of Table 5 presents the performance of CU-GCN with and without EPM using different adjacency matrices: distance-based (Dist), coherence-based (Coh), and random adjacency matrices. We observed that the inclusion of EP consistently yielded better results compared to other methods. Interestingly, traditional graph construction methods did not exhibit significant differences, suggesting that EEG signals in emotion recognition possess unique characteristics. Our method offers several distinct advantages: (a) It can be applied even when the physical locations of electrodes are unknown, providing flexibility in practical scenarios. (b) It captures dynamic brain connectivity patterns instead of relying solely on spatial sensor information, which is particularly desirable for emotion recognition tasks. (c) It enables adaptive selection of spectral features in EEG, facilitating precise categorization.

In our study, we conducted an investigation into the impact of using a learnable Bernoulli prior distribution compared to a constant edge dropout method. The learnable Bernoulli success rate $p$ in our CU-GCN model is adaptively optimized during training. The results presented in Table 6 provide compelling evidence for the effectiveness and necessity of this approach, as CU-GCN with the learnable Bernoulli prior consistently outperforms the constant $p$ value method in terms of classification accuracy. Moreover, it is worth noting that even without the learnable Bernoulli prior, our CU-GCN approach still achieves competitive performance, often surpassing the results obtained by state-of-the-art methods. This highlights the robustness and effectiveness of CU-GCN in capturing important graph structure information and achieving accurate emotion classification, even in the absence of the learnable prior distribution. Overall, these findings emphasize the significance of incorporating the learnable Bernoulli prior distribution in CU-GCN, as it contributes to improved classification accuracy and enhances the model’s ability to leverage graph structure for emotion recognition tasks.

In our CU-GCN framework, we employ the graph mixup data augmentation technique to enhance the probability confidence and facilitate the learning process in the presence of noisy labels. We conducted experiments to evaluate the effectiveness of this approach, and the results are presented in Table 4. The findings indicate that training with graph mixup leads to an improvement of 0.2$\sim$0.3% in performance compared to training with the original dataset alone. Furthermore, the confusion matrices depicted in Figure 4 for the SEED and SEEDIV datasets demonstrate the substantial enhancement achieved by our data augmentation methods in addressing the challenges posed by noisy labels.