VSSI-GGD: A Variation Sparse EEG Source Imaging Approach Based on Generalized Gaussian Distribution

Electroencephalographic (EEG) source imaging (ESI) is a powerful method for studying brain functions and surgical resection of epileptic foci. However, accurately estimating the location and extent of brain sources remains challenging due to noise and background interference in EEG signals. To reconstruct extended brain sources, we propose a new ESI method called Variation Sparse Source Imaging based on Generalized Gaussian Distribution (VSSI-GGD). VSSI-GGD uses the generalized Gaussian prior as a sparse constraint on the spatial variation domain and embeds it into the Bayesian framework for source estimation. Using a variational technique, we approximate the intractable true posterior with a Gaussian density. Through convex analysis, the Bayesian inference problem is transformed entirely into a series of regularized ${L}_{{2}{p}}$ -norm ( ${0} < {p} < {1}$ ) optimization problems, which are efficiently solved with the ADMM algorithm. Imaging results of numerical simulations and human experimental dataset analysis reveal the superior performance of VSSI-GGD, which provides higher spatial resolution with clear boundaries compared to benchmark algorithms. VSSI-GGD can potentially serve as an effective and robust spatiotemporal EEG source imaging method. The source code of VSSI-GGD is available at https://github.com/Mashirops/VSSI-GGD.git.


