Compressive Sensing of Multichannel Electroencephalogram Signals Based on Nonlocal Low-Rank and Cosparse Priors

Abstract

Recent studies have shown that by using channel-correlation and cosparsity in a centralized framework, the accuracy of reconstructing multichannel EEG signals can be improved. A single-channel electroencephalogram (EEG) signal is intrinsically non-sparse in both the converted and raw time domains, which presents a number of important issues. However, this is ignored by contemporary compressive sensing (CS) algorithms, resulting in less recovery quality than is ideal. To address these constraints, we provide a novel CS method that takes advantage of Nonlocal Low-Rank and Cosparse priors (NLRC). By utilizing low-rank approximations and block operations, our method aims to improve the CS recovery process and take advantage of channel correlations. The Alternating Direction Method of Multipliers (ADMM) are also used to efficiently solve the resulting non-convex optimization problem. The outcomes of the experiments unequivocally demonstrate that by using NLRC, the quality of signal reconstruction is significantly enhanced.

1. Introduction

Electroencephalogram (EEG) signals are generated by the bioelectrical activity of the brain and represent a critical tool in neuroscience, clinical medicine, and brain–computer interface research [1,2,3]. These signals are primarily generated by the collective electrical activity of neurons, specifically through postsynaptic currents resulting from synaptic transmission [4]. However, recent computational studies have highlighted that while postsynaptic currents are indeed a contributing factor, their individual amplitudes are relatively small. Consequently, the EEG signal is thought to be a composite of multiple electrical activities, including but not limited to postsynaptic currents, and reflects the summation of these activities across a large number of neurons.

EEG signals are non-invasive, and can be recorded from the scalp using electrode arrays. These signals carry rich information about brain activity, including cognitive processes, sensory and motor functions, and even sleep patterns. Understanding the generation and detection of EEG signals is crucial for advancing our knowledge of brain function and for developing new diagnostic and therapeutic tools.

In the practical development of modern medicine, in order to obtain and analyze the EEG signals of the human body, the electroencephalography and other equipment are often used to monitor and record the EEG activity. To ensure more accurate and reliable data, portable EEG monitors detector is generally multi-channel, long-term collection of EEG data, and the amount of data are very large. It requires high energy for transmission thus it cannot be applied in practical application scenarios [5,6,7].

Using a Compressive Sensing (CS) approach may be a potential solution to this problem, which is based on the sparsity of the signals [8,9]. It observes the signal through linear projection and adopts non-uniform random sampling. The number of samples needed is directly impacted by a signal’s sparsity. Because the data sampling and compression can be performed at the same time, compressive sensing reduces the amount of redundant information that builds up, making data easier to understand and more compressed.

Conventional sparse synthesis models [8,10,11] are predicated on the idea that signals are sparse within a specific transform domain, where past information is essential to reconstructing the computer system. Furthermore, the analysis sparsity model [9,12,13,14,15] assumes that sparsity is attained when an analysis operator is applied to the signal. Interestingly, studies have demonstrated that when it comes to reconstructing multi-channel electroencephalogram (EEG) signals, dense models perform better than sparse models. Recently, the low-rank prior [16,17,18,19], as a useful prior, has been widely used in multichannel EEG signals, and it is used to exploit the correlation between multiple channels. According to the singular value decomposition process of matrix rank, low rank regularization not only constrains the sparsity of each column (single channel signal), but also constrains the sparsity of each row (correlation between multiple channel signals). The benefits of combining channel correlation and joint influence sparsity into a single, unified strategy have been highlighted in recent research, which greatly enhances recovery outcomes.

However, the single-channel electroencephalogram (EEG) signal is essentially non-sparse in the converted domain, as well as in the original time domain, resulting in a low level of overall sparsity in the multichannel EEG signal. Although nearly all of the existing compressive sensing algorithms make use of the sparsity or similarity of the entire channel signal, they are unable to produce a recovery quality that is competitive. In the field of EEG signals denoising, nonlocal mean (NLM) methods [20], which calculate the weighted sum of a block, can efficiently exploit the accurate adjacent relationship between different blocks. In NLM, each point is given a weight based on how similar it is to two other points: one in its own block, and one in a nearby block. After that, the weighted total is calculated, which reduces noise in the process. Pursuing more precise recovery outcomes that can benefit from cosparsity prior, low-rank regularization, and nonlocal operation [7] is better.

