Solution of the EEG inverse problem by random dipole sampling

Electroencephalography (EEG) source imaging aims to reconstruct brain activity maps from the neuroelectric potential difference measured on the skull. To obtain the brain activity map, we need to solve an ill-posed and ill-conditioned inverse problem that requires regularization techniques to make the solution viable. When dealing with real-time applications, dimensionality reduction techniques can be used to reduce the computational load required to evaluate the numerical solution of the EEG inverse problem. To this end, in this paper we use the random dipole sampling method, in which a Monte Carlo technique is used to reduce the number of neural sources. This is equivalent to reducing the number of the unknowns in the inverse problem and can be seen as a first regularization step. Then, we solve the reduced EEG inverse problem with two popular inversion methods, the weighted Minimum Norm Estimate (wMNE) and the standardized LOw Resolution brain Electromagnetic TomogrAphy (sLORETA). The main result of this paper is the error estimates of the reconstructed activity map obtained with the randomized version of wMNE and sLORETA. Numerical experiments on synthetic EEG data demonstrate the effectiveness of the random dipole sampling method.

Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
The electroencephalography (EEG) inverse problem aims to reconstruct the neuroelectric activity map from the electric potential differences measured by a set of electrodes located on the skull.Since EEG is a non-invasive technique and EEG devices are quite cheap, EEG source imaging is widely used both in clinical applications and neuroscience studies.Furthermore, EEG measurements come to a high temporal resolution-on the order of milliseconds-making this brain imaging technique of primary interest in real-time applications.For details we refer the reader to [1,2] and references therein.
What makes EEG source imaging tricky is the fact that to obtain the brain activity map, we have to solve a severely ill-posed inverse problem and regularization techniques are needed to make the solution viable [1,2].Unfortunately, the numerical solution of the EEG inverse problem can be time and memory consuming, which makes it unattractive in real-time applications, especially when portable devices are employed.To reduce the computational cost, dimensionality reduction techniques can be used to extract the most significant features from the data or to select the brain regions of interest before applying an inversion method [3][4][5].In addition, dimensionality reduction can act as a first regularization step, since it reduces the degrees of freedom of the inverse problem.
In this paper, we consider a slim dimensionality reduction technique-the random dipole sampling method-where the dimensionality reduction is achieved by using a Monte Carlo technique to reduce the number of the unknowns-neural sources-without imposing any a priori information.This method was introduced in [6] to solve the magnetoencephalography (MEG) inverse problem, in which neuroelectric activity maps are reconstructed from the neuromagnetic field measured outside the head [1,7].Tests on synthetic and real magnetic data have shown that the method has a low computational load while keeping good accuracy [8,9].Here, we use the random dipole sampling method for reducing the dimensionality of the EEG inverse problem and show its effectiveness when using weighted Minimum Norm Estimate (wMNE) [10] and standardized LOw Resolution brain Electromagnetic TomogrAphy (sLORETA) [11] to reconstruct the brain activity map.The wMNE and sLORETA algorithms are two of the most popular inversion methods for EEG imaging and are also easy to implement.For this reason, they are particularly suitable for real-time applications.The main result of this paper is the error estimate for the reconstructed activity map obtained by the randomized version of wMNE and sLORETA.In the previous papers [6,8,9], the effectiveness of the method was shown through numerical tests but no error estimates were provided.The estimates we provide give an upper bound for the error and allow us to identify the optimal number of dipoles that need to be sampled to achieve good accuracy.Numerical experiments conducted on synthetic EEG data confirm our findings.
In literature, random methods are commonly used for dimensionality reduction of regression problems (see, for instance [12][13][14][15]), while their use for the solution of inverse problems is still limited and mainly based on the randomized singular value decomposition (randSVD) described in [16].In [17][18][19][20][21] randSVD was used to solve inverse problems by classical regularization methods, i.e.Tikhonov regularization and truncated SVD.A different approach was proposed in [22], where the authors introduced a sampled iterative method that converges to the Tikhonov solution.A randomized method for MEG/EEG imaging was introduced in [23] with the aim of simultaneously detecting cortical and subcortical brain activity [24].The method decomposes the source space into multiple randomized multiresolution levels, then solves the MEG/EEG inverse problem by a hierarchical Bayesian model implying Mimimum Norm Estimate (MNE) as a particular case.
The random dipole sampling method we consider was introduce in [6] and is related in some way to the random column sampling method described in [25], but its motivation stems from the neurophysiological assumption that brain activity is spatially sparse [26] and therefore only a small number of neural sources are sufficient to represent the neuroelectric current flowing in the brain.The key ingredient of the method is that the neural sources are chosen randomly by drawing from the uniform probability distribution, making the method easy to implement and efficient when dealing with real-time applications.Here, we use the method along with wMNE and sLORETA to solve the EEG inverse problem.
The paper is organized as follows.In section 2 some preliminaries on linear algebra and the EEG inverse problem are given as well as the description of the random dipole sampling method.The main results of the paper are given in section 3, where the error estimates are obtained.The results of the numerical experiments we conducted are displayed in section 4, while a discussion of our findings is given in section 5. Finally, in section 6 we draw some conclusions and outline future developments.

