A Calderon Regularized Symmetric Formulation for the Electroencephalography Forward Problem

TL;DR

This paper introduces a Calderon regularized symmetric formulation for the EEG forward problem, improving numerical stability and conditioning, especially with high discretization density and conductivity contrast.

Contribution

It proposes a novel regularized symmetric formulation using Calderon identities, resulting in a stable, well-conditioned second kind operator for EEG forward modeling.

Findings

The new formulation is stable under high discretization density.

It maintains compatibility with existing EEG software.

Numerical results validate improved stability and conditioning.

Abstract

The symmetric formulation of the electroencephalography (EEG) forward problem is a well-known and widespread equation thanks to the high level of accuracy that it delivers. However, this equation is first kind in nature and gives rise to ill-conditioned problems when the discretization density or the brain conductivity contrast increases, resulting in numerical instabilities and increasingly slow solutions. This work addresses and solves this problem by proposing a new regularized symmetric formulation. The new scheme is obtained by leveraging on Calderon identities which allow to introduce a dual symmetric equation that, combined with the standard one, results in a second kind operator which is both stable and well-conditioned under all the above mentioned conditions. The new formulation presented here can be easily integrated into existing EEG imaging packages since it can be obtained with the same computational technology required by the standard symmetric formulation. The performance of the new scheme is substantiated by both theoretical developments and numerical results which corroborate the theory and show the practical impact of the new technique.

Tables

Table 1. Table 1: Asymptotic behaviour of the different coefficients of Z q subscript 𝑍 𝑞 Z_{q} .

Coefficient	${\sigma }_i\to \infty$	${\sigma }_i\to 0$
$R_{ij}$	$0$	$0$
$P_i$	$1$	$0$
$P_i^{\ast }$	$0$	$1$

Table 2. Table 2: Memory and computation time information for computing a lead field matrix using the reciprocity method

Method	Memory (GB)	Operator Time (s)	Solution Time (s)	Lead Field Time (h)
DI Dense Symmetric	16.234	10845.67	54609.99	18.19
CGS Dense Symmetric	16.234	10845.67	7294.61	45.56
CGS Compressed Symmetric	1.254	1436.63	2322.89	13.95
CGS Compressed Calderon-Symmetric	2.542	7888.86	62.40	2.56

Equations

$$\sigma \Delta V=\nabla \cdot \mathbf{j}\text{in each Ωi}$$

$${\left[V\right]}_i=0\forall i\le N$$

$${\left[\sigma \mathbf{n}\cdot \nabla V\right]}_i=0\forall i\le N$$

$$G(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{1}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime }|}$$

$$S\Psi (r)$$

$$D\Phi (r)$$

$$D^{\ast }\Psi (r)$$

$$N\Phi (r)$$

$$\begin{gathered} {\sigma }_{i+1}^{-1}{\left(v_{i+1}\right)}_{{\Gamma }_i}-{\sigma }_i^{-1}{\left(v_i\right)}_{{\Gamma }_i}= \\ {D}_{i,i-1}{V}_{i-1}-2D_{ii}{V}_i+D_{i,i+1}{V}_{i+1}-{\sigma }_i^{-1}{S}_{i,i-1}{p}_{i-1}+\left({\sigma }_i^{-1}+{\sigma }_{i+1}^{-1}\right)S_{ii}{p}_i-{\sigma }_{i+1}^{-1}{S}_{i,i+1}{p}_{i+1} \end{gathered}$$

$$\begin{gathered} {\left({\partial }_n{v}_{i+1}\right)}_{{\Gamma }_i}-{\left({\partial }_n{v}_i\right)}_{{\Gamma }_i}= \\ {\sigma }_i{N}_{i,i-1}{V}_{i-1}-\left({\sigma }_i+{\sigma }_{i+1}\right)N_{ii}{V}_i+{\sigma }_{i+1}{N}_{i,i+1}{V}_{i+1}-D_{i,i-1}^{\ast }{p}_{i-1}+2D_{ii}^{\ast }{p}_i-D_{i,i+1}^{\ast }{p}_{i+1} \end{gathered}$$

$$\begin{array}{rl}& \underset{Z}{\underbrace{\left[\begin{array}{cccccccc}\left({\sigma }_1+{\sigma }_2\right)N_{11}& -2D_{11}^{\ast }& -{\sigma }_2{N}_{12}& D_{12}^{\ast }& & & & \\ -2D_{11}& \left({\sigma }_1^{-1}+{\sigma }_2^{-1}\right)S_{11}& D_{12}& -{\sigma }_2^{-1}{S}_{12}& & & & \\ -{\sigma }_2{N}_{21}& D_{21}^{\ast }& \left({\sigma }_2+{\sigma }_3\right)N_{22}& -2D_{22}^{\ast }& \cdots & & & \\ D_{21}& -{\sigma }_2^{-1}{S}_{21}& -2D_{22}& \left({\sigma }_2^{-1}+{\sigma }_3^{-1}\right)S_{22}& \cdots & & & \\ & & ⋮& ⋮& \ddots & & & \\ & & & & & \left({\sigma }_{N-1}+{\sigma }_N\right){\mathbf{N}}_{N-1,N-1}& -2{{\mathbf{D}}^{\ast }}_{N-1,N-1}& -{\sigma }_N{\mathbf{N}}_{N-1,N}\\ & & & & & -2D_{N-1,N-1}& \left({\sigma }_{N-1}^{-1}+{\sigma }_N^{-1}\right)S_{N-1,N-1}& D_{N-1,N}\\ & & & & & -{\sigma }_N{N}_{N,N-1}& D_{N,N-1}& {\sigma }_N{N}_{N,N}\end{array}\right]}}\\ & \underset{x}{\underbrace{\left[\begin{array}{c}{V}_1\\ p_1\\ V_2\\ p_2\\ ⋮\\ p_{N-1}\\ V_N\end{array}\right]}}=\underset{b}{\underbrace{\left[\begin{array}{c}{\left({\partial }_n{v}_1\right)}_{{\Gamma }_1}-{\left({\partial }_n{v}_2\right)}_{{\Gamma }_1}\\ {\sigma }_2^{-1}{\left(v_2\right)}_{{\Gamma }_1}-{\sigma }_1^{-1}{\left(v_1\right)}_{{\Gamma }_1}\\ {\left({\partial }_n{v}_2\right)}_{{\Gamma }_2}-{\left({\partial }_n{v}_3\right)}_{{\Gamma }_2}\\ {\sigma }_3^{-1}{\left(v_3\right)}_{{\Gamma }_2}-{\sigma }_2^{-1}{\left(v_2\right)}_{{\Gamma }_2}\\ ⋮\\ {\sigma }_N^{-1}{\left(v_N\right)}_{{\Gamma }_N-1}-{\sigma }_{N-1}^{-1}{\left(v_{N-1}\right)}_{{\Gamma }_N-1}\\ {\left({\partial }_n{v}_N\right)}_{{\Gamma }_N}\end{array}\right]}}\end{array}$$