I. INTRODUCTION
E LECTROENCEPHALOGRAPHY (EEG) is widely employed for brain imaging due to the millisecond-level time resolution [1], [2].However, EEG has limited spatial resolution because of the volume conduct effect [3].EEG source imaging (ESI) is widely recognized as an effective technique for improving the spatial resolution of scalp EEG signals.It is a technique that uses EEG signals and mathematical models of human head tissues to derive the image of intracranial neural activities, precisely reconstructing their locations and extents [3], [4].Previous studies on electrocorticography (ECoG) and functional magnetic resonance imaging (fMRI) have shown that brain activity exhibits spatial correlations, indicating that EEG signals may originate from widespread active regions of the cortex [5].In other words, it requires an activated brain area (usually no less than 6 cm 2 ) to generate a spike on the scalp [6].Reconstructing the spatial extent of cortical brain activity, along with the location, is not only possible but also holds significant scientific and medical importance [7], such as stroke recovery [8], [9], as well as detection and surgical resection of epileptic foci [10], [11].Additionally, ESI provides signals with a higher spatiotemporal resolution for brain-computer interfaces (BCI) and related applications [12], [13].
For ESI problems, the distributed current density (DCD) model is a widely employed source model [3], [14].The DCD model divides the cortex into several triangular grids, where each triangle represents a source with a fixed location.Typically, the number of sources is around several thousand and much larger than the number of scalp electrodes (usually 32 to 256 electrodes) [14], [15].Hence, ESI is a severely ill-posed inverse problem, where countless potential source configurations can fit the scalp measurements.It is similar to compressed sensing, as it relies on utilizing a limited quantity of observed data to reconstruct intricate signals.To obtain a unique and accurate solution, appropriate constraints need to be applied to narrow the solution space.
Conventional ESI methods employ the L 2 -norm regularize to constrain the potential sources.For example, the minimum norm estimation (MNE) applies L 2 -norm regularization to select the source configuration with the minimum energy that matches the scalp measurement [14].However, MNE is not sensitive to deep sources because the electric fields generated by superficial sources have less energy [3].To address this biased estimation, weighted MNE (wMNE) proposes to weight the regularization term with the norms of each column of the lead-field matrix [16].To obtain spatial coherence and smooth solutions, low-resolution electromagnetic tomography (LORETA) proposes to minimize the L 2 -norm of the second-order spatial derivative of the sources [17].Nonetheless, the imaging results of these L 2 -norm based methods are often blurry and cover most of the cortical regions [3], [18].To improve spatial resolution, many studies have employed sparse constraints such as the L 1 -norm further to narrow the solution space [19], [20].Other studies incorporate constraints as prior information into the Bayesian framework.One such method is sparse Bayesian learning (SBL) [21], [22].Under the SBL framework, the study in [23] combines the hyperparameters of the source variance along with the full-structure multivariate Gaussian noise for joint estimation, leading to promising outcomes.Cai et al. [24] utilize SBL to estimate the model data covariance and reconstruct potential sources with adaptive beamformers, effectively mitigating the effects of correlated brain sources and yielding good results.
The sparse constraint methods can accurately localize the activated brain sources with high spatial resolution.However, due to the sparsity enforcement on the source space, these methods usually generate only one or a few points within the active areas, thereby missing the extent information of source activities [14], [25].To remedy this issue, several studies applied sparse constraints in the transform domains, such as the variation domain [26], [27], [28], or a combination of multiple domains [29].These constraints can provide more accurate estimates of source location and extent.Nonetheless, most of these methods apply the L 1 -norm constraint in the spatial transform domains, which tends to overestimate the extent of the sources, especially for small ones [30].To reduce this bias, the study in [31] has pointed out that non-convex constraints like L p -norm (0 < p < 1) can acquire more accurate results in contrast to L 1 -norm-based methods, especially when measurements are noisy and restricted isometry property (RIP) is violated, which is always the case for ESI.Studies have indicated that the L p -norm is sparser than the L 1norm, as the inclusion of p allows for more flexible balancing between measurement data and sparsity [32].In addition, L p -norm-based methods require fewer measurements to yield reliable results [33].
In the Bayesian probabilistic framework, the solution of the typical L p -norm regularization methods is equivalent to the maximum a posterior (MAP) estimation with the generalized Gaussian prior.However, according to [34] and [35], this non-convex MAP reconstruction often encounters many shallow local optima.The posterior, which incorporates the generalized Gaussian prior (L p -norm constraint with 0 < p < 1), is asymmetric and exhibits multiple local modes.The MAP estimation is unable to account for the overall mass of the posterior and tends to become stuck in shallow optima.To address this issue, Bayesian integration can be employed to average out these shallow local optima and obtain the mean of the posterior distribution [35], [36].
Similar to the mixed L 21 -norm [25], [37], this work used the L 2 p -norm regularization to reconstruct the cortical extended sources.As a matrix form extension of L p -norm, we define this regularization operator to impose the L 2 -norm regularization in the temporal domain and the L p -norm regularization in the spatial variation domain.By formulating the ESI problem within the Bayesian framework, we introduced a new EEG source imaging algorithm called variation sparse source imaging based on generalized Gaussian distribution (VSSI-GGD).Using a variational technique, VSSI-GGD approximates the intractable posterior density with a Gaussian distribution.Through convex analysis, the approximated Gaussian posterior distribution is obtained by solving a series of regularized L 2 pnorm constraint problems.The alternating direction method of multipliers (ADMM) algorithm [38] is employed to solve each efficiently regularized L 2 p -norm constraint problem.By combining the generalized Gaussian prior and variational Bayesian inference technique, we anticipate that the proposed method will yield more accurate reconstructions of cortical sources compared to existing techniques.This is validated by simulations using synthetic and experimental EEG data.
The rest of this paper is summarized as follows.Section II provides a detailed description of VSSI-GGD.Section III presents the simulation design, benchmark algorithms, and performance metrics used to validate the performance of our proposed method.Section IV presents the results obtained using both simulated and real EEG data, followed by a discussion in Section V and a brief conclusion in Section VI.
Matrices and vectors in this paper are denoted using uppercase and lowercase boldface letters, respectively.N (µ, ) represents a Gaussian distribution with mean µ and covariance .f * and f −1 are the conjugate and inverse function of f , respectively.x i and x j,: denote the ith column and the jth row of matrix X, respectively.