Materials and methods
In this section, after recalling some basic facts from linear algebra (see section 2.1), we describe the EEG inverse problem and the inversion methods we are interested in (see section 2.2).The random dipole sampling method is described in section 2.3.

Linear algebra preliminaries
In this section we recall same definitions and properties that we will use in the following.For details, we refer to the classical books on linear algebra [27,28] and the review paper [15].
R m×n be a real matrix.The symbol ∥v∥ denotes the Euclidean norm of vectors while ∥A∥ denotes the usual spectral norm of matrices.We write ∥A∥ F for the Frobenius norm.We recall that where are the singular values of A and r = rank(A).It follows that Both norms are unitarily invariant, i.e. for any matrix where U ∈ R m×m and V ∈ R n×n are orthonormal matrices.
The condition number of a matrix A is We will make use of the following inequalities (see [27]).For any rectangular matrices A, δA ∈ R m×n it holds (2.6) In the following we will also make use of the following property of norms.For any pair of vectors v and u and matrices A and B it holds (2.7)

The EEG inverse problem
The EEG inverse problem aims to reconstruct the neuroelectric current flowing inside the brain once the electric difference potential on the scalp, generated by the activity of one or more neuroelectric sources, is given.The cortical surfaces are extracted from MRI images and discretized in a regular triangulation whose nodes form the source space.A current dipole is located in each point of the source space.The unknown is the current dipole moment vector Q = [q 1 , . . ., q n ], where q k = [q k x , q k y , q k z ] is the dipole moment of the kth dipole.We assume that the measurement vector b ∈ R m and the current dipole moment vector Q ∈ R 3n are linearly related, i.e.
where M = [M(r 1 ), . . ., M(r n )] ∈ R m×3n is the lead field matrix.The sub-matrix M(r k ) ∈ R m×3 , 1 ⩽ k ⩽ n, represents the electric potential produced by a unit current dipole located in r k .For details see [2,7].
It is well known that the EEG inverse problem is an underdetermined ill-posed and illconditioned inverse problem [1,2] for which regularization techniques are needed to obtain a feasible solution.Two of the most widely used inversion methods for solving the EEG inverse problem are the wMNE [29] and sLORETA [11] algorithms.
The estimated current distribution obtained by wMNE is where R ∈ R 3n×3n is the source covariance matrix and Σ ∈ R m×m is the noise covariance matrix.R is a diagonal matrix whose diagonal entries are where ρ is the depth weighting parameter.For the EEG inverse problem suitable values of ρ are in the interval (0, 5], depending on the source density, the dipole orientation and the SNR [29].Without loss of generality, we can assume the data has been whitened so that Σ = I.Thus, defining the inversion matrix the wMNE solution can be written as (2.12) Starting from the estimate Q wMNE , sLORETA produces a statistical measure of the brain activity normalizing the wMNE estimate w.r.t. the variance of the estimated sources.For any dipole, sLORETA produces a neural activity index z k as follows, where q k wMNE is the estimated activity of the dipole located in r k and Σ k QwMNE ∈ R 3×3 is the kth diagonal block of the variance of the estimated sources Σ QwMNE ∈ R 3n×3n .It can be shown that usually known as the wMNE resolution matrix.
When the cortical surface is known, the dipoles have fixed orientations since they are perpendicular to the cortical surface.In this case, the unknowns are the intensity of the n sources, i.e.Q = [∥q 1 ∥, . . ., ∥q n ∥], and the orientations can be plugged into the lead field M, which now is an m × n matrix.Similarly to (2.12), the wMNE solution Q wMNE = [∥q 1 wMNE ∥, . . ., ∥q n wMNE ∥] is given by [30] where and R ∈ R n×n has diagonal entries being M(r k ) a vector.In case of fixed orientation, the sLORETA index simplifies as where Σ k QwMNE is the kth entry of the n × n source covariance matrix Σ QwMNE .It can be shown [31] that (2.18) is equivalent to