$$[{\mathbf{D}}_{ij}{]}_{kl}={\int }_{t_k}{D}_{ij}(P_{1l})P_{0k}(r)dr$$

$$[{\mathbf{S}}_{ij}{]}_{kl}={\int }_{t_k}{S}_{ij}(P_{0l})P_{0k}(r)dr$$

$$[{\mathbf{N}}_{ij}{]}_{kl}={\int }_{{\mu }_{P_{1i}}}{N}_{ij}(P_{1l})P_{1k}(r)dr$$

$$[{{\mathbf{D}}^{\ast }}_{ij}{]}_{kl}={\int }_{{\mu }_{P_{1i}}}{D}_{ij}^{\ast }(P_{0l})P_{1k}(r)dr.$$

$$\begin{array}{rl}& \mathbf{Z}=\\ & \left[\begin{array}{cccccccc}\left({\sigma }_1+{\sigma }_2\right){\mathbf{N}}_{11}& -2{{\mathbf{D}}^{\ast }}_{11}& -{\sigma }_2{\mathbf{N}}_{12}& {{\mathbf{D}}^{\ast }}_{12}& & & & \\ -2{\mathbf{D}}_{11}& \left({\sigma }_1^{-1}+{\sigma }_2^{-1}\right){\mathbf{S}}_{11}& {\mathbf{D}}_{12}& -{\sigma }_2^{-1}{\mathbf{S}}_{12}& & & & \\ -{\sigma }_2{\mathbf{N}}_{21}& {{\mathbf{D}}^{\ast }}_{21}& \left({\sigma }_2+{\sigma }_3\right){\mathbf{N}}_{22}& -2{{\mathbf{D}}^{\ast }}_{22}& \cdots & & & \\ {\mathbf{D}}_{21}& -{\sigma }_2^{-1}{\mathbf{S}}_{21}& -2{\mathbf{D}}_{22}& \left({\sigma }_2^{-1}+{\sigma }_3^{-1}\right){\mathbf{S}}_{22}& \cdots & & & \\ & & ⋮& ⋮& \ddots & & & \\ & & & & & \left({\sigma }_{N-1}+{\sigma }_N\right){\mathbf{N}}_{N-1,N-1}& -2{{\mathbf{D}}^{\ast }}_{N-1,N-1}& -{\sigma }_N{\mathbf{N}}_{N-1,N}\\ & & & & & -2{\mathbf{D}}_{N-1,N-1}& \left({\sigma }_{N-1}^{-1}+{\sigma }_N^{-1}\right){\mathbf{S}}_{N-1,N-1}& {\mathbf{D}}_{N-1,N}\\ & & & & & -{\sigma }_N{\mathbf{N}}_{N,N-1}& {\mathbf{D}}_{N,N-1}& {\sigma }_N{\mathbf{N}}_{N,N}\end{array}\right]\end{array}$$

$$[\mathbf{x}{]}_{2l-1}=a_l$$

$$[\mathbf{x}{]}_{2l}=b_l$$

$$[\mathbf{b}{]}_{2k}={\int }_{{\mu }_k}\left({\sigma }_{k+1}^{-1}{v}_{k+1}-{\sigma }_k^{-1}{v}_k\right)P_{0k}dr$$

$$[\mathbf{b}{]}_{2k-1}={\int }_{t_k}\left({\partial }_n{v}_{k+1}-{\partial }_n{v}_k\right)P_{1k}dr.$$

$$S_{ii}{N}_{ii}=\frac{1}{4}I-D_{ii}^2$$

$$N_{ii}{S}_{ii}=\frac{1}{4}I-D_{ii}^{\ast 2}$$

$$D_{ii}{S}_{ii}=D_{ii}^{\ast }{S}_{ii}$$

$$D_{ii}^{\ast }{N}_{ii}=N_{ii}{D}_{ii}^{\ast }$$

$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}Q=\\ \begin{bmatrix}\frac{1}{\sqrt{max(\sigma_{1},\sigma_{2})}}&0&0&0&\cdots&0\\ 0&\sqrt{min(\sigma_{1},\sigma_{2})}&0&0&\cdots&0\\ 0&0&\frac{1}{\sqrt{max(\sigma_{2},\sigma_{3})}}&0&\cdots&0\\ 0&0&0&\sqrt{min(\sigma_{2},\sigma_{3})}&\cdots&0\\ \vdots&\vdots&\vdots&\vdots&\ddots&\cdots\\ 0&0&0&0&0&\frac{1}{\sqrt{\sigma_{N}}}\\ \end{bmatrix}}$$

$$\begin{array}{rl}& Z_q=QZQ=\\ & \left[\begin{array}{cccccccc}(1+R_{11})N_{11}& -2\sqrt{R_{11}}{D}_{11}^{\ast }& -P_2{N}_{12}& \sqrt{R_{21}}{D}_{12}^{\ast }& & & & \\ -2\sqrt{R_{11}}{D}_{11}& (1+R_{11})S_{11}& \sqrt{R_{12}}{D}_{12}& -P_2^{\ast }{S}_{12}& & & & \\ -P_2{N}_{21}& \sqrt{R_{12}}{D}_{21}^{\ast }& (1+R_{22})N_{22}& -2\sqrt{R_{22}}{D}_{22}^{\ast }& \cdots & & & \\ \sqrt{R_{21}}{D}_{21}& -P_2^{\ast }{S}_{21}& -2\sqrt{R_{22}}{D}_{22}& (1+R_{22})S_{22}& \cdots & & & \\ & & ⋮& ⋮& \ddots & & & \\ & & & & & (1+R_{N-1,N-1})N_{N-1,N-1}& -2\sqrt{R_{N-1,N-1}}{D}_{N-1,N-1}^{\ast }& -P_N{N}_{N-1,N}\\ & & & & & -2\sqrt{R_{N-1,N-1}}{D}_{N-1,N-1}& (1+R_{N-1,N-1})S_{N-1,N-1}& \sqrt{R_{N-1,N}}{D}_{N-1,N}\\ & & & & & -P_N{N}_{N,N-1}& \sqrt{R_{N-1,N}}{D}_{N,N-1}^{\ast }& N_{N,N}\end{array}\right]\end{array}$$

$$R_{ij}=\frac{\min ({\sigma }_i,{\sigma }_{i+1})}{\max ({\sigma }_j,{\sigma }_{j+1})},P_i=\frac{{\sigma }_i}{\sqrt{\max ({\sigma }_{i-1},{\sigma }_i)\max ({\sigma }_i,{\sigma }_{i+1})}},P_i^{\ast }=\frac{\sqrt{\min ({\sigma }_{i-1},{\sigma }_i)\min ({\sigma }_i,{\sigma }_{i+1})}}{{\sigma }_i}$$