A. Background
According to the quasi-static approximation of Maxwell's equations, the relationship of EEG signals and current sources is described as [3], [14], where B ∈ R N b ×N t is the EEG data collected from N b electrodes at N t time sampling points, and S ∈ R N s ×N t is the cortical source of N s potential brain sources at N t time points.L ∈ R N b ×N s is the lead-field matrix, which models how potential brain activities, S, are related to scalp EEG recordings.ε denotes the measurement noise which is assumed to be drawn from N (0, ε ) independently, and ε is the covariance obtained from pre-stimulus data.For generality, spatial whitening is applied to both the lead-field matrix and EEG measurements based on the estimated ε , resulting in ε ∼ N (0, I) [25].
Based on (1), ESI aims to estimate the appropriate potential source configuration S from the EEG measurement data B with a given lead-field matrix L. However, the number of suitable source configurations tends to approach infinity because the number of EEG electrodes is much less than the number of potential sources, i.e., N b ≪ N s .Thus, prior spatiotemporal constraints are needed to narrow the solution space for the illposed problem.For the temporal constraint, we assume that the source matrix S is a linear combination of K temporal basis functions (TBFs) [15], [25], is the kth TBF.By projecting EEG measurements B and sources S onto the subspace spanned by the TBFs, we can capture the temporal properties of brain sources through the projection coefficients B = B ⊤ ∈ R N b ×K , and S = S ⊤ ∈ R N s ×K .Hence, the data model in (1) can be written as where Here, we construct TBFs from singular value decomposition (SVD) of EEG signals and use the Kaiser criterion to determine the value of K [25].For orthogonal TBFs, we obtain εk ∼ N (0, I) [15].After obtaining S, we recover the potential source S as S = S .The use of TBFs not only reduces the dimension of the ESI problem but also helps stabilize the source reconstruction.For notation simplicity, we use S and B instead of S and B throughout the remainder of this paper.
For spatial constraint, to obtain locally smooth and globally clustered cortical sources [5], [15], we enforce source sparsity on the spatial variation domain and penalize the amplitude differences between adjacent dipoles.Specifically, we apply a linear transformation to the source through the variation operator V to obtain the variation source U. V is defined as [26], [39] where M is the number of triangular edges in the source model.Each row of V corresponds to two adjacent dipoles that share the same edge in the triangular grid on the cerebral cortex, indicated by the 1 and -1 values.Hence, each row of variation sources U = V S ∈ R M×K represents the amplitude difference between two adjacent sources.The use of the variation transformation enforces sparsity on the variation source, preventing the loss of spatial information caused by enforcing sparsity in the original domain [7].