In this work, we solve the problems caused by channel correlation and cosparsity in multichannel EEG signals by introducing a novel technique termed the Non-local Low-Rank and Analysis Sparsity Prior (NLRC). We recommend using low-rank approximation and block operations as regularization strategies to take advantage of the channel correlations in the context of compressive sensing (CS) recovery. Furthermore, we solve the resulting nonconvex optimization model by utilizing the alternating direction method of multipliers (ADMM) [21,22].

This paper is organized as follows. In Section 2, we briefly introduced sparse synthesis model, cosparse analysis model and cosparsity and low-rank property-based model. We present a thorough explanation of the NLRC approach in Section 3, emphasizing its application to cosparsity and accurate low-rank characteristics. Furthermore, we provide the suggested model’s general optimization algorithm. In Section 4, we conduct through numerical experiments to verify the algorithm’s efficacy, and finally make inferences based on our results.

2. Related Works

2.1. Sparse Synthesis Model and Cosparse Analysis Model

The sparse synthesis model assumes that a signal is sparse in a transform domain. Assuming that y is the random measurement and $\mathsf{\Phi }\in R^{M\times N}$ ($M<N$) is the sampling matrix, the sampling model is

$$y=\mathsf{\Phi }x$$

(1)

The signal $x\in R^N$ can be represented as $x=\mathsf{\Psi }s$, where $\mathsf{\Psi }$ is the transform domain and only contains few nonzero elements. We then proceed to design the model of l0 minimization.

$$\underset{s}{min}{∥s∥}_{0}s.t.y=\mathsf{\Phi }\mathsf{\Psi }s$$

(2)

where ${∥s∥}_0$ counts the number of nonzero elements of s. To tackle the aforementioned problems, a number of successful strategies [8,10,11] have shown to be highly effective. In contrast to the sparse synthesis model, the cosparse analysis model [9,12] assumes that a sparse vector is produced when a signal is multiplied by an analysis operator. The solving problem can be represented as

$$\underset{x}{min}{∥Ox∥}_{0}s.t.y=\mathsf{\Phi }x$$

(3)

where $O\in R^{K\times N}$ ($K\ge N$) is the analysis operator and $a=Ox$ is the cosparse vector. The above problem can be efficiently solved by many proposed methods [9,13,14]. The Analysis-by-Synthesis method (ABS) [13] is a convex relaxation that rewrites the problem in the form of a synthesis-based problem. The Analysis l1-minimization algorithm [14] amounts to using l1-norm instead of l0-norm. Greedy Analysis Pursuit (GAP) [9] is a greedy algorithm that is the analysis duality of Orthogonal Matching Pursuit.

The effectiveness of the cosparse analysis method outperforms the sparse synthesis strategy in compressive sensing recovery of multichannel EEG signals, as demonstrated in reference [16]. First, whereas the cosparse analysis model allows consistency between columns of the analysis operator, the sparse synthesis model requires imposing incoherence limitations on the sampling matrix. Second, the cosparse analysis model immediately estimates the signal, offering a more simple and efficient method than the sparse synthesis model, which estimates the sparse vector first, followed by the signal.

2.2. Cosparsity and Low-Rank Property Based Recovery

Motivated by the low rank regularization not only constraining the sparsity of each column (single channel signal), but also constraining the sparsity of each row (correlation between multiple channel signals), many researchers proposed the recovering model based on cosparsity and low-rank priors (SCLR) [16].

$$\underset{X}{min}{∥r\left(OX\right)∥}_0+\lambda rank\left(X\right)s.t.Y=\mathsf{\Phi }X$$

(4)

where $X\in R^{N\times C}$ (C is the number of the channels) and $r\left(·\right)$ sequentially puts all the columns into a column vector. Since the previously mentioned optimization issue is NP-hard, convex surrogate functions must be used in order to solve it successfully.

3. Nonlocal Low-Rank and Cosparse Priors for Multichannel EEG Signals