Random dipole sampling method
We can use the random sampling to select a small subset of source points so reducing the number of columns of the lead field matrix.Let I c = (k 1 , . . ., k c ) denote a subset of c indices of the dipoles drawn from the uniform probability distribution with replacement.We will denote by ν k , k ∈ I c , with c ⩾ ν k ⩾ 1, the number of times the kth dipole was drawn.The sampled lead field matrix is where the entries of the sampling matrix S ∈ R 3n×3c are all zeros except for the 3 × 3 blocks indicating that the kth source point was drawn at the jth drawing.Assuming c ≪ n, the sampled lead field matrix M S has far fewer columns than the full matrix M. The factor n c ensures that the singular values of M S are an unbiased estimate of the singular values of M [15].
We can obtain an approximate solution of the EEG inverse problem by using M S in either wMNE or sLORETA.Thus, the estimated randomized neuroelectric current obtained using wMNE as inverse solver, is given by where R S ∈ R 3c×3c can be obtained from R retaining just the entries corresponding to the drawn dipoles.Since R is a diagonal matrix, it follows that Similarly, the randomized sLORETA statistical measure is where Σk

QwMNE
is the variance of the estimated sources QwMNE .In case of fixed orientation, the randomized wMNE method estimates the intensity vector QwMNE = l[∥q k1 wMNE ∥, . . ., ∥q kc wMNE ∥ so that the sampled lead field matrix is where the only entries of the sampling matrix S ∈ R n×c different from zeros are the entries (k, j), k ∈ I c , 1 ⩽ j ⩽ c, which are equal to 1.
In this case, the randomized wMNE solution is given by with while the randomized sLORETA index becomes

Error estimates
In this section we define the local reconstruction error for the randomized version of both wMNE and sLORETA and obtain deterministic estimates for these errors.We note that in order for the estimates of theorems 3.3 and 3.6 to hold with high probability, we can repeat the random dipole sampling algorithm a suitable number of times [16].
For the sake of simplicity, in the following we consider the case of fixed orientation.The results can be easily extended to the general case (see section 4.2).

Randomized wMNE
For the randomized wMNE method we define the local reconstruction error as We note that In the following we will denote by S the matrix that selects the dipoles that have not been drawn, i.Let G wMNE and G wMNE be the full and the sampled inversion matrices, respectively.For the minimum singular value of G wMNE and G wMNE the following bounds hold ) The proof is in the appendix.
Proposition 3.2.Let Q wMNE and Q wMNE be the neuroelectric current estimates obtained with the exact wMNE and the randomized wMNE, respectively.The following estimates hold where ν = max 1⩽k⩽c ν k ⩾ 1, and The proof is in the appendix.
Theorem 3.3.For the reconstruction error E wMNE the following estimate holds ) Thus, we obtain the following bound for the relative error, Moreover, ∥ ST Q wMNE ∥ → 0 when c → n.Thus, as expected, in this case the reconstruction error E wMNE goes to 0 when c → n.
Remark 2. In case of Tikhonov regularization the solution can be written as where α > 0 is the regularization parameter.Q Tik is also known as minimum norm estimate.Its random version is Thus, the error estimate in theorem 3.3 can be easily adapted to the randomized Tikhonov method by substituting in (3.5) µ with α and κ(G wMNE ) with κ(M T M + αI).
Furthermore, we can obtain the error estimate for the randomized least squares solution by substituting µ with σ min (M T M) and κ(G wMNE ) with κ(M T M).