B. Inverse Solution
To reconstruct extended sources, we assume that the variation source U is sparse, and each row of U is independent of each other.We achieve this by incorporating the generalized Gaussian prior for each row of U [40] p(u i,: where λ and p are the scale and shape parameters of the generalized Gaussian distribution, and (•) is the Gamma function.Then the prior of U can be expressed as Based on the data model (1) and prior ( 5), the posterior distribution of S is derived as where Z is the model evidence.However, due the the non-Gaussian prior, the integration of Z is intractable.To efficiently leverage information from the posterior distribution, we utilize a variational technique to find a Gaussian density Q(S|B) to fit the true posterior p(S|B).
According to the convex analysis, the prior p(u i,: ) can be approximated as [41] p(u i,: ) = max where + C, and C is a constant.Therefore, we have in which u t denotes the tth column of U and = diag(γ ), where γ = [γ 1 , . . ., γ M ] ⊤ is a vector comprising M unknown nonnegative hyperparameters that govern the prior covari- With this Gaussian-form lower bound, the approximated Gaussian posterior are obtained as Given a value of γ , the mean µ and covariance of Q(S|B; γ ) are represented as follows As a result, we aim to determine the optimal value of γ that makes Gaussian density Q(S|B; γ ) the tightest to p(S|B).A natural way to fit p(S|B) by Q(S|B; γ ) is to maximize the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
The inference problem falls to maximize the right hand of Eq. ( 11) w.r.t.S and γ .Combining (8) and taking the negative logarithm of (11), we have For computationally convenience, we construct auxiliary functions using convex analysis to form an alternative objective function for log | A|.Since it is concave in γ , by Fenchel duality, it can be expressed as a minimum over upper-bounding hyperplanes [21], [39] log where ⊤ as auxiliary variables [21], [25], [39], which is obtained from simple geometric considerations as the slope at the current γ of log | A| In this work, we utilize the Lanczos algorithm to provide a scalable low-rank approximation of the covariance A −1 [42].
Combining (12) with ( 13) and removing the minimizations, the modified cost function is derived as The hyperparameter γ is obtained by minimizing the cost function ψ(γ , S, ξ ) with other variables fixed Taking the derivation of γ i , we have Note that g(x) is concave, so the inverse function of its concave conjugate is itself [41].Thus Considering the terms related to S in ψ(γ , S, ξ ), the source S is updated as Since we have, Hence, S is obtained by solving the following regularized optimization problem Eq. ( 22) can be fully solved by ADMM algorithm w.r.t. the augmented Lagrangian function of ( 22) where Z ∈ R M×K is the Lagrangian multiplier and ρ > 0 is penalty parameter.S, U, and Z are updated alternately as where S k , U k , and Z k are the corresponding output of the kth ADMM iteration, and Here, the optimization problem of U k+1 is solved using iteratively reweighted least squares (IRLS) algorithm [43] where Here, the scale parameter λ is chosen using cross-validation.The ADMM penalty parameter ρ is selected using adaptive residual balancing strategy [44], and the update for ρ is performed after every 10 ADMM iterations.To determine the value of p, according to prior studies [31], [32], [33]  and our numerical results provided in the supplementary materials, we experimentally set p = 0.8.The proposed algorithm was conducted on a standard PC (Corei9-10980XE CPU 3 GHz and 128GB RAM).During our numerical simulation, VSSI-GGD converges after ten cycles, where each cycle consists of 400 ADMM iterations which takes about 30 seconds.A schema of VSSI-GGD above is presented in Algorithm 1, and the code of it is offered at https://github.com/Mashirops/VSSI-GGD.git.update S k , U k , Z k with (24).

III. SIMULATION PROTOCOL AND PERFORMANCE METRICS
This section provides a detailed description of the benchmark algorithms, Monte Carlo numerical simulations, and performance metrics used to validate the performance of the proposed method.

B. Simulation Design
We conducted a series of Monte Carlo numerical simulations to verify the performance of these ESI algorithms due to the lack of ground truth.In this work, we utilized Brainstorm with the default ICBM 152 magnetic resonance imaging (MRI) to generate the head model [45].The cortex surface was divided into 6006 triangular grids, with each triangle representing a potential brain source.The lead-field matrix was calculated based on the 64-channel Neuroscan Quik-cap sensor system and the OpenMEEG [46] (with two non-EEG channels abandoned).To construct a simulated cortical source, a seed point was randomly selected on the cortical triangular grids, and adjacent triangles were added to the cluster until the source region reached a specified value.The time series of the simulated ground truth, denoted as S r eal , was divided into pre-and post-stimulus parts.The post-stimulus parts were simulated by applying several Gaussian-damped sinusoidal signals to the activated source region.Interference noise, labeled as S noise , was added to mimic background source activities.Specifically, five clusters outside the real sources were randomly chosen as interference sources with an extent of 3-5 cm 2 .Pink noise with 1/ f power was applied to these interference sources.The final simulated source signals S were obtained by the weighted sum of normalized S r eal and S noise under a specific signal-to-noise/interference ratio (SNIR), Here, ϕ is the ratio of real brain source signals to interference noise, and SNIR is defined as 10 lg 1 ϕ [15].Gaussian white noise was added as the measurement noise under a specified SNR defined as SNR = 10 lg ∥L S∥ F ∥ε∥ F .Consequently, the simulated observed EEG data, denoted as B, was obtained through the forward model.
The performance of ESI methods was compared under various scenarios 1) Various extents -we tested the influence of sources with different extents on ESI.Specifically, we considered source extents of 2 cm 2 , 5 cm 2 , 10 cm 2 , 18 cm 2 , and 32 cm 2 for a single source; Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
-similarly, we evaluated the performance of the proposed methods under different measurement noise levels (0, 3, 5, and 10 dB) with a single patch source of 8 cm 2 ; 3) Various numbers of channels -we investigated the effect of varying amounts of data on the algorithm results by considering different numbers of channels with SNR and SNIR set to 5 dB.Specifically, we down-sampled the EEG channels and analyzed cases with 100% (62 channels), 75% (46 channels), 50% (32 channels), and 25% (16 channels) of the data.In each simulation, channels were manually selected for non-complete cases (refer to the supplement document for detailed information).In the aforementioned scenarios, unless specified otherwise, the SNR and SNIR were set to 5 dB with a single source of 8 cm 2 .In the supplement document, we also present the results of different ESI methods with various SNIR levels.Each case involved 50 Monte Carlo simulations.