$$\begin{gathered} Z_q= \\ \left[\begin{array}{cccccccc}{c}_{11}{N}_{11}& c_{12}{D}_{11}^{\ast }& c_{13}{N}_{12}& c_{14}{D}_{12}^{\ast }& & & & \\ c_{21}{D}_{11}& c_{22}{S}_{11}& c_{23}{D}_{12}& c_{24}{S}_{12}& & & & \\ c_{31}{N}_{21}& c_{32}{D}_{21}^{\ast }& c_{33}{N}_{22}& c_{34}{D}_{22}^{\ast }& \cdots & & & \\ c_{41}{D}_{21}& c_{42}{S}_{21}& c_{43}{D}_{22}& c_{44}{S}_{22}& \cdots & & & \\ & & ⋮& ⋮& \ddots & & & \\ & & & & & c_{2N-3,2N-3}{N}_{N-1,N-1}& c_{2N-3,2N-2}{D}_{N-1,N-1}^{\ast }& c_{2N-3,2N-1}{N}_{N-1,N}\\ & & & & & c_{2N-2,2N-3}{D}_{N-1,N-1}& c_{2N-2,2N-2}{S}_{N-1,N-1}& c_{2N-2,2N-1}{D}_{N-1,N}\\ & & & & & c_{2N-1,2N-3}{N}_{N,N-1}& c_{2N-1,2N-2}{D}_{N,N-1}& c_{2N-1,2N-1}{N}_{N,N}\end{array}\right] \end{gathered}$$

$$c_{ij}=c_{ji},c_{2i-1,2i-1}=c_{2i,2i}$$

$$C_q=\left[\begin{array}{cccccccc}{c}_{11}{S}_{11}& -c_{12}{D}_{11}& c_{13}{S}_{12}& c_{14}{D}_{12}& & & & \\ -c_{21}{D}_{11}^{\ast }& c_{22}{N}_{11}& c_{23}{D}_{12}^{\ast }& c_{24}{N}_{12}& & & & \\ c_{31}{S}_{21}& c_{32}{D}_{21}& c_{33}{S}_{22}& -c_{34}{D}_{22}& \cdots & & & \\ c_{41}{D}_{21}^{\ast }& c_{42}{N}_{21}& -c_{43}{D}_{22}^{\ast }& c_{44}{N}_{22}& \cdots & & & \\ & & ⋮& ⋮& \ddots & & & \\ & & & & & c_{N-2,N-2}{S}_{N-1,N-1}& -c_{N-2,N-1}{D}_{N-1,N-1}& c_{N-2,N}{S}_{N-1,N}\\ & & & & & -c_{N-1,N-2}{D}_{N-1,N-1}^{\ast }& c_{N-1,N-1}{N}_{N-1,N-1}& c_{N-1,N}{D}_{N-1,N}^{\ast }\\ & & & & & c_{N,N-2}{S}_{N,N-1}& c_{N,N-1}{D}_{N,N-1}& c_{N,N}{S}_{N,N}\end{array}\right]$$

$$C_q{Z}_q=\left[\begin{array}{cccc}{c}_{11}^2{S}_{11}{N}_{11}+K_{11}& K_{12}& K_{13}& \cdots \\ K_{21}& c_{22}^2{N}_{11}{S}_{11}+K_{22}& K_{23}& \cdots \\ K_{31}& K_{32}& c_{33}^2{S}_{22}{N}_{22}+K_{33}& \cdots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right]$$

$$\begin{array}{rl}{K}_{2n-1,2m-1}=& \sum _{i=n-1}^{n+1}{\chi }_m(i)(c_{2n-1,2i-1}{c}_{2i-1,2m-1}{S}_{ni}{N}_{im}+(-1{)}^{{\delta }_{ni}}{\chi }_i{c}_{2n-1,2i}{c}_{2i,2m-1}{D}_{ni}{D}_{im})\\ & -c_{2n-1,2m-1}{S}_{nm}{N}_{nm}{\delta }_{nm}\end{array}$$

$$\begin{array}{rl}{K}_{2n,2m}=& \sum _{i=n-1}^{n+1}{\chi }_m(i)((-1{)}^{{\delta }_{ni}}{c}_{2n,2i-1}{c}_{2i-1,2m}{D}_{ni}^{\ast }{D}_{im}^{\ast }+{\chi }_i{c}_{2n,2i}{c}_{2i,2m}{N}_{ni}{S}_{im})\\ & -c_{2n,2m}{N}_{nm}{S}_{nm}{\delta }_{nm}\end{array}$$

$$K_{2n-1,2m}=\sum _{i=n-1}^{n+1}{\chi }_m(i)(c_{2n-1,2i-1}{c}_{2i-1,2m}{S}_{ni}{D}_{im}^{\ast }+(-1{)}^{{\delta }_{ni}}{\chi }_i{c}_{2n-1,2i}{c}_{2i,2m}{D}_{ni}{S}_{im})$$

$$K_{2n,2m-1}=\sum _{i=n-1}^{n+1}{\chi }_m(i)((-1{)}^{{\delta }_{ni}}{c}_{2n,2i-1}{c}_{2i-1,2m-1}{D}_{ni}^{\ast }{N}_{im}+{\chi }_i{c}_{2n,2i}{c}_{2i,2m-1}{N}_{ni}{D}_{im})$$