Randomized sLORETA
In this section we give an upper bound for the reconstruction error when using randomized sLORETA to solve the EEG inverse problem.Recalling (2. 19) and (2.27), we define the local error as where and We note that L(r k ), L S (r k ) ∈ R 3×3 .
Lemma 3.4.For the difference of the reciprocals of ∥L(r k )∥ F and ∥L S (r k )∥ F the following bound holds where ξ = min 1⩽k⩽n σ min (L(r k )).
The proof is in the appendix.
Proposition 3.5.Let Q wMNE be the neuroelectric current estimate obtained with the randomized wMNE.It holds ) The proof is in the appendix.
Proof.When k ∈ I c we have The first term in the r.h.s. is related to the wMNE reconstruction error while the second term is related to the error produced by the sampling procedure applied to the lead field matrix.From (2.7) it follows ∥q k wMNE ∥ − ∥q k wMNE ∥ ⩽ ∥q k wMNE − qk wMNE ∥.Thus, we can use the results in theorem 3.3 to bound the first term.
As for the second term, since ∥q k wMNE ∥ ⩽ ∥ Q wMNE ∥, combining lemma 3.4 and proposition 3.5 we get When k / ∈ I c , the bound easily follows recalling that ∥q k wMNE ∥ ⩽ ∥Q wMNE ∥, and combining lemma 3.4 and proposition 3.2.

Results
To show the effectiveness of randomized wMNE and randomized sLORETA methods we conducted some numerical experiments on synthetic EEG data.The latter have been generated by solving the forward EEG problem in OpenMEEG software [32] using a source space obtained by an average subject [33] containing 25 000 points and an EEG helmet with 65 sensors at an average distance of 10 cm from the origin of the head coordinate system.One single dipole was randomly selected from the source space and activated by a unitary current dipole moment with fixed orientation.The measurement vector b was then obtained by applying (2.8) and by adding Gaussian noise with SNR = 10 dB and SNR = 20 dB.A coarser source space containing 15 000 points was used to solve the inverse problem, that is the lead field matrix is M ∈ R m×n , with n = 15000 and m = 65.The inverse problem was solved by wMNE and sLORETA methods implemented in Brainstorm software [34].To improve the accuracy, the inverse problem was solved N = 10 times, each time using a different sampling matrix.In our experiments, no significant improvement has been observed by increasing N.
In section 4.1, we present a comprehensive set of numerical tests in the case of fixed orientation.In section 4.2, we present some numerical tests for the free orientation case.

Fixed orientation
Since the error estimates E wMNE and E sL depend on the condition number κ(G wMNE ) and on the quantities ξ and µ (see theorems 3.3 and 3.6), first of all we computed their values for different values of the regularization parameter α.We note that in Brainstorm the regularization parameter is included in the matrix R; the default value is α = 3.The computed values of κ(G wMNE ) and ξ are listed in table 1.As expected, the condition number κ(G wMNE ) approaches 1 while the regularization parameter increases.Conversely, ξ decreases.Thus, the error estimates E wMNE and E sL decrease for increasing α.Table 2 shows how µ (see proposition 3.1) changes with respect to the depth weighting parameter ρ.We note that the default in Brainstorm is ρ between 0 and 1.As expected, the values follow an exponential growth.Then, we studied how the reconstruction errors E wMNE and E sL , as defined in (3.6) and (3.10), depend on the number of sampled columns c.The numerical experiments have been carried out for different values of c, ranging from a minimum of c = 500 columns to a maximum of c = 15000 columns.In this test we used the default values for α and ρ.The averaged absolute errors ∥E wMNE ∥ and ∥E sL ∥ are shown in figure 1 while the relative error ∥EwMNE∥ ∥QwMNE∥ is shown in figure 2. To reduce the computational cost of the dipole sampling method, we sampled the columns with replacement.For comparison, we evaluated the errors also in case of sampling without replacement.
Finally, since our aim is to show that the random sampling method is efficient in localizing neural sources, we evaluated the distance localization error (DLE), i.e. the distance between the real source and the center of mass of the reconstructed current distribution.The DLE is a common metric used in the MEG/EEG community to test the accuracy of the reconstructed neural current [35].For the random dipole sampling method, we compute the center of mass of the source localizations obtained in 10 different runs of the method.The number of  sampled columns c ranges from 100 to 2000 and the results are averaged over 1000 simulations.Figure 3 shows the boxplots for superficial, middle and deep sources when SNR = 10 dB and SNR = 20 dB while the average DLE and the standard deviation are listed in table 3.For comparison, the results obtained using the full source space are also shown.
To analyze the effect of the source space size, we performed numerical tests for n = 2000, 4000, 8000.Figures 4 and 5 show the boxplots of DLE across 100 simulations for superficial, middle and deep sources when SNR = 10 dB and SNR = 20 dB.The average DLE and the standard deviation are listed in tables 4-6.
In figure 6 the median of the condition number κ(G wMNE ) across 100 simulations for different values of the number of sampled columns c is plotted.First column refers to superficial sources (at least 50 mm from the origin of the head coordinate system), second column to middle sources (between 30 and 50 mm from the origin of the head coordinate system), third column refers to deep sources (at less then 30 mm from the origin of the head coordinate system).
Table 3.The average DLE and the standard deviation σ (mm) for n = 15 000 and different source depth (LOC = S (superficial sources at least 50 mm from the origin of the head coordinate system), M (middle sources between 30 and 50 mm from the origin of the head coordinate system), D (deep sources at less then 30 mm from the origin of the head coordinate system)).First column refers to superficial sources (at least 50 mm from the origin of the head coordinate system), second column refers to middle sources (between 30 and 50 mm from the origin of the head coordinate system), third column refers to deep sources (at less then 30 mm from the origin of the head coordinate system).