C. Performance Metrics
Four performance metrics are employed to evaluate the performance of ESI algorithms.
1) The area under the receiver operating characteristic (ROC) curve (AUC) [4], [47] evaluates the detection accuracy and sensitivity of the reconstructed sources.2) Spatial dispersion (SD) [29], [48] measures the spatial blurring of the estimated sources compared with the ground truth.
3) The distance of localization error (DLE) [29], [48] measures the localization error of the reconstructed sources.4) The shape error (SE) [25], [39] quantifies the relative squared error between the normalized reconstructed and simulated sources.For the details of each performance metric, please refer to [4] and [39].In general, higher AUC values and lower SD, DLE, and SE values indicate better performance of the ESI methods.To assess the statistical significance, we utilized the Kruskal-Wallis test.If the statistic from this test is significant, we then conducted Wilcoxon rank sum tests, which are adjusted by the false discovery rate (FDR), within VSSI-GGD against each benchmark algorithm.Moreover, we employed Otsu's threshold to visualize the imaging results as in [25].

A. Results of Monte Carlo Simulations 1) Influence of Various Extents:
We first evaluated the influence of different source extents (2, 5, 10, 18, and 32 cm 2 ) for ESI methods with a single patch.The performance metrics of the results are displayed in Fig. 1.For very small size sources (i.e., 2 cm 2 ), all methods except SBL show poor performance with lower AUC and higher SD, DLE, and SE values.SBL demonstrates exceptional performance, indicated by nearly one AUC and zero SD, DLE, and SE values.This is attributed to SBL over-enforcing sparsity, which always produces point estimations.As the source extents increase, the SD and DLE values of SBL consistently remain low, but the AUC values significantly decrease.On the contrary, LORETA often generates diffuse estimations, as evidenced by the high SD values.For source extents greater than 18 cm 2 , LORETA shows the highest AUC values, indicating it is highly proficient at detecting large-scale sources, which aligns with the findings of [39] and [47].In most cases, VSSI-GGD achieves larger AUC values and smaller SD (except that of SBL), DLE, and SE values than the compared methods.The results demonstrate the superior performance of VSSI-GGD in reconstructing the location and extent of extended sources.
One example of source imaging of different extents is presented in Fig. 2, in which all imaging results are thresholded by Otsu's method.The results of both wMNE and LORETA are diffused across all extents.While the results of SBL always shrink to several points, in the case of extended cortical activities, SBL loses most of the extent information.Four VSSI-based methods provide appropriate extents estimation, while VSSI-L 2 p , VSSI-CM, and VSSI-MAP generate some spurious sources around the ground truth.Overall, the estimated results from VSSI-GGD are closer to the simulated sources in terms of localization and extent among seven ESI algorithms, which is also supported by the performance metrics provided in the supplementary document.
2) Influence of Various SNRs: We next assessed the performance of each ESI method against different SNR levels (0 dB, 3 dB, 5 dB, 10 dB), and the performance metrics are shown in Fig. 3.As SNR increases, all algorithms show improved performance, as evidenced by increased AUC and decreased SD, DLE, and SE values.Since SBL always produces focal estimations, the SD and SE values of SBL are not sensitive to noise.However, as found in the previous scenarios, SBL provides little information on source extents.VSSI-GGD outperforms the other benchmark methods under various measurement noise levels, indicated by higher AUC and lower SD and SE values.

