A multigrid solver in EEG source analysis applying
the finite difference method

B. Vanrumste¹, M. Mohr² and P. J. Bones^1
1Department of Electrical and Computer Engineering, University of Canterbury,
Private Bag 4800, Christchurch, New Zealand.
²System Simulation Group of the Computer Science Department, Fredrich-Alexander-University,
Cauerstraße  6, D-91058 Erlangen, Germany.

Abstract: The Algebraic Multigrid method for a discretisation based on cell-centered finite difference is applied to solve the forward problem in EEG source analysis.  We found that, depending on the computer platform, the method performs 1.8-5.6 times faster than the successive over relaxation (SOR) method often employed in this context.

INTRODUCTION

It is becoming standard nowadays to base EEG source analysis on the individual geometry of a patient's head, as can be obtained from magnetic resonance imaging (MRI). The MRI is segmented into different compartments. Each compartment is allocated an electrical conductivity depending on the tissue it represents. In these realistically shaped volume conductor models the potentials at the computational points are calculated for a given electrical source. Numerical methods, such as the Boundary Element Method (BEM), the Finite Element Method (FEM) and the Finite Difference Method (FDM), are applied for this purpose. These methods transfer the differential equation of Poisson into a linear equation at the computational points. Hence we obtain a linear system Ax=b which needs to be solved, where x represents the potentials at the computational points, and b corresponds to the electrically active sources.

We will further focus on the FDM which applies a volume based discretisation technique. FDM gives us the ability to distinguish a large number of conducting compartments. It also allows an easy introduction of anisotropic conductivity. Due to the volume based discretisation in FDM a large number of computational points (> 500 000) is involved.  Since the system matrix A is sparse, iterative solvers are considered to solve the linear systems. These solvers compute x only for the given right hand side. Subsequently, for another right hand side, the iterative solver has to be re-applied. However they only require the storage of non-zero matrix entries, and, in general, the amount of work per iteration is dominated by the cost of one or two matrix vector multiplications.

Multigrid methods in general are known to be very efficient iterative solvers for the linear system associated with the equation of Poisson. The basic idea is to notionally split the error of an approximate solution into smooth and rough components, representing low- and high-frequency Fourier modes. The rough components are reduced in magnitude on the original (fine) grid by applying a limited number of steps of some iteration scheme, like e.g. Gauss-Seidel or SOR. In the remaining error the smooth components are dominant. For these components a correction is computed by transferring the problem to a coarser grid with a larger mesh size. The correction is then interpolated back to the fine grid and added to the approximation. This scheme is repeated until convergence. The problem on the coarser grid is solved recursively by application of the same principle.

The difficulty in multigrid is finding the proper inter-grid transfer functions. In the case of complex geometries and/or jumping coeficients, this is not trivial. Therefore the idea of algebraic multigrid methods (AMG), as illustrated  in [1] ,  is again attracting increased attention. Here a grid hierarchy and inter-grid transfer functions are derived automatically from the problem matrix A.

We compare on different platforms the run time of the AMG method with that of the SOR method. 

METHODS

The model of the human head is constructed from a MRI. In a segmentation process each voxel is assigned to one of the four compartment types scalp, skull, brain, and air. Constant isotropic conductivities are assumed in each of them. The conductivity ratios applied are  σ_scalp : σ_skull : σ_brain = 16 : 1 : 16, and σ_air = 0. Figure 1 presents a coronal slice of the volume conductor model. An overview of the properties of the data set and of the resulting linear system is given in table 1.

 TABLE 1. Data set details.

The amount of computational work involved with EEG source analysis can be drastically reduced with the help of the reciprocity theorem of Helmholtz [2] . It allows a forward evaluation by interpolation from potential distributions computed in a preparatory step. This step consists of the solution of N-1 equations of Poisson corresponding with the N electrodes used. Each solution is generated by solving Poisson’s equation with one electrode chosen as unit sink and the other one as unit source. Figure 1 presents the resulting potential fields for a current sink at C3 and a current source at T4. We use 27 electrodes hence 26 Poisson’s equations need to be solved.