Free orientation
In case of free orientation we define the local reconstruction error for wMNE as the estimate in theorem 3.3 holds also in this case.
For sLORETA the local reconstruction error is defined as where z k and zk are given in (2.13) and (2.23), respectively.It can be shown that The first term depends on the wMNE local reconstruction error while the second term is related to the error of the covariance matrix.Since ∥q k wMNE ∥ ⩽ ∥Q wMNE ∥ and ∥q k wMNE ∥ ⩽ ∥ QwMNE ∥, estimates in proposition 3.5 holds.Moreover, Table 4.The average DLE and the standard deviation σ (mm) for n = 2000 and different source depth (LOC = S (superficial sources at least 50 mm from the origin of the head coordinate system), M (middle sources between 30 and 50 mm from the origin of the head coordinate system), D (deep sources at less then 30 mm from the origin of the head coordinate system)).The average DLE and the standard deviation σ (mm) for n = 4000 and different source depth (LOC = S (superficial sources at least 50 mm from the origin of the head coordinate system), M (middle sources between 30 and 50 mm from the origin of the head coordinate system), D (deep sources at less then 30 mm from the origin of the head coordinate system)).

SNR LOC inversion method
where ξ = min l(σ min (Σ k QwMNE ), σ min ( Σk QwMNE ) .Thus, the sLORETA local error behaves as the wMNE local error.Table 6.The average DLE and the standard deviation σ (mm) for n = 8000 and different source depth (LOC = S (superficial sources at least 50 mm from the origin of the head coordinate system), M (middle sources between 30 and 50 mm from the origin of the head coordinate system), D (deep sources at less then 30 mm from the origin of the head coordinate system)).To show the effectiveness of the random dipole sampling method even in the free orientation case, we performed numerical tests with n = 15 000 and SNR = 20 dB.The boxplots of the DLE in figure 7 show that the DLE is comparable to that obtained in the fixed orientation case.

Related methods
Monte Carlo methods have a long history dating back to the middle of the twentieth century [36].Since then, these methods have been widely used in scientific computing, for instance for integrating differential equations or probability distributions.In numerical linear algebra randomness is commonly used to initialize iterative methods [28].The first papers in which randomized methods were used to construct low-rank approximations of matrices are from the late 1990s [37,38].Nowadays, these methods are successfully used in several fields, such as computer science, numerical linear algebra, regression problems.For a comprehensive review on randomized methods we refer the reader to [15,16].
There are many ways randomness comes into play in the solution of the MEG/EEG inverse problem.In the Bayesian framework Markov Chain Monte Carlo methods are usually used to sample posterior densities (see [3,39] and references therein).In parametric methods the number of sources is modeled as a random variable while their locations and orientations are modeled as random vectors (see [31] and references therein).In the randomized multiresolution scanning method [23,24] randomness is used to generate multiple randomized decompositions of the source space at different resolution levels.
The random dipole sampling method described in section 2.3 uses a different strategy.It samples a small set of dipoles in the source space using the uniform probability distribution.Then, it solves an inverse problem where the unknowns are the dipole moments on the drawn (fixed) locations.The rationale behind the method is that we can assume the neuroelectric current J(r) to be spatially sparse so that it can be represented using just few sources [40,41], i.e.
The uniform distribution is a good choice when no prior information on the spatial distribution of the neuroelectric current is available.On the other hand, using different probability distributions would increase the computational cost.Now, sampling c dipoles means to sample 3c columns in the lead field matrix since we select only the c sub-matrices M(r k ) corresponding to the dipoles with indices k ∈ I c .Thus, the random dipole sampling can be seen as a randomized column sampling procedure and it reduces to the column sampling introduced in [25,42] in case of fixed orientation.For a given matrix A ∈ R m×n the column sampling is optimal when the columns are sampled using the probability distribution where A (k) denotes the kth column of A [15].This is the same as dividing each column of the matrix A by p k and then sampling the columns using the uniform distribution.We note that usually in wMNE the columns of the lead field matrix are normalized by a depth weighting that is related to the norm of the columns.This is evident when writing the wMNE solution (2.12) in the equivalent form [2] Q wMNE = RM T MRM T + I −1 b, and justifies the use of the uniform distribution in the random dipole sampling method described in section 2.3.

Our findings
As shown in figures 1 and 2, the error with replacement and without replacement are comparable for lower values of c and the difference only gets relevant if c > 5000.In fact, in case of sampling without replacement the error goes to zero as the number of sampled columns approaches the maximum number of columns n, while in case of sampling with replacement the error decreases to a minimum located around c = 6000 for wMNE and c = 9000 for sLORETA.For higher values of c it grows slightly.This does not pose a strict restriction on the use of the random dipole sampling method, especially in real-time applications, as this method is aimed to significantly reduce the dimensionality of the EEG inverse problem, thus the goal is to use it for c ⩽ 2000.The behavior of the error (see figures 1 and 2) is in agreement with the estimates given in section 3.In fact, when the number of sampled columns c is much smaller than n, due to factor max( n c ν − 1, 1) the error decreases rapidly as c increases reaching a minimum approximately for n c ν = 2.Then, the error increases slowly, roughly as ν.Instead, in case of sampling without replacement the error decreases monotonically as c increases and is equal to zero for c = n.
The numerical tests show that the choice 500 ⩽ c ⩽ 2000 is a good compromise between accuracy and computational cost.This is more evident looking at the results shown in figures 3-5 and tables 3-6.When n = 15000 the DLE obtained by the random dipole sampling method with c = 500, 1000, 2000 is comparable to the error obtained using the full source space and is quite often even lower (see figure 3 and table 3).In general, sLORETA performs better than wMNE, regardless of source depth, in agreement with known results in the literature [43].This is true for both SNR = 10 dB and SNR = 20 dB, causing the higher noise a slightly larger error, but still comparable with that of the full source space case.We note that even in the case of a smaller source space, the DLE values for the two noise levels are comparable.
When the size of the source space is reduced (see figures 4 and 5 and tables 4-6 where n = 2000, 4000, 8000) the DLE for sLORETA is lower when the number of sampled columns is c = 2000, i.e. c is in the order of n.The DLE for wMNE is lower both when n = 2000, for all values of c and even in case of full source space, and when c = 100, for all values of n.This means that wMNE is able to localize sources well not only when the source space has a low dimension, but also when few sources of a large source space are sampled by the random dipole sampling method.It is known that a low-dimensional source space can improve the localization of sparse neural sources [44] but this comes at the expense of low resolution.The random dipole sampling method has the same effect, with the advantage that the resolution can be easily improved by averaging the results obtained using different sampled source spaces.Interestingly, the reconstruction error is lower especially for deep and middle sources, as shown in figure 8, where the norm of the local reconstruction error, averaged for different source depths, is displayed for both wMNE and sLORETA.This suggests that the sparsity induced by the random dipole sampling method is more effective for deep sources, according to the results in the literature (see [24,44]).
Finally, we note that the behavior of the localization error can be related to the conditioning of the matrix G wMNE .In fact, the values of κ(G wMNE ) displayed in figure 6 show that the condition number increases when c increases, being smaller when n = 2000 or when c ⩽ 2000.

Conclusion
In this paper, we analyzed in detail how to apply the random dipole sampling method for solving the inverse EEG problem.To the best of our knowledge, this is the first time that the random dipole sampling method has been used along with wMNE and sLORETA, two inversion methods well-known in the literature.The numerical tests have shown that the localization of neural sources produced by randomized wMNE and randomized sLORETA is comparable to that obtained when using the full source space.This suggests that the randomized version of wMNE and sLORETA can be successfully used in real-time applications, possibly with wearable systems, for which it is important to reduce the computational cost of the inversion step.
To this end, it is important to analyze the performance of the random dipole sampling method in solving the inverse EEG problem in case of real data.In addition, the reconstruction error estimates we provided are deterministic.It would be interesting to obtain probabilistic error estimates.These topics will be the subject of future work.
The last inequality is a consequence of the fact that M T S b ∈ R c while M T b R n .The inequality for ∥Q wMNE ∥ can be proved in a similar way.

. 6 ) 1 .
Remark In case of uniform sampling without replacement we have ∥S T ∥ = 1 and ∥I − n c I S ∥ ⩽ ∥I − n c I S ∥ F = 2(n − c) so that the estimate in (3.3) becomes