3) Influence of Various Number of Channels:
To evaluate the influence of data volume, we down-sampled the 62 EEG channels to 46, 32, and 16 channels.In the supplement document, we provide the electrode distribution map in 4 cases.Most methods show a declined performance as the number of channels decreases, indicated by the decreased AUC and increased SD, DLE, and SE values.However, the performance metrics of SBL remain almost unchanged.This stability can be attributed to its focal estimations.Except for the SD values of SBL, VSSI-GGD achieves larger AUC, and lower SD, DLE, and SE values than the compared methods in all conditions, suggesting the superior performance and robustness of VSSI-GGD in the presence of missing data.

B. Results of Experimental EEG Data Analysis
We further apply these ESI algorithms on two real human EEG datasets to assess the practical efficacy of each method.
1) Analysis of Localize-MI Dataset: The Localize-MI dataset is an open dataset that includes high-density EEG data recorded from 61 sessions of 7 subjects with drug-resistant epilepsy to evaluate source location methods [49], which is available at.This dataset includes the spatial locations of stimulating contacts within the brain.These contacts were used to precisely deliver electrical currents during the presurgical evaluation of patients with drug-resistant epilepsy.Hence, Localize-MI provides the stimulation site as the ground truth.The EEG data were recorded at a sampling rate of 8000 Hz across 256 channels.We re-referenced the data to the average of all good channels and clipped the time window from −2 to 2 ms.The dataset also provides the surface model with 8196 sources and a corresponding lead-field matrix.We applied those ESI methods to the averaged EEG data of each session to reconstruct the potential sources and compare them with the ground truth.The MNE is used to visualize the results [50].
The source reconstruction results for the average time course of session 1 of subject 01 are depicted in Fig. 5. Imaging results of other subjects are provided in the supplementary material due to space limitations.In each imaging result, the green point indicates the ground truth.As shown in Fig. 5, the reconstructed source of VSSI-MAP misses the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.ground truth point, whereas SBL estimates some point sources around the ground truth.Among the other six methods, VSSI-GGD provides more compact estimated sources with clear boundaries.To quantitatively assess the performance of each algorithm, we computed the spatial dispersion defined in [49] for all 61 sessions, and results are presented in Table II.VSSI-GGD provides the lowest SD values across all subjects in line with the imaging results.

TABLE II SPATIAL DISPERSION (MM) OF SEVEN ESI ALGORITHMS FOR ALL SESSIONS OF EACH SUBJECT (MEAN±STD). PERFORMANCE STATISTICAL SIGNIFICANCE LEVEL IS DENOTED
2) Analysis of Public Epilepsy Dataset: We further applied the proposed method to the public epilepsy EEG dataset that belongs to the Brainstorm tutorial.The dataset is available at https://neuroimage.usc.edu/brainstorm/DatasetEpilepsy.The EEG data was recorded from a patient suffering from focal epilepsy and various types of seizures.Besides, this dataset was fully analyzed and reported in [27] and [51].Here, we follow the Brainstorm tutorial to derive the head model, lead-field matrix, and averaged EEG data for ESI analysis.Source estimation results at the peak of averaged EEG data (0 ms) are presented in Fig. 6.All methods locate the source on the left frontal lobe, while wMNE and LORETA are widely spread, covering a large area.In contrast, SBL provides a focal estimation that consists of only one point source.Results of four VSSI-based methods are highly similar and consistent with the clinical findings in [51] and the results in [27].Among all, VSSI-GGD provides the clearest boundary between active and non-active sources.