Figure 1. A coronal slice of the head model with the reciprocal current sources and the equipotential lines.

Poisson's equation with pure Neumann boundary conditions possesses a unique solution only up to an additive constant. The FDM leads to a singular matrix A. One can either apply the iterative methods directly to the singular system, in which case they will converge to one of the infinitely many solutions, or introduce an additional constraint that fixes one element by introducing a reference potential, which yields a regular problem and then apply the methods. We found that applying the singular matrix yields equivalent or smaller run times compared to applying the regular matrix, see [3] .

Both tested algorithms depend on parameters. In case of the SOR method we have to specify the over-relaxation parameter ω.  The optimal value is found to be 1.935. In the AMG algorithm there is quite a number of parameters and algorithmic variants, that have to be fixed, before the method can be applied. These include for example: cycle type, number of smoothing steps, and so on. We have used a V-cycle with one pre- and one post-smoothing step. The smoother was a hybrid Gauss-Seidel / Jacobi method. Construction of the grid hierarchy was performed with Ruge-Stüben coarsening, see e.g. [1] . Here the threshold parameter α, which allows for the distinction of weak and strong inter-node connections, has been chosen as α = 0.05.

The property that is of primary interest to the user is run time. Besides the numerical characteristics of an algorithm, also its suitability for modern computer architectures influences the run time Therefore the run times of both  algorithms have been tested on three different architectures, a 700MHz AMD Athlon, a 500MHz Alpha A21264 and a 1500MHz Pentium IV. We have implemented the SOR method ourselves and used the BoomerAMG implementation from [4]. The measured run times consist of the time spent in the setup phase and for the solution of 26 forward problems. A problem was considered solved, as soon as after i iterations the residual (rⁱ=b-Axⁱ) of the approximate solution measured in the Euclidean norm became smaller than 10^-8.  The setup times are negligible for the SOR, but in the AMG case, the grid hierarchy adds considerably to the total cost. This varies between 5 and 11 % of the total time depending on architecture and problem size.

RESULTS

Table 2 presents the run times in seconds. We see, that the AMG solver performs the task fasted for the three platforms.  The SOR performs the task 1.8 – 5.6 times slower than the AMG method, depending on the platform used.  We can also see that the run times of a specific solver are platform dependent.

TABLE 2. Run times in seconds for the complete problem.
This includes the setup phase and the solution of 26 singular forward problems.

DISCUSSION

A important group of iterative solvers are Conjugate Gradient (CG) type methods. As was shown in [3] AMG outperforms the CG method, and CG method preconditioned with symmetric SOR. It remains to be seen however, if this superiority of AMG also holds for more sophisticated preconditioners.

A multigrid application for electrical nerve stimulation is  presented by Hoekema et al. [5] . It was shown that with this technique the computing time was reduced significantly compared to the conventional iterative methods. These results are in concordance with our findings.

REFERENCES

[1] J. W. Ruge and K. Stüben, "Algebraic multigrid (AMG)," in Multigrid Methods, vol. 3, Frontiers in Applied Mathematics, S. F. McCormick, Ed.: SIAM, 1987, pp. 73-130.

[2] B. Vanrumste, G.Van Hoey, R. V. d. Walle, M. D'Havé, I. Lemahieu, and P. Boon, "The validation of the finite difference method and reciprocity for solving   the inverse problem in EEG dipole source analysis," Brain Topography, vol. 14, pp. 83-92, 2001.

[3] M. Mohr, Comparison of solvers for a bioelectric field problem,"   Lehrstuhl für Informatik 10 (Systemsimulation),   Friedrich-Alexander-Universitat Erlangen-Nurnberg,  Technical Report 01-2 2001 2001.

[4] V. E. Henson, U. Meier Yang, "BoomerAMG: A parallel algebraic multigrid solver and preconditioner", Applied Numerical Mathematics, vol. 41, pp.155-177, 2002.

[5] R. Hoekema, K. Venner, J. J. Struijk, and J. Holsheimer, "Multigrid solution of the potential field in modeling electrical nerve   stimulation," Computers and Biomedical Research, vol. 31, pp. 348-362, 1998.