In this article, we present a new technique for regularizing covariances and inter-channel correlations in multichannel EEG signals. We present our technique, which we call the Non-local Low-Rank and Sparse (NLRC) prior. It employs inter-channel correlations to improve compressive sensing (CS) recovery by utilizing block operations and low-rank approximations.

Using the same methods, we first divide each channel signal into a concatenation of non-overlapping blocks. Then, for each i-th block $x_i\in R^{d⁕1}$, we search k-nearest similar blocks throughout all channel signals. Then, for each i-th block, we search k-nearest similar blocks satisfying ${∥x_i-x_{i_{\mathrm{c}}}∥}_2<\epsilon$ ($\epsilon$ is a predefined threshold) throughout all channel signals. Low-rank qualities are basic on the final data matrix $X_i\in R^{d⁕k}$ if these signal blocks have similar topologies (see Figure 1).

Next, we design the following optimization model for multichannel EEG signal compressive sensing recovery,

$$\underset{X}{min}{∥r\left(OX\right)∥}_0+\lambda \sum _i\mathrm{rank}\left(P_{i}X\right)s.t.Y=\mathsf{\Phi }X$$

(5)

We introduce the lq norm and weighted Schatten-p norm as appropriate non-convex surrogate functions to address the difficult NP-hard optimization problem. Since the model involving the l0 functions of norm and matrix rank is non-convex, it does not have a global optimal solution, and only a local optimal solution. In addition, the optimization problem is a combinatorial optimization problem that requires finding the optimal solution among all possible sparse solutions, which leads to a very high computational complexity. We take inspiration from the success of the lq norm in improving the accuracy of sparse synthesis models [23,24] and the efficiency of the weighted Schatten-p norm in denoising images via low-rank matrix recovery [25,26].

Then, the problem can be described as

$$\underset{X}{min}{∥r\left(OX\right)∥}_q^q+\lambda \sum _i{∥P_{i}X∥}_{W,S_p}^{p}s.t.Y=\mathsf{\Phi }X$$

(6)

where ${∥·∥}_q^q$ norm computers the total of the absolute values of the entries to the power of q, and ${∥·∥}_{S_p}^p$ sums all the singular values of X to the power of p.

The constrained optimization problem above is difficult to solve, so we introduce ADMM, a widely used method in compressive sensing [21,27,28], to divide the problem into some sub-problems which are easier to address. By adding a set of auxiliary variables $\left\{A,B\right\}$, the recovery problem can be reformulated as

$$\underset{X}{min}{∥r\left(A\right)∥}_q^q+\lambda \sum _i{∥B_i∥}_{W,S_p}^{p}s.t.Y=\mathsf{\Phi }X,A=OX,P_{i}X=B_i$$

(7)

Its augmented Lagrangian form [28] is as follows:

$$\begin{array}{c}\underset{X,A,B}{min}{∥r\left(A\right)∥}_q^q+\lambda {\displaystyle \sum _i}{∥B_i∥}_{W,S_p}^p+{\displaystyle \frac{{\beta }_1}{2}}{∥\mathsf{\Phi }X-Y+{\displaystyle \frac{f_1}{{\beta }_1}}∥}_F^2+\hfill \\ {\displaystyle \frac{{\beta }_2}{2}}{∥OX-A+{\displaystyle \frac{f_2}{{\beta }_2}}∥}_F^2+{\displaystyle \frac{{\beta }_3}{2}}{\displaystyle \sum _i}{∥P_{i}X-B_i+{\displaystyle \frac{g_i}{{\beta }_3}}∥}_F^2\hfill \end{array}$$

(8)

where $f_1$, $f_2$ and $g_i$ are Lagrangian multipliers. In this paper, we employ ADMM to solve the resulting minimization problem.

$$\begin{array}{c}{X}^{k+1}=\underset{X}{min}{\displaystyle \frac{{\beta }_1}{2}}{∥\mathsf{\Phi }X-Y+{\displaystyle \frac{f_1}{{\beta }_1}}∥}_F^2+{\displaystyle \frac{{\beta }_2}{2}}{∥O^{k}X-A^k+{\displaystyle \frac{f_2}{{\beta }_2}}∥}_F^2\hfill \\ +{\displaystyle \frac{{\beta }_3}{2}}{\displaystyle \sum _i}{∥P_i^{k}X-B^k+{\displaystyle \frac{g_i}{{\beta }_3}}∥}_F^2\hfill \end{array}$$