V. DISCUSSIONS
We proposed an EEG extended source imaging algorithm, VSSI-GGD, within the Bayesian framework.VSSI-GGD employs the generalized Gaussian prior in the variation domain to enforce the sparsity of brain sources.Using variational inference, we obtained the approximated posterior by solving a series of regularized L 2 p -norm regularization problems, which can be efficiently solved using the ADMM method.Both Monte Carlo numerical simulations and real EEG data analysis have revealed the superior performance of VSSI-GGD in estimating the location and extent of brain sources.
ESI plays a crucial role in stroke rehabilitation and resection of epileptic foci, as previous studies have confirmed the capability of EEG to capture potential neural activities [1], [7].However, it is limited by its low spatial resolution due to the volume conduction effect.Furthermore, the number of potential sources is much larger than the number of EEG electrodes, making this problem inherently ill-posed.To improve the spatial resolution, appropriate constraints are necessary to narrow the source space.
The L 2 -norm constraint is a widely used method in source estimation techniques such as wMNE and LORETA due to its computational convenience [3], [14], [30].However, one drawback of this method is its tendency to produce diffuse estimations, as shown in Fig. 2, where the results of wMNE and LORETA cover a large portion of the cortical areas with low spatial resolution.Ideally, an activated potential brain source is expected to be a patch of moderate size with clear boundaries.Unfortunately, the information on source location and extent provided by wMNE and LORETA is often too blurred.To improve spatial resolution, some studies have employed sparse constraint methods such as L 1 -norm and sparse Bayesian learning [19], [22], [24].However, as depicted in Fig. 2, the results of these methods often collapse to a single point within the ground truth due to the direct imposition of sparse constraints in the original domain [7], [25], [39].
Several studies have shown that potential brain sources are sparse in certain transform domains, such as the variation domain, rather than the original domain [19], [25], [27], [39].Therefore, enforcing sparsity in the original domain can lead to information loss.Keeping these factors in mind, in this work, we project the sources onto the variation domain to obtain sparse variation sources, following the strategy as prior studies [25], [26], [39].Previous studies have demonstrated that variation sparsity-based methods are effective in providing clear information about source location and extent [7], [25], [26], [27], [39].For VSSI-GGD, we assume that each row of the variation source follows the generalized Gaussian distribution.This statistical model more flexible and suitable for non-Gaussian distribution variables like EEG signals, compared to the common Gaussian model [40].The source sparsity is ensured by limiting the shape parameter p range from 0 to 1.The generalized Gaussian prior is then embedded into the Bayesian framework for statistical inference.It's important to note that our objective is not to directly maximize the posterior distribution, but rather to obtain the posterior's mean by maximizing the model evidence's lower bound.This is because the MAP estimation only gives the mode of p(S|B) and blind to the mass, especially for the multimodal asymmetric posterior distribution [39], The Bayesian framework provides a flexible framework where many ESI algorithms can be derived and explained.
There is a direct connection between many classical distributions and norms within this framework.For instance, when the prior is assumed to be a Gaussian distribution, the L 2 -norm based method MNE corresponds to the MAP estimation of this prior [21].Compared to regularization methods, the Bayesian framework is powerful because it allows for explicit modeling of prior assumptions using probability distributions and hyper-parameters in a fully data-driven manner.In the case of VSSI-GGD, the generalized Gaussian prior we use corresponds to the L p -norm (0 < p < 1).The L p -norm is also a sparse constraint similar to the L 1 -norm.Both norms are approximations of the NP-hard L 0 -norm.However, as suggested by [30], methods based on the L 1 -norm on the variation domain tend to overestimate the extent of sources, particularly for small-sized sources.On the other hand, the more sparse mathematical nature of the L p -norm allows it to better fit the measurement and source signals.This is also confirmed in the study of [31], where they compared these two sparse constraint-based methods and revealed that L p -norm can mitigate the bias caused by L 1 -norm.
VSSI-GGD is novel as we transform the generalized Gaussian prior source estimation problem into a regularized L 2 p -norm optimization problem through convex analysis.The L 2 p -norm combines the L 2 -norm in the temporal domain and the L p -norm in the spatial domain, similar to the L 21norm.This allows for effective localization of active sources and improved estimation of their spatial extent compared to traditional single-norm methods.However, the L 21 -norm also suffers from the mathematical properties of the L 1 -norm we mentioned above.To overcome this limitation, the L 2 pnorm employs the L p -norm to constrain the spatial domain of the potential sources for better spatial sparsity.VSSI-GGD outperforms VSSI-CM in reconstructing extended sources, as suggested by the higher AUC values and lower SD, DLE, and SE values achieved by VSSI-GGD shown in Fig. 1.Similarly, this conclusion holds for VSSI-L 2 p and VSSI-MAP.
In general, the proposed VSSI-GGD narrows the solution space through temporal and spatial constraints in order to obtain a unique solution like most state-of-the-art methods.Specifically, VSSI-GGD utilizes a better sparse generalized Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Gaussian inference within the Bayesian framework, transforming intractable posterior estimation problem into a regularized optimization problem.Besides, efficient solutions can be achieved using existing methods such as ADMM and IRLS.The results of simulation experiments and real data analysis show that it is precisely because of the sparser constraint and Bayesian variation inference that, for extended source reconstruction, VSSI-GGD outperforms other ESI methods in terms of accuracy in localization and precision in extent estimation.Nevertheless, there are certain areas where improvements can be made to this method.Firstly, similar to other optimization methods, the computation time of VSSI-GGD is not short.It takes several minutes to complete the source estimation using our method, which does not meet the of real-time systems, where millisecond-level response is expected.Secondly, determining an appropriate value of the shape parameter (referred to as the value of p) of the generalized Gaussian prior is crucial for accurately fitting the measurements in our proposed ESI algorithm.Based on our experiments (result is provided in the supplementary material) and previous studies [31], [32], [33], we finally set p = 0.8.An advanced approach to handle this problem is to utilize the algorithm unrolling theory [52] to build a deep neural network (DNN) using the iteration framework of VSSI-GGD.Algorithm unrolling is a powerful tool to connect existing iterative algorithms with DNN, where all parameters in the original algorithms can be learned through end-to-end training.Our future work will combine the Bayesian update rules and the algorithm unrolling technique to derive a full data-driven interpretable neural network for ESI.Hence, both the value of p and λ can be determined in a fully data-driven way.Moreover, by changing the iterative methods into DNN form, we can derive the output at the millisecond level, which is crucial in real-time source imaging applications.