Figure 1 .
Figure 1.(a) The absolute error ∥E wMNE ∥ versus the number of sampled columns c, in the case of sampling with replacement and sampling without replacement.(b) The absolute error ∥E sL ∥ versus the number of sampled columns c, in the case of sampling with replacement and sampling without replacement.

Figure 2 .
Figure 2. The relative error ∥EwMNE∥ ∥QwMNE∥ versus the number of sampled columns c, in the case of sampling with replacement and sampling without replacement.

Figure 3 .
Figure 3. Boxplots of DLE across 1000 simulations for n = 15 000 when SNR = 10 dB (first row) and SNR = 20 dB (second row).First column refers to superficial sources (at least 50 mm from the origin of the head coordinate system), second column to middle sources (between 30 and 50 mm from the origin of the head coordinate system), third column refers to deep sources (at less then 30 mm from the origin of the head coordinate system).

Figure 4 .
Figure 4. Boxplots of DLE across 100 simulations when SNR = 10 dB for n = 2000 (first row), n = 4000 (second row), n = 8000 (third row).First column refers to superficial sources (at least 50 mm from the origin of the head coordinate system), second column refers to middle sources (between 30 and 50 mm from the origin of the head coordinate system), third column refers to deep sources (at less then 30 mm from the origin of the head coordinate system).

Figure 5 .
Figure 5. Boxplots of DLE across 100 simulations when SNR = 20 dB for n = 2000 (first row), n = 4000 (second row), n = 8000 (third row).First column refers to superficial sources (at least 50 mm from the origin of the head coordinate system), second column refers to middle sources (between 30 and 50 mm from the origin of the head coordinate system), third column refers to deep sources (at less then 30 mm from the origin the head coordinate system).

Figure 6 .
Figure 6.The median of the condition number κ(G wMNE ) across 100 simulations versus the number of sampled columns c for different values of n.

Figure 7 .
Figure 7. Boxplots of DLE across 100 simulations with n = 15 000 and SNR = 20 dB in the free orientation case for superficial (a), middle (b), and deep (c) sources.

Figure 8 .
Figure 8.The averaged absolute error ∥E wMNE ∥ for wMNE (a) and sLORETA (b) versus the number of sampled columns c, in the case of sampling with replacement for superficial sources, middle sources and deep sources.

Table 1 .
Condition number κ(G wMNE ) and the quantity ξ (see lemma 3.4) for different values of the regularization parameter α.

Table 2 .
The quantity µ (see proposition 3.1) for different values of the depth weighting parameter ρ.