(9)

$$A^{k+1}=\underset{A}{min}{∥r\left(A\right)∥}_q^q+{\displaystyle \frac{{\beta }_2}{2}}{∥O^{k+1}{X}^{k+1}-A+{\displaystyle \frac{f_2}{{\beta }_2}}∥}_F^2$$

(10)

$$B_i^{k+1}=\underset{B_i}{min}\lambda {∥B_i∥}_{W,S_p}^p+{\displaystyle \frac{{\beta }_3}{2}}{∥P_i^{k+1}{X}^{k+1}-B_i+{\displaystyle \frac{f_3}{{\beta }_3}}∥}_F^2$$

(11)

3.1. X Sub-Problem

The X sub-problem is a quadratic optimization problem admitting a closed-form solution.

$$\begin{array}{c}{X}^{k+1}={\left({\beta }_1{\mathsf{\Phi }}^T\mathsf{\Phi }+{\beta }_2{\left(O^k\right)}^T{O}^k+{\beta }_3{\displaystyle \sum _i}{\left(P_i^k\right)}^T{P}_i^k\right)}^{-1}\hfill \\ \left({\beta }_1{\mathsf{\Phi }}^{T}Y+{\beta }_2{\left(O^k\right)}^T{A}^k+{\beta }_3{\displaystyle \sum _i}{\left(P_i^k\right)}^T{B}_i^k-{\mathsf{\Phi }}^T{f}_1-{\left(O^k\right)}^T{f}_2-{\displaystyle \sum _i}{\left(P_i^k\right)}^T{g}_i\right)\hfill \end{array}$$

(12)

3.2. A Sub-Problem

The A subproblem itself is non-convex, which complicates the search for a global minimizer. However, we may overcome this difficulty by using the iterative reweighting technique outlined in [24]. Let us assume that

$$A^{k+1,t+1}=\underset{A}{min}\sum _i{u}_i^t{a}_i+{\displaystyle \frac{{\beta }_2}{2}}{∥r(O{X}^{k+1}-A+{\displaystyle \frac{f_2}{{\beta }_2}})∥}_2^2$$

(13)

where $u_i^t$ is the weight $q{\left(\left|a_i^t\right|\right)}^{q-1}$ that is computed from the previous iterative $A^{k+1,t}$, and $a_i$ is the i-th value of $r\left(A\right)$. By applying the mathematical frameworks described in [25,26], the solution to this issue can be found directly.

$$a^{k+1,t+1}=max\left(r\left(O{X}^{k+1}+{\displaystyle \frac{f_2}{{\beta }_2}}\right)-{\displaystyle \frac{1}{{\beta }_2}}{u}^t,0\right)$$

(14)

and when $a^{k+1,t+1}$ satisfies the convergence condition, we set $a^{k+1}=a^{k+1,t+1}$.

3.3. B Sub-Problem

The B sub-problem is also a non-convex problem. We use weighted singular value shrinkage [27] to solve this problem, as shown below.

$$B_i^{k+1,c+1}=Umax\left\{\Delta -{\displaystyle \frac{\lambda }{{\beta }_3}}diag\left(w_i^c\right),0\right\}{V}^T$$

(15)

where $U\Delta V^T$ is the Singular Value Decomposition (SVD) of $\left(P_i{X}^{k+1}+{\displaystyle \frac{f_3}{{\beta }_3}}\right)$ and $w_{i,j}^c=p{\delta }_j^{p-1}\left(B_i^{k+1,c}\right)$. ${\delta }_j\left(B^{k+1,c}\right)$ is the j-th singular value of $B^{k+1,c}$. When $B^{k+1,c+1}$ satisfies the convergence condition, we set $B^{k+1}=B^{k+1,c+1}$.

The procedure is shown in Algorithm 1.

Algorithm 1 Compressive sensing of multichannel EEG signals based on nonlocal low-rank and cosparse priors

