The EEG simulation of a multivariate case can be described at a time instant as (1):

where y is an 𝑚 𝖷 1 vector describing the measurements of 𝑚 EEG electrodes at a time instant; 𝑥, an 𝑛 𝖷 1 vector describing the brain activity at 𝑛 possible locations; and 𝑀, an 𝑚 𝖷 𝑛 matrix called the lead-field matrix, which relates the neural activity inside the brain with the EEG signals. (1) can be rewritten as a sub-set of block matrices (2):

Where 𝑀^𝑗 is the 𝑗 -𝑡ℎ matrix block of dimension 𝑚 𝖷 𝑛₁, being if 𝑛 can be written as an integer number 𝑁 of 𝑛₁

A parallel version of FISTA optimization, which can be used since the structure of the estimated brain activity 𝑥_𝑘 is sparse, is presented in ^[12]. The parallel FISTA algorithm is described by the following set of equations (3), (4), (5), (6):

Where 𝑘 = 1, … ,𝐾 is the iteration index; , the 𝑗 -𝑡ℎ block activity computed in parallel processes, and according to ^[11]; and the proximal operator in terms of the norm and the regularization parameter 𝜆 ^[13] defined as (7):

The resulting estimated activity is computed at the last iteration as (8):

This work proposes a parallel method for a Tikhonov-like solution as an extension of the parallel FISTA algorithm with a proximal operator based on the norm ^[11], as follows (8), (9), (10), (11), (12):

Where 𝑘 = 1, … ,𝐾 is the iteration index; , the 𝑗 -𝑡ℎ block activity computed in parallel processes; and , the proximal operator in terms of the norm.

Two high-resolution head models were studied in this work for evaluating the performance of the algorithms in regular and parallel implementations. The first model, named the New York head model, was used with 𝑛 = 2004 (2 k), 𝑛 = 10016 (10 k), and 𝑛 = 74382 (75 k) sources and 𝑚 = 230 EEG channels, as described in^[14]. The second model, named the Belgian head model, was used with 𝑛 = 8000 (8 k), 𝑛 = 20000 (20 k), and 𝑛 = 70000 (70 k) sources and 𝑚 = 69 EEG channels. A detailed description of the second model can be found in ^{[14], [15]}. A visual representation of both head models with their different numbers of sources is shown in Fig. 1.

Fig. 1 Source: Authors.

Noticeably, as the neural activity is estimated at each vertex of the brain model, a higher spatial resolution (high 𝑛 values) can reduce the spatial error that occurs when a low-resolution brain model is considered as a reference for an inverse problem solution.

The computer processor used for the experiments is a 12-core Intel Xeon Silver 4116 with 2.1 GHz of speed and 64 GB of RAM. The algorithms were implemented employing the Open Message Passing Interface (Open MPI). To this end, Equations (4) and (8) were solved by means of the MPI-ALLreduce function with the sum operator. It is worth mentioning that the GNU Scientific Library (GSL) was used for the calculations.

In order to evaluate the performance of the algorithms, we propose the following simulation benchmark: A static active source placed in a random brain cortex location is simulated in both high-resolution models with different numbers of sources, as shown in the first column of Fig. 2. For the sake of simplicity in the parallel implementations, here we only simulated a single time instant. Then, we reconstructed the brain activity using the FISTA and Tikhonov-like methods described in the previous sections with one, two, four, six, eight, ten, and twelve simultaneous processes. The number of parallel processes N was selected from one to twelve, which is the maximum number of cores of the computer processor used for the calculations.

Fig. 2 Source: Authors.

We used two performance measurements:

i) Computational time (t): It helps us to quantify the computational improvement achieved by the parallel solutions. Here, we measured the estimation time in seconds for the 500 iterations cycle used in Tikhonov, including the communication time among processes. In addition, a comparison of the algorithms, using the speedup measurement, is also included.

ii) Relative error (𝑒): It evaluates if the parallel solutions worsen the brain activity estimation. This measurement is computed as (13):

Where 𝑥 and are the 𝑛 𝖷 1 vector comprising the simulated and estimated brain activity.

The performance of the algorithms, in terms of computational time, is shown in Fig. 3. In the latter, we notice a monotonic computational time reduction if the number of parallel processes is increased, regardless of the type of algorithm used to solve the inverse problem (namely, FISTA or Tikhonov-like). Moreover, we observe that the difference in the computational times with the three n values (considered when using a single process) is significantly reduced as the number of processes is increased.

Fig. 3 Source: Authors.

Specifically, with a single process, using n ≈ 70 k requires about four times the computational time needed for n ≈ 2 k, and about three times that for n ≈ 10 k. However, with 12 parallel

processes, using n ≈ 70 k requires two times the computational time used for n ≈ 2 k or n ≈ 10 k.

Furthermore, noticeably, the computational time required to solve the inverse problem using parallel processes for high-resolution models (namely, n ≈ 70 k) is equal to or even lower than that needed to solve the optimization with a single process for a reduced head model (namely, n ≈ 10 k). For example, the same reduced computational time is required to solve the inverse problem of a model with 20 k sources implementing a single process or that of a model with 70 k sources implementing 4 or more parallel processes.

This is a remarkable result since a high-resolution model of 70 k sources can be used instead of a model of 20 k sources in the same or even less computational time. In order to validate this result, an additional calculation, using the speedup measurement, is presented in Fig. 4.

Fig. 4 Source: Authors.

Moreover, to validate whether the parallel solutions deteriorate the source reconstruction quality, Fig. 5 shows the relative error for all the models considered in this work. From such Figure, we can note that the relative error does not depend on the number of parallel processes but on the number of sources of the head model.

Fig. 5 Source: Authors.

An additional example of the localization error reduction that can be obtained using high-resolution models can be seen in Fig. 6. The latter compares parallel FISTA and Tikhonov solutions for the Belgian head model with 8 k, 20 k, and 70 k sources. By increasing the number of sources of the head model, the localization error is reduced since the distance between sources is lower in high-resolution models.

Fig. 6 Source: Authors.

This property is specially observed in the results of the parallel FISTA method in Fig. 3 and Fig. 5, where it is noticeable that, by using parallel processes for the solution of optimization problems, a reduction in computational load can be achieved.

Therefore, the same computational time is required to solve a low-resolution model implementing one process or a high-resolution model implementing several parallel processes; and the brain activity estimation error is not increased.

The speedup measurement was used in order to validate the computational time results, especially for the models with 70 k and 75 k sources. Fig. 4 presents the speedup measurement of the FISTA algorithm. It can be seen that the speedup increases along with number of parallel processes. It is worth mentioning that, in order to compute the speedup measurement, 45 % of the code is parallelized.