Iterative solution of dense linear systems arising from the electrostatic integral equation in MEG

15 Pages·2003·0.18 MB·English

INSTITUTEOFPHYSICSPUBLISHING PHYSICSINMEDICINEANDBIOLOGY
Phys.Med.Biol.47(2002)961–975 PII:S0031-9155(02)29872-3
Iterative solution of dense linear systems arising from
the electrostatic integral equation in MEG
JussiRahola1,3andSatuTissari2,3
1SimulintuOy,Niittykummuntie6A8,02200Espoo,Finland
2CSC—ScientificComputingLtd,POBox405,02101Espoo,Finland
E-mail:satu.tissari@csc.fi
Received19October2001
Published1March2002
Onlineatstacks.iop.org/PMB/47/961
Abstract
We study the iterative solution of dense linear systems that arise from
boundary element discretizations of the electrostatic integral equation in
magnetoencephalography (MEG). We show that modern iterative methods
can be used to decrease the total computation time by avoiding the time-
consuming computation of the LU decomposition of the coefficient matrix.
More importantly, the modern iterative methods make it possible to avoid
the explicit formation of the coefficient matrix which is needed when a
large number of unknowns are used. To study the convergence of iterative
solvers we examine the eigenvalue distributions of the coefficient matrices.
For the sphere we show how the eigenvalues of the integral operator
are approximated by the eigenvalues of the coefficient matrix when the
collocation and Galerkin methods are used as discretization methods. The
collocation method approximates the eigenvalues of the integral operator
directly. The Galerkin method produces a coefficient matrix that needs to
be preconditioned in order to maintain optimal convergence speed. With
the ILU(0) preconditioner iterative methods converge fast and independent
of the number of discretization points for both the collocation and Galerkin
approaches. The preconditioner has no significant effect on the total
computationaltime.
1. Introduction
In the inverse problem of magnetoencephalography (MEG) (Ha¨ma¨la¨inenetal 1993) the
sourceofelectricalactivityinthebrainistobelocalizedfrommagneticfieldmeasurements
recorded outside the head. When the same method is applied to the heart, it
3 PartsoftheworkforthispaperwerecarriedoutintheParallelAlgorithmsProject,CERFACS,42Av.GCoriolis,
31057ToulouseCedex01,France.
0031-9155/02/060961+15$30.00 ©2002IOPPublishingLtd PrintedintheUK 961
962 JRaholaandSTissari
is called magnetocardiography (MCG). The related electrical recording techniques,
electroencephalography(EEG)andelectrocardiography(ECG),areinwideclinicaluse. The
sourcelocalizationusingEEGisamuchmorechallengingproblemthanthatofMEGbecause
ofthe spreadingof the electricpotentialon the scalp causedbythe insulatingskull. In this
paperwewillbeconcentratingontheMEG,butsomeoftheresultsareapplicabletotheother
measurementtechniques.
TosolvetheinverseproblemusingMEG,asourcedistribution,suchasacurrentdipole,
is assumed for the electric activity in the brain. In addition to this, a conductor model of
the brain is needed. In the source localization, a series of forward problemsis solved, i.e.,
the magneticfield outsidethe headcaused by a knownsource is calculated. The computed
and measured magnetic fields are then compared and the source parameters (e.g., location,
orientationandstrengthofthedipole)areadjustedusinganonlinearoptimizationmethod.
In the current paper we use a homogeneous single-compartment conductor model for
the brain. It is sufficient to solve the integral equation for the electric potential on the
surface of the brain as the skull insulates the surrounding tissues (Ha¨ma¨la¨inenandSarvas
1989). For the MEG, the magnetic field outside the head can be obtained from the
computed potential. For the accurate localization of the electric activity we need to
have a realistically shaped model for the surface of the brain given by a triangular mesh.
For the EEG source localization problems it is necessary to use a multi-compartment
integral equation where piecewise homogeneous conductivity is assumed. The multi-
compartment integral equations and their iterative solution will not be discussed in this
paper.
Inearlierstudiesofthebioelectromagneticinverseproblems,adirectsolverhasbeenused
forthe dense linear systems arising from the boundaryelementdiscretizations. Techniques
basedonLUfactorizationsofthecoefficientmatrixarepracticalbecauseaseriesofproblems
has to be solved with the same coefficientmatrix. On the other hand, the LU factorization
basedonGaussianeliminationisratherexpensivetocomputeandthusformsacomputational
bottlenecktogetherwiththeexplicitformationofthecoefficientmatrix.
Modern iterative techniques for dense linear systems have been applied in many
applicationareas. However,inthebioelectromagneticliterature,therehavebeenonlysome
references to the use of the simple Gauss–Seidel method for the inverse problem in MEG
(Ha¨ma¨la¨inenetal1993). WeapplymodernKrylov-subspaceiterativemethodstothisproblem
andpresentconvergenceandperformanceresultsfortheiterativesolvers.
Theconvergenceofiterativesolversdependsonthedistributionoftheeigenvaluesofthe
coefficientmatrices. Westudythetheoreticaleigenvaluesoftheelectrostaticintegraloperator
in the case of a spherical source region. We describe how the theoretical eigenvalues are
mappedtotheeigenvaluesofthecoefficientmatrixwhentheintegralequationisdiscretized
with the collocation and Galerkin methods. We study the eigenvalues of the coefficient
matricesalsofortherealisticbrain-shapedconductormodels.
The convergence of iterative solvers can be accelerated with the use of proper
preconditioningtechniques.Fromtheeigenvalueanalysisweproposeasimplepreconditioning
strategyandexplainwhythisstrategyshouldbeanefficientone.
Themostimportantreasonforourstudyofthemoderniterativemethodspresentedinthe
present paper is not to replace the LU decomposition with the studied iterative methods as
suchbuttocombineamoderniterativemethodwithamethodwhichmakesitpossibletosolve
thematrixequationwithoutexplicitlyformingthematrix. Inanotherpaperwehavestudied
a precorrected FFT method combined with a modern iterative method (TissariandRahola
2001). Thecombinationofthetwomethodsgivesaveryefficienttoolforthesolutionofthe
forwardproblemwithalargenumberofunknowns.
IterativesolutionofdenselinearsystemsinMEG 963
Inthetraditionalmethod,thebottleneckistheformationofthecoefficientmatrixandits
LUdecomposition,whichtakesahugeamountoftime. Thecomplexitiesofthephasesare
O(N2) forthe matrixformation, O(N3) forthe computationof the LU decomposition,and
O(N2)forthesolutionusingtheLUdecomposition(withasmallconstantofproportionality).
Incontrast,thewholecomplexityfortheiterativemethodandtheprecorrectedFFTmethodis
onlyO(NlogN). Thusforlargeenoughproblems,thewholeiterativesolutionisfasterthan
theapplicationoftheLUdecompositioninthetraditionalmethod. Inaddition,thememory
requirementsare significantly diminished. A large number of unknownsare needed in the
evaluation of the improvements in the source localization accuracy due to the realistically
shapedconductormodels.
2. Boundaryelementmethodformulation
In this section we introduce the potential integral equation and discuss the properties of
the corresponding integral operators. We introduce some notation that is used when the
eigenvaluesoftheintegraloperatorarecomputed. Wealsobrieflydiscussthediscretization
schemesandthedeflationoftheresultingmatrixequation.
2.1. Potentialintegralequationandoperator
IntheforwardproblemofMEGweneedtocomputethepotentialV(r)onthesurfaceofthe
brainduetothesource. Thebrainismodelledasahomogeneousconductorofconductivityσ
1
thatoccupiesthevolumeboundedbythesurface(cid:6). Thepotentialintegralequationisgiven
byHa¨ma¨la¨inenetal(1993):
(cid:1)
σ 1 (r(cid:3)−r)·n(r(cid:3))
V(r)=2 0V∞(r)+ V(r(cid:3)) dS(cid:3) r ∈(cid:6) (1)
σ 2π |r(cid:3)−r|3
1 (cid:6)
whereV∞(r)isthepotentialcausedbythesourceinaninfinitemediumofunitconductivity
σ ,nistheoutwardunitnormalanddS(cid:3) istheelementofarea.
0
The above integral equation is well known in potential theory (see, e.g., Sloan 1992)
andisrelatedtotheexteriorDirichletproblemfortheLaplaceequation. Theequationisan
integralequationofthesecondkindthatusestheso-calleddouble-layerrepresentationofthe
potential. Theoperatorintheequationiscalledtheelectrostaticintegraloperator.
LetusdefinetheintegraloperatorK actingonthepotentialV(r):
(cid:1) 0
1 (r(cid:3)−r)·n(r(cid:3))
K V(r)= V(r(cid:3)) dS(cid:3). (2)
0 2π |r(cid:3)−r|3
(cid:6)
ThefunctionK V(r)satisfiestheLaplaceequationforallfunctionsV(r)andallpointsrthat
0
arenotonthesurface(cid:6)(Sloan1992).
On the surface the integral operator becomes singular. When r is a point outside the
domainboundedby (cid:6) and approachesa pointr on (cid:6), the so-called jump relations(Sloan
0
1992)yield
lim K V(r)=KV(r )−V(r ) (3)
0 0 0
r→r0∈(cid:6)
wheretheelectrostaticintegraloperatorKisdefinedasanimproperintegralby
(cid:1)
1 (r(cid:3)−r)·n(r(cid:3))
KV(r)= lim V(r(cid:3)) dS(cid:3). (4)
(cid:10)→02π |r(cid:3)−r|3
(cid:6),|r(cid:3)−r|>(cid:10)
Intheaboveintegralasmallenvironmentofr isexcludedandsubsequentlyalimitistaken
wherethesizeoftheenvironmentdecreasestozero(deMunck1992).
964 JRaholaandSTissari
Theintegralinequation(1)mustbeunderstoodinthesenseofaCauchyprincipalvalue,
andtheintegralequationisthusequivalentto
σ
V(r)=2 0V∞(r)+KV(r). (5)
σ
1
2.2. Discretizationofthepotentialintegralequation
We will now discuss the discretization of the integral equation (5). We will express the
potentialasalinearcombinationofNglobalbasisfunctionsh (r):
j
(cid:2)N
V(r)= V h (r). (6)
j j
j=1
The geometry is described by a triangular mesh of the surface of the brain. We use either
piecewise constant basis functions, the unknownsbeing the values of the potential on each
triangle,orpiecewiselinearbasisfunctions,wheretheunknownsarethevaluesofthepotential
attheverticesofthetriangularmesh. Inthelattercase,eachglobalbasisfunctionisnonzero
onthetrianglesthatsharethegivenvertex.
To form the coefficient matrix from the integral equation we use either the collocation
or the Galerkin method of weighted residuals. In the collocation approach we require that
theintegralequationbesatisfiedattheso-calledcollocationpointsc . Thenumberofthese
i
pointsisthesameasthenumberofunknowns. Forthepiecewiseconstantbasisfunctionsthe
collocationpointsareatthecentresofthetriangles. Forthepiecewiselinearbasisfunctions
thecollocationpointsareattheverticesofthetriangles.
Wecannowinsertthepotentialexpansion(6)intoequation(5)andapplythecollocation
method, noting that a basis functionh (r) is one at the collocation pointc and zero at all
j j
othercollocationpoints. ThisleadstothesystemoflinearequationsfortheunknownsV
i
(cid:2)N
σ
V =2 0V∞(c )+ V Kh (c ) i =1,...,N. (7)
i i j j i
σ
1 j=1
Inmatrixnotationthisreadsas
(I −G)x =b (8)
wherexisthevectorofunknownsV ,bistheright-handsidevectorgivenby
i
σ
b =2 0V∞(c ) (9)
i i
σ
1
andtheelementsofthematrixGaregivenby
g =Kh (c ). (10)
ij j i
TheGalerkinapproachmeansthatafterthepotentialexpansion(6)issubstitutedintothe
integralequation(5),theequationismultipliedbyeachofthebasisfunctionsh(r)initsturn
i
andafterthattheresultisintegratedoverthesurface(cid:6). TheGalerkinmethodwasstudiedby
Mosheretal(1999)andindependentlybythepresentauthors(TissariandRahola1998).
Tosimplifythenotationweintroduceaninnerproductdefinedby
(cid:1)
(cid:9)f(r),g(r)(cid:10)= f(r)g(r)dS. (11)
(cid:6)
NowthelinearequationsarisingfromtheGalerkinapproachforNtriangles(i = 1,...,N)
canbewrittenas
(cid:2)N (cid:2)N
σ
V (cid:9)h (r),h (r)(cid:10)=2 0(cid:9)h (r),V∞(r)(cid:10)+ V (cid:9)h (r),Kh (r)(cid:10). (12)
j i j i j i j
σ
j=1 1 j=1
IterativesolutionofdenselinearsystemsinMEG 965
Inmatrixnotationthiscanbewrittenas
(C−F)x =b (13)
wheretheelementsofmatrixCaregivenby
c =(cid:9)h (r),h (r)(cid:10) (14)
ij i j
andthoseofFby
f =(cid:9)h (r),Kh (r)(cid:10) (15)
ij i j
whiletheright-handsidebisgivenby
σ
b =2 0(cid:9)h (r),V∞(r)(cid:10). (16)
i i
σ
1
Whenpiecewiselinearbasisfunctionsareused,thematrixCisasparsematrixwithc
ij
beingnonzeroonlyifverticesiandj shareanedgeinthetriangularmesh. ThematrixFisa
densematrix. Eachglobalbasisfunctionh (r)isnonzeroonthetrianglessharingthevertexi.
i
Inpractice,weuselocalelementwisebasisfunctionsthatarenonzeroonlywithinonetriangle,
wheretheir value is one at one of the verticesand zero at the two others. The matrices are
assembledtrianglebytriangleusinglocalbasisfunctionsandutilizinganalyticallyintegrated
elements(deMunck1992).
2.3. Deflation
Theintegralequationdeterminesthepotentialonlyuptoanadditiveconstant. Inotherwords,
ifthepotentialV(r)satisfiestheintegralequation,sodoesthepotentialV(r)+v . Interms
1
ofthecoefficientmatricesthisimpliesthatthelinearsystems(8)and(13)aresingular.
We can remove the singularity by requiring that the sum of the potentials V is zero,
i
thereby√fixingtheconstantv1. Inmatrixterms,wewilladdthenormalizedvectorofallones,
e = 1/ N(11... 1)T, multiplied by a constant γ, to all equations, so that the modified
coefficientmatrixwillbe
A+γeeT (17)
where A refers to the original coefficient matrix. This procedure is called a deflation
(Barnardetal1967).
The form of the coefficient matrix implies that e is the normalized eigenvector of the
matrix corresponding to the eigenvalue zero: Ae = 0. The deflation has the effect that
the simple eigenvalue zero is moved to the value γ and the rest of the eigenvalues remain
unchanged. Wewilldiscussthechoiceofγ inthenextsection.
3. Eigenvalues
3.1. Eigenvaluesoftheintegraloperatoranddifferentdiscretizations
Inthissectionwestudytheeigenvaluesofthepotentialintegraloperatorandthecorresponding
coefficient matrices. The distribution of the eigenvalues governs the convergence of the
iterativesolvers,andsuggestshowtodeflatethelinearsystem.
First,wewillshowhowthecollocationandGalerkinmethodswillapproximatethetrue
eigenvalues of the integral operator 1−K appearing in the integral equation (5). For this
operator,wewillwritetheeigenvalueequation
(1−K)V(r)=λV(r) (18)
966 JRaholaandSTissari
whereλisaneigenvalueandV(r)isaneigenfunction. Tosolvethisnumerically,weusethe
potentialexpansion(6). Whenweapplythecollocationmethodtothepreviousequation,we
obtainthealgebraiceigenvalueproblem
(I −G)x =λx (19)
where the matrix G and the vector x are the same as in the previous section. Thus the
eigenvalues of the coefficient matrix (1 − G) in equation (8) for the collocation method
approximatedirectlytheeigenvaluesoftheintegralequation.
FortheGalerkinmethodthediscretizedeigenvalueequationisageneralizedeigenvalue
problem
(C−F)x =λCx. (20)
ThisequationcanbepremultipliedbythematrixC−1totransformittotheordinaryeigenvalue
problem
(I −C−1F)x =λx. (21)
Thus,theeigenvaluesoftheGalerkincoefficientmatrixC−Ffromequation(13)arenotthe
sameastheeigenvaluesofthematrixI−C−1F ofequation(21)aboveandhencetheydonot
approximatetheeigenvaluesoftheintegralequationdirectly.
Forasphere,theeigenvaluesoftheintegraloperator1−Kappearinginequation(5)canbe
foundwithastraightforwardcalculationusingthesphericalharmonics(AhnerandArenstorf
1986). Theresultisgiveninthefollowingtheorem.
Theorem 3.1. When the integral operator K is applied on the surface (cid:6) of a sphere, the
eigenvaluesandeigenvectorsof1−Karegivenby
2n
(1−K)Ym(θ,φ)= Ym(θ,φ) n=0,1,... m=−n,...,n
n 2n+1 n
whereYm(θ,φ)arethesphericalharmonicsasdefinedinArfken(1985).
n
3.2. Deflation
Wecannowreturntothequestionofdeflatingthecoefficientmatrix,i.e.,thechoiceofγ in
thetransformationA → A+γeeT. Fromtheorem3.1weseethattheanalyticaleigenvalues
for the sphere convergetowards one and thus to movethe zero eigenvalueto one would be
ideal in terms of the solution of the linear systems. Therefore, we have set γ = 1 for the
collocationapproach.
In the Galerkinapproachwe would like to deflate the matrix I −C−1F with the value
γ = 1, because we know the eigenvalue distribution of this matrix for the sphere. This is
equivalenttochangingtheoriginalcoefficientmatrixC−FtoC−F +CeeT.
3.3. Comparisonofanalyticalandnumericaleigenvalues
Forasphere,wepresentacomparisonofthenumericallycomputedandanalyticeigenvalues
showingthequalityoftheeigenvalueapproximationsfromdifferentboundaryelementmethod
formulations. The numerical eigenvalues were computed using the LAPACK library. The
firsttestswererunusingadiscretizationofthespherewith720trianglesand362vertices. In
figures1and2 we showthe eigenvaluesofthecoefficientmatrixofthe collocationmethod
using the piecewise constant and linear basis functions, respectively, together with the
analytically computedeigenvaluesfor the sphere. In the plots we omit the zero eigenvalue
whichhasbeendeflatedfromthecoefficientmatrix.
IterativesolutionofdenselinearsystemsinMEG 967
x 10−8
3
2
1
λ)
g( 0
a
m
I
−1
−2
−3
0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 1.05
Real(λ)
Figure1. Numericaleigenvalues (pluses)forthecollocation methodusingpiecewise constant
basisfunctionstogetherwiththeanalyticaleigenvalues(circles)whentheunitsphereisdiscretized
with720triangles. Notethedifferentscalesfortherealandimaginaryaxes.
1
0.8
0.6
0.4
0.2
λ)
g( 0
a
m
I−0.2
−0.4
−0.6
−0.8
−1
0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 1.05
Real(λ)
Figure2.Numericaleigenvalues(pluses)forthecollocationmethodusingpiecewiselinearbasis
functionstogetherwiththeanalyticaleigenvalues(circles)whentheunitsphereisdiscretizedwith
720triangles.
Itcanbeseenfromtheplotsthatthequalityoftheeigenvaluesisapproximatelythesame
for the two collocation methods. The piecewise constant basis functions generate a larger
matrixfromthesamemesh, asthenumberofunknownsisthenumberoftriangles. Forthe
piecewiselinearbasisfunctions,thenumberofunknownsisthenumberofvertices.
968 JRaholaandSTissari
x 10−8
8
6
4
2
λ)
g( 0
a
m
I
−2
−4
−6
−8
0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 1.05
Real(λ)
Figure3. Numericaleigenvalues (pluses)fortheGalerkinmethodusingpiecewiselinearbasis
functionstogetherwiththeanalyticaleigenvalues(circles)whentheunitsphereisdiscretizedwith
720triangles. Notethedifferentscalesfortherealandimaginaryaxes.
Infigure3weshowtheanalyticaleigenvaluestogetherwiththeeigenvaluesofthematrix
I −C−1F,whereC−FisthecoefficientmatrixfortheGalerkinapproachusingpiecewise
linearbasisfunctions. NotethattheaccuracyoftheeigenvaluesismuchbetterfortheGalerkin
methodthanforthecorrespondingcollocationapproach.
To get a more quantitative measure for the accuracy of the computed eigenvalues, we
computedtheabsolutevalueoftheerrorinthesmallestnonzeroeigenvalue(2/3)forthethree
discretizationschemesandforfourdifferentdiscretizationsoftheunitsphere. Theresultsare
showninfigure4. ItcanbeseenthattheGalerkinmethodhasamuchsmallererrorthanthetwo
collocationapproachesandthatfortheGalerkinmethodtheerrordecreasesmuchfasterwhen
thediscretizationisrefined. AlsointermsofthesolutionoftheintegralequationtheGalerkin
method givesmore accurate results than the collocation methods (Mosheretal 1999). The
faster convergencerate of the Galerkin method can also be expected from theoretical error
analysesofintegralequations(Sloan1992).
Our final eigenvalue experiment consists of finding the eigenvalues for a brain-shaped
meshof886triangles. Figure5showstheeigenvaluesofthematrixI−C−1F ofthepiecewise
linearGalerkinapproachforthisrealisticallyshapedmesh. Theeigenvaluesaremorespread
out than for the sphere. Deflating the zero eigenvalue to 1 is still a good choice for the
realisticallyshapedmodel,asplentyofeigenvaluesclusteraroundthispoint.
4. Iterativemethods
4.1. Moderniterativemethods
For solving a system of linear equations, iterative methods offer an efficient alternative to
direct methods, such as the LU decomposition. Starting from an initial guess the iterative
methodsproduceiterates and the iteration can be stoppedwhen a user-definedconvergence
IterativesolutionofdenselinearsystemsinMEG 969
10−1
ue10−2
al
v
n
e
g
ei
est 10−3
all
m
s
n
Error i10−4
10−5
101 102 103 104
Number of unknowns
Figure4. Errorinthesmallestnonzeroeigenvalueasafunctionofthenumberofunknownsfor
thepiecewise constantcollocation approach (pluses), thepiecewise linearcollocation approach
(stars)andthepiecewiselinearGalerkinapproach(circles).
x 10−3
3
2
1
λ)
g( 0
a
m
I
−1
−2
−3
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
Real(λ)
Figure5. Computedeigenvalues (pluses)fortheGalerkinmethodusingpiecewiselinearbasis
functions for abrain-shaped meshof 886triangles. Note the different scales for the real and
imaginaryaxes.
criterionismet. Iterativemethodscanbeveryeffectiveifthenumberofiterationsneededfor
convergenceissmall.
Historically, the first iterative method was stationary in the sense that the iterates were
updated in the same way in each iteration. Stationary iterative methods such as the Jacobi
or Gauss–Seidel methods are not used much nowadays because of their slow convergence.
InmodernnonstationaryKrylovsubspacemethodsthe iteratesare updatedineachiteration
970 JRaholaandSTissari
based on the convergence history. Such methods access the matrix only by matrix–vector
products. Thusthewholematrixneednotbeformedexplicitly,onlyitsproductwithagiven
vectorneedstobecomputed.
For symmetric linear systems the method of choice is the conjugate gradient method
because at each iteration the error between its iterates and the true solution is the smallest
possible in a certain norm (Barretetal 1993). In addition, the new iterate of the conjugate
gradientmethodiseasytocomputefromthepreviousone.
For nonsymmetric linear systems the situation is more complicated. The generalized
minimal residual method (GMRES) produces iterates x such that the norm of the residual
k
(cid:12)Ax − b(cid:12) is the smallest possible at each iteration (SaadandSchultz 1986). However,
k
GMRES needs all the previousiterates in the computation of the new iterate and therefore
thecostofthemethodincreaseswiththeiterationnumber. Alessexpensivebutsuboptimal
method(GMRES(k))isobtainedifGMRESisrestartedeverykiterations.
Other efficient methods for nonsymmetric linear systems include the quasi-minimal
residualmethod(QMR)(FreundandNachtigal1991)andthe bi-conjugategradientmethod
stabilized(Bi-CGSTAB)(vanderVorst1992). Boththesemethodsaresimpletoapplyasthey
onlyrequirethestorageofonepreviousiterate. However,bothmethodsproduceiteratesthat
haveonlyanapproximateminimizationproperty.
4.2. Convergenceofiterativemethods
Theconvergenceofiterativesolversisgovernedbytheeigenvaluedistributionofthecoefficient
matrices. Iterativesolversconvergequicklyiftheeigenvaluesareclusteredandareallaway
from zero. In contrast, if there are both small and large eigenvalues, iterative solvers may
convergeveryslowly.
Theeigenvaluesoftheelectrostaticintegraloperatorareideallysuitedforiterativesolvers.
Ifthedeflationiscarriedoutusingthevalueγ =1,theeigenvaluesofasphericalconductor
are located on the positive real axis on the interval [2/3, 1]. The coefficient matrices for
thecollocationmethodinheritthiseigenvaluedistributiondirectly,andthustheconvergence
should be very fast. The number of iterations can be expected to be independent of the
numberofunknowns,becauseasthediscretizationisrefined,thedistributionofeigenvalues
of the coefficient matrices and thus the behaviour of iterative solvers remains practically
constant. Forthe realistically shapedbrainmodel, theeigenvaluesare morespreadout, but
thedistributionisstillverygoodforiterativesolvers.
4.3. Iterativemethodsforrealisticbrain-shapedmodel
To choose an iterative method for realistic conductor models, we studied the performance
ofunpreconditionediterativesolvers. As a stoppingcriterionand errormeasurewe use the
normwisebackwarderror(Chaitin-ChatelinandFraysse´1996)
(cid:12)b−Ax (cid:12)
η(x )= m 2 (22)
m α(cid:12)x (cid:12) +(cid:12)b(cid:12)
m 2 2
where the constant α should approximate the norm of the matrix A. In the code we use
the value α = 1, which is the two-norm of the coefficient matrix in the case of a spherical
conductormodel. The normwise backwarderror measuresthe relative error in the solution
andisindependentofthescalingoftheright-handsidevectorb.
In figure 6 we compare the convergenceof QMR, Bi-CGSTAB and restarted GMRES
foramatrixarisingfromarealisticallyshapedmeshof886trianglesandusingthepiecewise
constant collocation technique. For GMRES we used the restart parameters 5, 10 and 20.