Input:$\mathsf{\Phi }$, $O\in R^{K\times N}\left(K=2\times N\right)$, ${\beta }_1={\beta }_2={\beta }_3=1$$,\lambda =1$$,\tau =0.05$$, f_1=zeros\left(M,C\right)$,          $f_2=zeros\left(K,C\right)$, $f_3=zeros\left(N,C\right)$;    while stopping criteria unsatisfied do       (a) Compute X via Equation (12);       (b) Compute A by computing Equation (14);       (c) Compute B via Equation (15);       (d) Update Lagragian multipliers:             $f_1\leftarrow f_1-\tau {\beta }_1\left(Y-\mathsf{\Phi }{X}^{k+1}\right)$;             $f_2\leftarrow f_2-\tau {\beta }_2\left(A^{k+1}-O^{k+1}{X}^{k+1}\right)$;             $g_i\leftarrow g_3-\tau {\beta }_3\left(B_i^{k+1}-P_i^{k+1}{X}^{k+1}\right)$;    end while Output: final reconstructed signal $\stackrel{⌢}{X}$;

4. Experimental Results

We give numerical tests in this section to assess the NLRC method’s performance. We contrast NLRC with a number of cutting-edge techniques, such as the BSBL approach [7], ADMM-based SCLR (SCLR-A) [16], SCLR based on the interior point method (SCLR-I) [16], and SCLR-based lq norm and Schatten-pnorm (LQSP) [29]. One well-known use of the BSBL approach is the recovery of EEG signals in the sparse synthesis CS model. In compressed sensing reconstruction, all compared approaches take advantage of cosparsity and channel correlation: SCLR-I, SCLR-A, and NLRC. LQSP has shown the competitive reconstruction results at present. The LQSP method directly exploited the overall cosparsity of the multichannel EEG signal and low-rank property, and took lq norm and Schatten-p norm as surrogate functions for l0 norm and matrix rank.

The BCI III dataset 1 [30] and the CHB-MIT scalp EEG database, which is a component of the Physiobank database [31], are the two databases used in the study’s trials. The BCI competition dataset was actually used for binary classification problems. Each channel signal in the BCI III dataset 1 consists of 128 samples over 64 channels, and the dimension of a second test signal is $128\times 64$. In the Physiobank database, all used datasets consist of 23-channel EEG recordings, which were sampled at 256 samples per second with 16-bit resolution. The international 10–20 system of EEG electrode positions and nomenclature was used for these recordings. We firstly use the ‘chb01_31.edf’ signal from this database as our testing data to evaluate the performance.

As is the case with many other approaches [15,16,32], we use the metrics Mean Squared Error (MSE) and Matthews Correlation Coefficient (MCC) to assess the quality of the recovery results. The average squared deviation (MSE) between the true and estimated values quantifies the degree of dispersion in the estimation mistakes. Below is an outline of how MSE is expressed mathematically. $MSE={\displaystyle \sum _{g=1}^G}{\displaystyle \frac{{∥\stackrel{⌢}{X}-X∥}_F^2}{GNC}}$. $MCC={\displaystyle \sum _{g=1}^G}{\displaystyle \frac{vec{\left(X\right)}^{T}vec(\stackrel{⌢}{X})}{G{∥X∥}_F{∥\stackrel{⌢}{X}∥}_F}}$ measures the similarity of two waveforms G = 10 is the number of the experiments. Both X(X = $\frac{X}{{∥X∥}_F}$) and $\stackrel{⌢}{X}$ are normalized by their Frobenius norms, respectively.

The parameter setting of NLRC is as follows: the uniform random matrix as the sampling matrix; the number of compressive measurements is denoted by $rate=M/N$; $t=10$, $c=10$, $d=8$ and $k=30$; q and p are changed from 0.1 to 0.5 by step 0.1; ${\beta }_1=1$, ${\beta }_2=1$, ${\beta }_3=1$ on the first test data; ${\beta }_1=0.1$, ${\beta }_2=10$, ${\beta }_3=10$ on the second test data.

4.1. MSE and MCC Results on Test Data