$$\chi_{m}(i)=\left\{\begin{array}[]{ll}1&\mbox{if }|i-m|<2\\ 0&otherwise\end{array},\right.$$

$$\chi_{i}=\left\{\begin{array}[]{ll}1&\mbox{if }i<N\\ 0&otherwise\end{array},\right.$$

$$\delta_{ij}=\left\{\begin{array}[]{ll}1&\mbox{if }i=j\\ 0&\mbox{otherwise}\end{array}.\right.$$

Full text

a Computational Electromagnetics Research Laboratory, IMT Atlantique, Brest, France

b Politecnico di Torino, Turin, Italy

K****eywords EEG  $\cdot$ BEM  $\cdot$ Calderon Solvers

1 Introduction

Functional brain imaging based on high-resolution scalp Electroencephalographies (EEGs) is characterized by a high temporal resolution and, as such, it provides an unmatched overview on the underlying brain activity [1, 2, 3, 4]. This technique relies on the key task, referred to as the EEG inverse problem, of recovering the brain electric current sources responsible for a measured potential at the EEG scalp electrodes [5, 6]. The EEG inverse problem requires multiple solutions of the EEG forward problem, i.e. the computation of the scalp potential starting from the source currents [7, 8]. It has been widely studied and reported that the accuracy of EEG forward problem solvers has a direct impact on EEG inverse solution procedures [9, 10, 11, 12]. For this reason, any advancement of the state of the art in EEG forward solution technologies will have a direct impact on the overall high-resolution EEG imaging process.

When realistic head models are used [11, 13, 14, 15, 16, 17], the solution of the EEG forward problem can only be obtained numerically. Classical strategies to obtain this numerical solution are the Finite Element Method (FEM), the Finite Difference Method (FDM) or the Boundary Element Method (BEM) [18]. Both FDM and FEM leverage on a volume discretization of the considered head model. This allows these methods to account for the inhomogeneity and anisotropy of the head’s conductivity at the cost, however, of a higher computational demand. Previous works have shown that by using transfer matrices the computation time of the FEM formulations can be reduced [19, 20]. The use of these transfer matrices in FEM formulations yield similar computational times as BEM formulations for comparable accuracies. [21].

When the conductivity of the head is modelled as piecewise homogeneous, BEM can be easily used to compute the solution of the EEG forward problem. In other words, the main limitation of the BEM formulation resides in its inability to model anisotropies. However, this method has the advantage that it requires only the discretization of the interface between regions with different conductivities [18, 22]. Several studies (for example [12, 23, 24]) focused on the impact of the head model simplifications in recovering the electric brain sources from the measurement of scalp potential. In particular, when computing the EEG forward problem, [25, 26, 27] have shown the importance of modelling correctly the skull anisotropy. However, when the anisotropic conductivity values are not known, it can be preferable to model this region as isotropic, as explained in [25]. Moreover, the anisotropic conductivity of the skull is due to its layered structure, a cancellous bone between two compact bones. This means that when those three layers are available, the skull can accurately be modelled with three isotropic layers instead of one anisotropic layer as [28] shows.

The relevant computational savings which the use of BEM strategies can lead to, explain the attention the technique has received by the community, resulting in a continuous series of advances [18, 29, 30, 31, 32, 33]. Among them a method, published in [34] and referred to as the “symmetric formulation”, became quite popular and impacted several EEG based imaging tools [35, 36, 37, 38]. The peculiarity of the BEM method proposed in [34] is the quite higher level of accuracy that it can achieve when compared to previously existing schemes. However these beneficial properties are obtained at the cost of using a first kind formulation (while the majority of standard strategies relies on second kind formulations). The computational consequence of this fact is that, when the “symmetric formulation” is discretized to be solved numerically, the condition number of the resulting BEM matrix (the ratio of the largest over the smallest singular value of the matrix) will grow as a function of the discretization density (the number of boundary elements used to discretize the structure) [39]. Similarly, a condition number growth is observed in the symmetric formulation also when the conductivity contrast between two regions of the head is increased (a case of practical interest given that the conductivity of the skull is often modeled with a much smaller value with respect to the conductivity of the brain [11, 40, 41, 42]). In several cases, especially when handling models issued of high resolution Magnetic Resonance Imaging (MRI) [43], the solution of the EEG forward problem is obtained iteratively [18, 44]. A low and stable condition number is desirable since, on the one hand, the number of iterations of an iterative solver is growing with the condition numbers [45] and, on the other hand, the condition number controls the amplification in the solution of any initial error in the sources [45]. In other words, the higher the condition number, the longer the time needed to compute the solution, and the less correct the solution will be.

The purpose of this work is to address the ill-conditioning problems of the symmetric formulation. Given the favor that the formulation has found in the community and the fact that it is already implemented in several neuroimaging packages, a particular attention will be devoted to develop a solution strategy that will be conservative, in the sense that will not require the change of previous implementations of the symmetric formulation but will just require the addition of some extra steps to it. This will be achieved by developing a purely multiplicative preconditioner based on Calderon formulas, i.e. we will design a preconditioning matrix that is spectrally equivalent to the inverse of the symmetric formulation. After left multiplication of this matrix with the symmetric formulation matrix, the resulting linear system will, on the one hand, keep the accuracy the symmetric formulation is well known for and, on the other hand, will provide a stable condition number both when the mesh is refined and when the conductivity contrast between two adjacent domains increases. The reader should notice that Calderón strategies have been successfully applied to the regularization of other integral operators in the context of full-wave vector electromagnetic problems scattered by metallic [46, 47, 48, 49, 50, 51, 52, 53] and penetrable [54, 55, 56, 57] objects. The regularization presented here, however, does not automatically follow from none of the strategies above due to the peculiar nature, both in terms of frequency and operator structure/properties, of the symmetric formulation under consideration in this work. It should also be noted that the problem of simulating high contrast dielectric materials has been addressed for high frequency formulations with Calderon strategies both in [55] and [57]. Both of these two clever approaches, designed for the vector case, are not applicable in our scalar scenario where the spectral high-contrast scaling expresses itself as a block scaling issue due to the decoupling of electric and magnetic quantities characterizing a static problem.

This paper is organized as follows: Section 2 provides the reader with some necessary background material and notation used in the following developments. Section 3 presents the new Calderon preconditioner proposed in this work, while Section 4 focuses on its discretization and on the solution of the preconditioned symmetric formulation system. Section 5 complements the paper’s theoretical developments with numerical results which will show the efficiency and effectiveness of the new approach. Partial results from this work has been presented in the conference contribution [58].

2 Background on the EEG Forward Problem

This section will briefly review the relevant formulations, currently available in literature, used to solve the EEG forward problem. The treatment will be synthetic and for the sole purpose of setting up the notation. The reader interested in a more profuse treatment should refer, for example to [34, 59] and to references therein.

2.1 The EEG Problem

Let $\Omega=\bigcup_{i=1}^{N}\Omega_{i}$ be a nested domain with Lipchitz boundaries $\partial\Omega_{i}=$ $\left(\bar{\Omega}_{i-1}\bigcap\bar{\Omega}_{i}\right)\bigcup\left(\bar{\Omega}_{i}\bigcap\bar{\Omega}_{i+1}\right)$ as in Fig. 1. We denote with $\boldsymbol{n}_{i}$ the outward going normal to the surface $\Gamma_{i}$, where $\Gamma_{i}=\bar{\Omega}_{i}\bigcap\bar{\Omega}_{i+1}$.

Solving the EEG forward problem amounts to computing the potential $V$ at given electrodes’ positions when the active brain current sources are known. Under quasi-static assumptions and isotropic conductivity, the EEG forward problem reads [60]:

[TABLE]

where $\Delta=\nabla\cdot\nabla$ is the Laplace operator, $\sigma$ is the conductivity and $\boldsymbol{j}$ the current sources. The conductivity is assumed to be piecewise isotropic and homogeneous: in $\Omega_{i}$, $\sigma=\sigma_{i}$. In the exterior domain, the conductivity is assumed to be 0. The current sources, as it is customary in literature [61], are assumed to be dipolar in nature. Hence, denoting with $f_{i}=\nabla\cdot\boldsymbol{j}$ the electric source in $\Omega_{i}$, we have $f_{i}=q_{i}\cdot\nabla\delta_{r_{i}}$ with $q_{i}$ the electric dipole moment and $\boldsymbol{r}_{i}$ its position. Furthermore, the symbol $[g]_{i}=g^{-}-g^{+}$, will refer to the jump of the function $g$ at the interface $\Gamma_{i}$, with $g^{\mp}$ the inner and outer trace of $g$ at $\Gamma_{i}$ respectively. Then, the solvability of (1) is assured under the following boundary conditions [60]:

[TABLE]

that enforce the continuity of the potential and the current between the different layers of the domain $\Omega$.

The Green’s function associated with (1) reads [39]

[TABLE]

for which we can derive Green’s representation theorem using the integral operators [39]

[TABLE]

In the above equation and $p.v.$ and $f.p$ stand for Cauchy principal value and Hadamard finite part respectively. In the following, we will denote with $L_{ij}$ the operator $L$ when $r\in\Gamma_{i}$ and $r^{\prime}\in\Gamma_{j}$ with $L=D,S,N\text{ or }D^{*}$.

2.2 The Symmetric Formulation for the EEG Forward Problem

Several BEM formulations have been proposed to solve the EEG forward problem [18, 34]. Among them, the symmetric formulation [34] is quite popular and known for providing high levels of accuracy [36]. In solving the EEG forward problem, an efficient strategy is to build the unknown potential $V$ starting from two functions, a function $u$ harmonic in $R^{3}$ and a function $v$ that takes into account the source term. The starting point of the symmetric formulation is to build $u_{i}$ in each domain such that $u_{i}=V-v_{i}/\sigma_{i}$ in $\Omega_{i}$ and $u=-v_{i}/\sigma_{i}$ in $R^{3}\setminus\Omega_{i}$, with $v_{i}$ the solution of (1) in an unbounded medium: $v_{i}(r)=\int_{\Omega_{i}}f_{i}(r^{\prime})G(r,r^{\prime})dr^{\prime}$. In this fashion, $u_{i}$ is harmonic in $R^{3}\setminus\partial\bar{\Omega_{i}}=\Gamma_{i-1}\cup\Gamma_{i}$. Using the boundary conditions (2a) and (2b) as well as the representation theorem [39], two integral equations for the potential and its derivative can be obtained [34]. They read:

[TABLE]

[TABLE]

with $V_{i}$ the potential on the surface $\Gamma_{i}$ and $p_{i}=\sigma_{i}\left[\boldsymbol{n}\cdot\nabla V\right]_{i}$. Equations (5) and (6) are obtained by applying the boundary conditions on the surface $\Gamma_{i}$. In a nested domain, it only involves the computation of the operators with functions defined in the surrounding surfaces $\Gamma_{i-1}$, $\Gamma_{i}$ and $\Gamma_{i+1}$.To clarify the ideas, equations (5) and (6) have been rewritten in matrix form in equation (7). For a more detailed explanation on the symmetric formulation, the reader is referred to [34].

[TABLE]

2.3 Discretization of the Operators

The numerical solution of an integral equation is often obtained by using a Boundary Element Method (BEM). Following a well-established strategy, the surface $\Gamma$ is tessellated into $N_{t}$ triangular cells $t_{k}$ of area $A_{k}$ and average length $h$. The set of vertices of the tessellation will be denoted by $\{v_{k}\}_{k=1}^{N_{v}}$. Cells and vertices will form a mesh denoted by $M_{\Gamma}$. The number of triangles (respectively, vertices) of the surface $\Gamma_{i}$ will be denoted $N_{t_{i}}$ (respectively, $N_{v_{i}}$ ). To discretize the unknowns and to test the equations the following standard basis functions will be used. The piecewise constant functions in $P_{0}$ are defined such that, $P_{0k}=1/A_{k}$ in $t_{k}$ and 0 elsewhere. The piecewise linear functions are the set $P_{1}=span\{P_{1k}\}_{k=1}^{N_{v}}$. The support of $P_{1k}$, denoted by $\mu_{P_{1k}}$, is the set of triangles around $v_{k}$ such that $P_{1k}=1$ in $v_{k}$ and [math] on all other vertices. $P_{0}$ and $P_{1}$ functions are shown Fig. 3a and 3b respectively.

Following [34], on the surface $\Gamma_{i}$ we respectively discretize the two unknowns $V_{i}$ and $p_{i}$ with $P_{1}$ and $P_{0}$ basis functions such that $V_{i}=\sum_{l}^{N_{v}}a_{k}P_{1l}$ and $p_{i}=\sum_{l}^{N_{t}}b_{l}P_{0l}$. In order to obtain the system matrix, the integral equations (5) and (6) are tested with $P_{0}$ and $P_{1}$ basis functions respectively. The operator matrices arising are then given by

[TABLE]

Consequently, the system to be solved reads $\mathbf{Z}\mathbf{x}=\mathbf{b}$ with $\mathbf{Z}$ given by

[TABLE]

where

[TABLE]

and where

[TABLE]

In the following, the operator associated with the matrix $\mathbf{Z}$ (obtained by replacing in (9) $\mathbf{D}$, $\mathbf{S}$, $\mathbf{N}$, and $\mathbf{D^{*}}$ with ${D}$, ${S}$, ${N}$, and ${D^{*}}$) will be denoted by $Z$.

3 A Calderon Preconditioner for the Symmetric Formulation

The high accuracy of the symmetric BEM formulation [34] has made of it a very popular tool for solving the EEG forward problem. However its system matrix suffers from ill-conditioning that can lead to the non-convergence of the employed iterative solver used to compute the solution [45], which has to be used for high refined meshes where the direct methods are impractical. Indeed, the operator $S$ is compact [39]. This means that its spectrum will accumulate at zero when the mesh is refined and it will therefore have a condition number increasing inversely proportional to the average mesh length $h$. Moreover, the hypersingular operator $N$ is an unbounded operator [39]. This implies that its condition number will also grow with $1/h$. Since these operators are the diagonal blocks of the matrix $\mathbf{Z}$ in (9) and the off-diagonal blocks of the matrix are smoothers [62], it follows that the overall conditioning of $\mathbf{Z}$ will increase when the mesh discretization will increase ($h\to 0$).

By leveraging on the Calderon identities, it is possible to build a preconditioner for the system matrix $\mathbf{Z}$. The rationale behind our strategy can be understood by considering the continuous operators first. The Calderon identities that for our approach read [63]

[TABLE]

where the symbol $I$ stands for the identity operator. The spectral analysis of (12a) and (12b) shows that the operators $S_{ii}N_{ii}$ and $N_{ii}S_{ii}$ are well conditioned. Indeed, given that $D$ and $D^{*}$ are compact operators, then $D_{ii}^{2}$ and $D^{*2}_{ii}$ are also compact operators as a product of two compact operators. Then the spectrum of $D_{ii}^{2}$ and $D^{*2}_{ii}$ is bounded above and accumulates at zero. However, the presence of the identity operator in (12a) and (12b) guarantees that the spectrum of the operators $S_{ii}N_{ii}$ and $N_{ii}S_{ii}$ will be approaching $1/4$. In other words, $S_{ii}N_{ii}$ and $N_{ii}S_{ii}$ are second kind operators whose spectrum accumulates at $1/4$. This property can be exploited to build a left preconditioner for the symmetric operator $Z$.

Before introducing the Calderon preconditioner, we present a regularization matrix for the coefficients of the symmetric operator, which is necessary for the stability of the matrix condition number with respect to the conductivity ratio and for the correct behaviour of the subsequent mesh refinement preconditioner. Indeed, the symmetric operator is unstable when there is a high conductivity contrast between two adjacent domains ,$\Omega_{i}$ and $\Omega_{j}$, due to the conductivity factors in the diagonal blocks of the matrix. These conductivity factors are given by $a_{2l-1,2l-1}=\sigma_{l}+\sigma_{l+1}$ and $a_{2l,2l}=\sigma_{l}^{-1}+\sigma_{l+1}^{-1}$. As a consequence, in the case of high conductivity contrast between two adjacent domains, i.e. asymptotically, when $\frac{\sigma_{i}}{\sigma_{j}}\rightarrow\infty$, that is $\left(\sigma_{i}+\sigma_{j}\right)\rightarrow\infty$ when $\sigma_{i}\rightarrow\infty$ or $\left(\sigma_{i}^{-1}+\sigma_{j}^{-1}\right)\rightarrow\infty$ when $\sigma_{j}\rightarrow 0$, the condition number of the system matrix will grow as a function of the conductivity ratio $CR_{ij}=\frac{\max(\sigma_{i},\sigma_{j})}{\min(\sigma_{i},\sigma_{j})}$. Because of high conductivity contrast between the brain and the skull [40, 41], this undesirable situation is likely to appear when solving the EEG forward problem. To solve this problem, we will rescale with respect to the conductivity the symmetric operator $Z$ using a diagonal operator $Q$ given by

[TABLE]

We then define a rescaled symmetric operator $Z_{q}$ as

[TABLE]

with

[TABLE]

and $R_{ij},P_{i},P^{*}_{i}\leq 1$. The asymptotic behaviour of those coefficients is presented in table 1 , from which we can see that when $CR_{ij}\to\infty$, none of the blocks of $Z_{q}$ tends to infinity. Furthermore, the diagonal blocks are bounded between 1 and 2. Therefore the instability with respect to high conductivity ratio of the symmetric operator have been solved.

In order to simplify the notation for the following part, we write the regularized symmetric operator as

[TABLE]

where the coefficients $c_{ij}$ have the properties

[TABLE]

Now we present the Calderon preconditioner for the regularized symmetric operator $Z_{q}$ which is based on the Calderon identities (12). We denote such preconditioning operator by $C_{q}$. Its definition is given in (18), where the constant coefficients $c_{ij}$ are the same as in (16).

[TABLE]

We highlight the minus signs present in the blocks $D^{*}_{ii}$, which are necessary for the compactness of the off diagonal blocks when the preconditioner is applied. A minus sign is added in the $D_{ii}$ blocks as well for symmetry.

Then, as desired, $C_{q}Z_{q}$ is a block operator exhibiting the Calderon identities (12a) and (12b) in its diagonal up to the multiplicative factor $c_{ii}^{2}$ as can be seen in (19).

[TABLE]

The terms denoted with $K_{ij}$, contain linear combinations of the compact operators $D_{ij}D_{jk}$, $S_{ij}N_{jk}$, $S_{ij}D^{*}_{jk}$, $D_{ij}S_{jk}$, $D^{*}_{ij}N_{jk}$, $N_{ij}D_{jk}$, $N_{ij}S_{jk}$ and $D^{*}_{ij}D^{*}_{jk}$. They read

[TABLE]

where the symbols $\chi_{i}$, $\chi_{m}(i)$ are given by

[TABLE]

[TABLE]

and where $\delta_{nm}$ is the Kronecker’s delta

[TABLE]

As shown previously, the terms $S_{ii}N_{ii}$ and $N_{ii}S_{ii}$ are well conditioned second kind operator matrices, as a consequence the terms $c_{ii}S_{ii}N_{ii}$ $c_{ii}N_{ii}S_{ii}$ will also be well conditioned with respect to the mesh parameter $h$ since $1<c_{ii}<2$.

In order to show the compactness of the operators $K_{ij}$, we analyze each operator product in them. First, in the diagonal blocks $K_{2n-1,2n-1}$ and $K_{2n,2n}$ we have the following products: $S_{ni}N_{in}$, $N_{ni}S_{in}$, $D^{*}_{ni}D^{*}_{in}$, $D_{ni}D_{in}$. The operators $D$ and $D^{*}$ are compact, therefore the products $D_{ni}D_{in}$ and $D^{*}_{ni}D^{*}_{in}$ are compact as well [63]. The products $S_{ni}N_{in}$, $N_{ni}S_{in}$ are present when $i\neq n$, then the operators $N_{in}$ and $N_{ni}$ have a regular kernel due to the lack of singularity. Hence, the products with the compact operators $S_{ni}$ and $S_{in}$ are also compact operators [39, 64, 65]. Overall, in the diagonal block we have the sum of compact operators. In the off diagonal blocks $K_{2n-1,2m-1}$, $K_{2n,2m}$, when $n\neq m$, we have the products $S_{ni}N_{im}$, $D_{ni}D_{im}$, $D^{*}_{ni}D^{*}_{im}$, $N_{ni}S_{im}$. As stated before, the product of the $D^{*}$ and $D$ are compact. Moreover, in the products $S_{ni}N_{im}$ and $N_{ni}S_{im}$ at least one of the operators has a regular kernel, since $n\neq m$. Hence, these products yield compact operators [64, 65]. Finally, for the off diagonal blocks $K_{2n-1,2m}$, $K_{2n,2m-1}$ we have the following products: $S_{ni}D^{*}_{im}$ , $D_{ni}S_{im}$, $D^{*}_{ni}N_{im}$, $N_{ni}D_{im}$. The products $S_{ni}D^{*}_{im}$ and $D_{ni}S_{im}$ are clearly compact, since all the operators involved are compact. The only not-compact terms are present in the products $D^{*}_{ni}N_{im}$ and $N_{ni}D_{im}$ when $n=i=m$. However, recalling the Calderon identity (12d) and the properties (17), we have

[TABLE]

which is true for the continues operators. When the operators are discretized, this cancellation is not exact, however the residue tends to zero.

Then, we can write

[TABLE]

with

[TABLE]

and $B$ such that $[B]_{ij}=K_{ij}$, $C_{q}Z_{q}$ can be seen as the sum of the well conditioned matrix $A$ and a compact perturbation $B$ (as the operators $K_{ij}$ are compact operators). We can therefore expect the operator $C_{q}Z_{q}$ to be well conditioned with respect to the mesh parameter.

For now on, we will refer to the product

[TABLE]

as the Calderon-Symmetric operator, which is well conditioned for both mesh refinement and conductivity ratio between adjacent surfaces.

4 Discretization of the Calderon Preconditioner and Solution of the Preconditioned Symmetric Formulation

In order to solve the preconditioned symmetric integral equation, the proposed multiplicative preconditioner $C_{q}$ has to be discretized. This discretization should be carried out with care. In fact, the preconditioned operator in (25) will contain operator products which will not directly translate into matrix products in the general case. A suitable choice of basis functions should be made to guarantee that this could instead be the case here. To fix the ideas, we could consider the discretization of the operator product $N_{11}S_{11}$ appearing in the top-left block of $Z_{c}$. The matrix $\mathbf{S_{11}}$ is obtained by using test and trial functions in $P_{0}$ while the matrix $\mathbf{N_{11}}$ is obtained by using test and trial functions in $P_{1}$. Yet the number of vertices $N_{t}$ and the number of cells that defines the dimensions of the space $P_{1}$ and $P_{0}$ are different. As a consequence, the blocks $\mathbf{N_{11}}$ and $\mathbf{S_{11}}$ do not have compatible shapes and cannot be multiplied. Furthermore, the basis functions used for discretizing $N_{11}$ and $S_{11}$ must satisfy appropriate inf-sup conditions with respect to the duality pairing $\langle v,w\rangle:H^{1/2}(\Gamma_{1})\times H^{-1/2}(\Gamma_{1})\rightarrow\mathbb{R}$ (the reader should refer to [66] for further technical details on this topic). This condition enables to get a stable condition number for the Gram matrices, necessarily present in the discretized system as they orthonormalize the two chosen basis and testing functions sets. To properly take care of this fact, we propose to discretize the preconditioner $C_{q}$ on the dual mesh $M_{\Omega}^{*}$ of the standard mesh $M_{\Omega}$ and to leverage on the dual basis functions introduced in [67] on such a mesh. In the dual mesh, each vertex corresponds to a cell of the standard mesh and vice-versa. This means that in $M_{\Omega}^{*}$ we can build a discrete space in $H^{1/2}(\Gamma)$ which has the same dimension as the discrete space associated with $H^{-1/2}(\Gamma)$ in $M_{\Omega}$ and vice-versa. Moreover, the dual basis functions introduced by [67] abide by the inf-sup conditions required to obtain stable discrete products [66]. As a consequence, the discretization of the preconditioner operator $C_{q}$ by using these basis functions enables to perform the matrix multiplication associated with the operator multiplication $C_{q}Z_{q}$ in such a way that the spectral bounds holding for the continuous operator products will translate in well-conditioned matrix products.

The dual mesh $M_{\Omega}^{*}$ can be obtained by barycentric refinement of the standard mesh $M_{\Omega}$ by dividing the triangles $t_{k}$ into six smaller triangles $t_{bk}$ whose edges are built by tracing the medians of the standard triangles $t_{k}$ [67]. The cells $c_{k}$ of $M_{\Omega}^{*}$ are defined as the set of triangle $t_{bk}$ sharing a common vertex $v_{k}$ of $M_{\Omega}$. The vertices of $M_{\Omega}^{*}$ are the barycenters $b_{k}$ of the triangles $t_{k}$ in $M_{\Omega}$. The reader should refer to Fig.  2 for an example of such a refinement.

In $M_{\Omega}^{*}$, we define the dual piecewise linear functions $\tilde{P_{1}}=span\{\tilde{P}_{1_{k}}\}_{k=1}^{N_{t}}$ obtained with a linear combination of $P_{1}$ functions built on the barycentric refined mesh [67]. A dual piecewise linear function is shown Fig. 3d. The coefficients of the linear combination are shown Fig. 2. The support of $\tilde{P}_{1_{i}}$ is denoted $\mu_{\tilde{P}_{1_{i}}}$. The dual piecewise constant functions in $\tilde{P_{0}}$, denoted $\tilde{P}_{0k}$, are the constant functions equal to $1/A_{c_{k}}$ on the cell $c_{k}$, whose area is $A_{c_{k}}$, of $M_{\Omega}^{*}$ and equal to zero elsewhere. A dual piecewise constant function is shown in Fig. 3c. An extended explanation of the dual mesh and dual basis functions is given by [67].

Then the operators matrices used to build $\mathbf{C_{q}}$ are given by:

[TABLE]

and the discretized preconditioner $\mathbf{\tilde{C}_{q}}$ is made explicit as

[TABLE]

The final preconditioner is then obtained by introducing the necessary rescaling to obtain a uniform conditioning with respect to the conductivity profiles (according to the analysis of the previous section).

In order to perform the multiplication of matrices with two different discretizations, a Gram matrix $\mathbf{G}$ to link them is necessary. This Gram matrix is computed by taking the scalar product between the trial functions of one operator and the test functions of the other operator. Hence, the computation of the Gram matrix does not require the evaluation of any operator. Additionally, it is almost diagonal, therefore the computational cost is very low. This matrix is obtained as

[TABLE]

Finally, the discretization of the Calderon symmetric operator is given by

[TABLE]

where $\mathbf{Z_{q}}=\mathbf{Q}\mathbf{Z_{q}}\mathbf{Q}$. The solution of the preconditioned symmetric formulation is then obtained by solving the following system $\mathbf{Z}_{c}\mathbf{y}=\mathbf{C}_{q}\mathbf{G}^{-1}\mathbf{Q}\mathbf{b}$ and $\mathbf{x}$ is obtained with $\mathbf{x}=\mathbf{Q}\mathbf{y}$.

Summarizing, the Calderon preconditioning strategy is multiplicative in nature. Its aim is to build a preconditioning operator spectrally equivalent to the inverse of the original operator. Thus, once this operator is built, multiplying the ill-conditioned operator with it yields an operator spectrally equivalent to an identity. The preconditioning operator is built on a dual mesh in order to allow matrix multiplication and stability. Moreover, regularization matrices are added in order to get a condition number independent of the conductivity ratio. In a nutshell, the steps are:

Compute the standard symmetric system matrix $\mathbf{Z}$; 2. 2.

Compute the Calderon preconditioning matrix $\mathbf{\tilde{C}_{q}}$ on the dual mesh; 3. 3.

Compute the Gram matrices linking the dual and standard discretization, known as Gram matrices $\mathbf{G}$; 4. 4.

Normalize the operator $\mathbf{Z}$ with the regularization matrices $\mathbf{Q}$; 5. 5.

Perform the multiplication $\mathbf{Z}_{c}=\mathbf{\tilde{C}_{q}}\mathbf{G}^{-1}\mathbf{Z_{q}}$; 6. 6.

The right hand side b must be modified accordingly : compute $\mathbf{b}_{c}=\mathbf{\tilde{C}_{q}}\mathbf{G}^{-1}\mathbf{Q}\mathbf{b}$; 7. 7.

Solve the system $\mathbf{Z}_{c}\mathbf{y}=\mathbf{b}_{c}$; 8. 8.

Get the solution using $\mathbf{x}=\mathbf{Q}\mathbf{y}$.

5 Numerical Results

The new Calderon regularized symmetric formulation proposed in this work has been first tested on the canonical scenario of three homogeneous and concentric spheres of radii 0.8, 0.9, and 1 in normalized units respectively. Such model is shown in Fig. 4a. Indeed, in the case of homogeneous nested spheres, an analytical solution is available as a reference [68, 69], this solution will be denoted with $V_{ref}$. In these simulations, a single dipole source is placed in (0, 0, 0.5) with a dipole moment of (0, 0, 1) and the normalized conductivity of the layers is chosen to be 1, 0.0125 and 1 starting from the inner domain, which are the conductivity values used in the standard symmetric formulation in [34]. As a complement to these results, and to validate the new formulation on a real case scenario, the new formulation has been tested also on a realistic head model obtained from MRI data , which is presented in Fig. 4b.

5.1 Assessments on accuracy and condition number

The first test conducted aimed at verifying that applying the proposed preconditioner to the symmetric formulation does not modify its accuracy. The assessment parameter is the relative error computed as $RE=\frac{||V_{num}-V_{ref}||}{||V_{ref}||}$ where $V_{num}$ refers to the numerical solution. In Fig. 5a and Fig. 5b it is shown that the Calderon preconditioned symmetric formulation and the non preconditioned symmetric formulation provide exactly the same accuracy for different mesh refinement levels and different conductivity ratios respectively. This means that the proposed preconditioner does not alter the accuracy of the initial formulation. These figures also confirm the higher level of accuracy that the symmetric formulation can reach with respect to two others existing BEM formulations, namely the adjoint double layer formulation and the double layer formulation.

In the second test we compared the condition numbers of the preconditioned and the non preconditioned symmetric formulation system matrices . The variation of the condition number with respect to the mesh refinement is shown in Fig. 6a. We see that the condition number of the symmetric formulation grows rapidly with the mesh refinement parameter $1/h$ while the condition number of the proposed formulation stays constant as expected by the theory. In Fig. 6b is presented the condition number when the conductivity ratio increases. Here, we can see that the new formulation has a condition number independent of the conductivity contrast.

The third test aimed at showing the efficiency in terms of number of iterations of the proposed preconditioner. For different spectral indices $1/h$, we present in Fig. 7a the number of iterations needed for a Conjugate Gradient Square (CGS) solver to reach a relative residual error of $10^{-6}$. In this figure, it is clear that only with the proposed Calderon preconditioner this number of iterations stays constant with the mesh refinement. The number of iterations for increasing conductivity ratio is shown in Fig. 7b. In this case, all the preconditioners lead to a constant number of iterations. Still, the best performance is given by the Calderon preconditioner.

5.2 Assessments on time

Since for small number of unknowns direct solvers can be used to solve the obtained system of equations, a preconditioner is only useful when the number of unknowns becomes too high for direct inversion, that is when it is necessary to use iterative solvers. Indeed, if $N$ denotes the number of unknowns, then the complexity of a direct solver is $O(N^{3})$ while the complexity of an iterative solver is $O(kN^{2})$, where $k$ is the number of iterations needed to reach a desired accuracy. However, for high number of unknowns, given that BEM matrices are dense, the time needed to compute the system matrix is also increasing. To deal with this issue and to show the benefits of the proposed Calderon preconditioner in a high number of unknowns context we coupled it with an Adaptive Cross Approximation (ACA) algorithm [70], that provides a compressed version of the system matrix.

The time necessary to compute the solution for different numbers of unknown is shown Fig.8. In this test, we compare the performance of direct and iterative solutions with or without the proposed Calderon preconditioner. It is clear that direct inversion (DI) has, as expected, a time complexity that increases rapidly with $N$. This prevents the use of this solver for very detailed models. The iterative solver used is the CGS algorithm. We see that without using any preconditioner the time necessary to obtain the solution with this iterative solver is also increasing with the number of unknowns even if this choice of solver is faster than DI. This is due to the fact the number of iterations increases with the mesh refinement parameter (due to the ill-conditioning of the matrices). However, employing the proposed Calderon preconditioner solves this issue.

The time for computing the dense and compressed Calderon-Symmetric operator is presented in Fig. 9. It can be seen that the use of the ACA yields in a linear time complexity

5.3 Assessments on a MRI-obtained head model

Finally, we seek to assess the performance of the proposed preconditioner in a realistic scenario that consists in computing a leadfield matrix using a head model obtained from MRI data. This matrix that provides the propagation model between known brain electric current sources and electrodes situated at the surface of the head of a patient is a key element in distributed inverse solution. For this purpose, we constructed a three layer mesh using [35]. These layers model the brain, the skull, and the scalp. They contain 11850, 11616 and 22948 triangular cells respectively. The potential generated on the scalp by a single dipole situated in the brain is presented Fig. 10a. This figure also shows the position of the 21 electrodes for which we computed the leadfield matrix. The error when the Calderon-Symmetric formulation is solved with the ACA compared with the original Symmetric formulation solved with direct inversion is displayed Fig. 10b. It can be seen that the error is never greater than $0.05\%$. To fill-in the leadfield matrix, we placed 1500 unitary dipole source in the brain layer, each having an orientation orthogonal to the brain surface. Using reciprocity [71], the forward model is then solved at each electrode position. We compared in Table 2, for four different cases, the time needed to compute the operator, the time needed to solve the forward system once, and the time needed to compute the full leadfield matrix. Hence, the total time needed to get the leadfield matrix is given by the sum of the computation for obtaining the operator and the leadfield matrix. It can be seen that even if the time necessary to compute the compressed Calderon-Symmetric operator is greater than the compressed Symmetric operator, the fast convergence of the new method allows to compute the complete leadfield matrix in 2.56 hours, that is almost 10 times faster than without the proposed preconditioner. This compensates largely the computation overload in computing the preconditioning operator.

6 Conclusion

In this work, leveraging on Calderon identities, we have proposed a Calderon preconditioned symmetric formulation for the EEG forward problem. When compared to the standard symmetric equation, the proposed formulation has the advantage of showing constant condition numbers both as a function of the mesh refinement and of the conductivity contrast. Employing this preconditioner does not degrade the accuracy of the symmetric formulation. This means that, for a given relative accuracy of the solution, the proposed formulation converges substantially faster than the standard one. Moreover, the integration of the proposed approach into existing symmetric formulation implementations is achieved only at the cost of computing, on a barycentric refined mesh, the preconditioning operators with already existing tools. Numerical results have substantiated the theoretical claims and have shown the practical impact of the newly proposed scheme.

Acknowledgement

This work was supported in part by CominLabs project SABRE under the reference ANR-10-LABX-07-01, in part by the ANSES project ECLAIR under the grant EST-2016-2 RF-23, and in part by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 724846, project 321).