VI. CONCLUSION
This paper presents a novel algorithm for EEG extended source imaging within the Bayesian framework.Our method exhibits robustness against noise and interference from numerical simulation perspectives, delivering high spatial resolution results for reconstructing extended sources, which outperforms other benchmark algorithms.Moreover, the results of the human epilepsy dataset also suggest the potential of our method in real human EEG processing.In addition, in future work, our method can be further elucidated in neural networks through algorithm unrolling techniques to generate fully data-driven interpretable neural networks.

Fig. 1 .
Fig. 1.Performance metrics of various extents.This figure presents the Mean ± SEM of the results for 50 Monte-Carlo simulations.* and ** indicate that the corresponding performance metric of VSSI-GGD is significantly different from those of the compared methods with p correct < 0.05 and p correct < 0.01 respectively.

Fig. 2 .
Fig. 2. One imaging example of the simulated and estimated sources by ESI algorithms under various extents, and the performance metrics of the results corresponding to this instance are listed below.The thresholds of each estimated source are determined by Otsu's method.

Fig. 3 .
Fig. 3. Performance metrics of various SNRs.This figure shows the Mean ± SEM of the results for 50 Monte-Carlo simulations.* and ** indicate that the corresponding performance metric of VSSI-GGD is significantly different from those of the compared methods with p correct < 0.05 and p correct < 0.01 respectively.

Fig. 4 .
Fig. 4. Performance metrics of various number of channels.This figure shows the Mean ± SEM of the results for 50 Monte-Carlo simulations.* and ** indicate that the corresponding performance metric of VSSI-GGD is significantly different from those of the compared methods with p correct < 0.05 and p correct < 0.01 respectively.

Fig. 5 .
Fig. 5. Source reconstruction results for the first session of subject 01 from the Localize-MI dataset.The green point is the sEEG simulation site.

TABLE I LIST
OF SYMBOLS IN THIS PAPER