We carefully adjusted the parameters of the various approaches to maximize the quality of the reconstructed signals in order to provide an impartial and equitable comparison between them. The Matthews Correlation Coefficient (MCC) and Mean Squared Error (MSE) values obtained using various methods at various sampling rates are shown in Figure 2, Figure 3, Figure 4 and Figure 5. It can be found that: (1) since SCLR-I and SCLR-A simultaneously exploits the cosparsity and low-rank property, they can obtain better results than BSBL; and (2) the NLRC approach uses block operations and low-rank approximation to great advantage, which makes it perform better than SCLR-I and SCLR-A. When NLRC uses these techniques, it may take advantage of the natural correlation that exists across various channels, which leads to better recovery outcomes than when using other methods. As a result, our NLRC methodology is unique in the reconstruction of multi-channel EEG signals, offering a more accurate and computationally efficient solution.

4.2. Influence of Variable p and q on Recovery Results

Moreover, we thoroughly assessed how various q and p variables affected the results of the reconstruction procedure.

Table 1. The optimal values of p and q on the first test data at different sensing rates.

Value	Rate
	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8
p	0.2	0.2	0.2	0.1	0.5	0.4	0.5	0.5
q	0.3	0.2	0.4	0.1	0.1	0.5	0.5	0.1

and

Table 2. The optimal values of p and q on the second test data at different sensing rates.

Value	Rate
	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8
p	0.3	0.1	0.1	0.5	0.1	0.2	0.4	0.1
q	0.3	0.2	0.2	0.5	0.3	0.3	0.1	0.2

list all the optimal values of q and p at different sensing rates. In Figure 6, we set $rate=0.1$ and $p=0.2$, and compare the MSE values under different q values on the first test data. In Figure 7, we set $rate=0.1$ and $q=0.3$, and compare the MSE values under different p values on the second test data. Figure 6 and Figure 7 show that NLRC with appropriate q and p can show better performance, while with inappropriate q and p is not. All in all, the selections of q and p are critical to NLRC. So, we take the same processing method as [33] to change the q and p values in our experiments in order to obtain the best performance.

4.3. Complexity Analysis

We estimate the time costs associated with our NLRC approach in order to evaluate model complexity. Using a Windows 10 operating system on a PC with an Intel Core i7-8550U CPU and 16 GB of RAM, the experiments are carried out on Matlab 2018a. As

Table 3. The time costing of different methods on the first test data at different sensing rates.

Method	Rate
	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8
BSBL	10.0	7.0	5.5	5.7	5.9	8.0	9.2	13.7
SCLR-I	10.5	12.0	14.1	16.2	19.0	22.0	25.0	26.0
SCLR-A	3.7	4.3	5.8	7.0	8.0	8.5	10.7	11.0
LQSP	3.7	4.3	5.8	7.0	8.0	8.6	10.7	11.0
NLRC	61.1	61.6	62.8	63.4	63.7	63.9	64.2	64.9

and

Table 4. The time costing of different methods on the second test data at different sensing rates.

Method	Rate
	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8
BSBL	8.3	7.7	7.6	7.7	7.3	10.0	11.4	17.0
SCLR-I	13.1	14.9	17.5	19.1	20.6	21.1	22.2	32.3
SCLR-A	4.6	5.3	7.2	8.7	10.0	10.6	13.3	13.7
LQSP	4.6	5.4	7.2	8.7	10.1	10.6	13.3	13.8
NLRC	76.7	77.9	78.3	76.8	74.9	77.2	76.4	74.4

illustrate, the NLRC technique needs a longer computing time than other approaches. This is mostly because every iteration requires an extra group matrix approximation. Nevertheless, NLRC justifies a notable gain in reconstruction quality in MSE and MCC.

5. Conclusions

This paper presents a new recovery method for multichannel EEG signals that combines cosparsity and nonlocal low-rank property into a single theoretical framework. Especially, block operations and nonlocal low-rank approximation methods are used to accurately extract interaction among different channels. Moreover, the analytical operator chooses second-order difference matrix to ensure cosparsity. To efficiently solve the resulting nonconvex optimization problem, we use the weighted Schatten-p norm and $l_q$ norm as proxy functions. Using the Alternating Direction Method of Multipliers (ADMM), the problem can be addressed. Based on the same measurement levels, experimental data show that the NLRC method performs much better than other competing reconstruction methodologies. However, the algorithm complexity of our method is slightly higher, and the next step of work we will consider improving the computational speed.