NeuroImage 49 (2010) 641–655

Contents lists available at ScienceDirect

NeuroImagej o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / y n i m g

Robust Bayesian estimation of the location, orientation, and time course of multiple
correlated neural sources using MEG
David P. Wipf a, Julia P. Owen a, Hagai T. Attias b, Kensuke Sekihara c, Srikantan S. Nagarajan a,⁎abcDepartment of Radiology and Biomedical Imaging, University of California, San Francisco, 513 Parnassus Avenue, S362, San Francisco, CA 94143, USA
Goldenmetallic Inc., San Francisco, USA
Tokyo Metropolitan University, Tokyo, Japana r t i c l e

i n f oArticle history:
Received 21 January 2009
Revised 11 June 2009
Accepted 20 June 2009
Available online 10 July 2009a b s t r a c t
The synchronous brain activity measured via MEG (or EEG) can be interpreted as arising from a collection(possibly large) of current dipoles or sources located throughout the cortex. Estimating the number, location, andtime course of these sources remains a challenging task, one that is signiﬁcantly compounded by the effects ofsource correlations and unknown orientations and by the presence of interference from spontaneous brain
activity, sensor noise, and other artifacts. This paper derives an empirical Bayesian method for addressing each ofthese issues in a principled fashion. The resulting algorithm guarantees descent of a cost function uniquely
designed to handle unknown orientations and arbitrary correlations. Robust interference suppression is also easily
incorporated. In a restricted setting, the proposed method is shown to produce theoretically zero reconstruction
error estimating multiple dipoles even in the presence of strong correlations and unknown orientations, unlike avariety of existing Bayesian localization methods or common signal processing techniques such as beamforming
and sLORETA. Empirical results on both simulated and real data sets verify the efﬁcacy of this approach.
© 2009 Elsevier Inc. All rights reserved.

Introduction
Magnetoencephalography (MEG) and related electroencephalography(EEG) use an array of sensors to take electromagnetic ﬁeld (or voltage)
measurements from on or near the scalp surface with excellent temporal
resolution. In both MEG and EEG, the observed ﬁeld can in many cases be
explained by synchronous, compact current sources located within thebrain. Although useful for research and clinical purposes, accurately determining the spatial distribution of these unknown sources is a challenging inverse problem. The relevant estimation problem can be posed asfollows: the measured electromagnetic signal is BaRdb × dt , where db
equals the number of sensors and dt is the number of time points at whichmeasurements are made. Each unknown source Si aRdc × dt is a dcdimensional neural current dipole, at dt time points, projecting from the ith (discretized) voxel or candidate location distributed throughout thebrain. These candidate locations can be obtained by segmenting a
structural MR scan of a human subject and tesselating the brain volume
with a set of vertices. B and each Si are related by the likelihood model
B=

ds
X

Li Si + e;

ð1Þ

i=1

where ds is the number of voxels under consideration and Li aRdb × dc isthe so-called lead-ﬁeld matrix for the i-th voxel. The k-th column of Li⁎ Corresponding author. Fax: +1 415 502 4302.
E-mail address: sri@radiology.ucsf.edu (S.S. Nagarajan).
1053-8119/$ – see front matter © 2009 Elsevier Inc. All rights reserved.
doi:10.1016/j.neuroimage.2009.06.083

represents the signal vector that would be observed at the scalp given aunit current source/dipole at the i-th vertex with a ﬁxed orientation in
the k-th direction. It is common to assume dc = 2 (for MEG) or dc = 3(for EEG), which allows ﬂexible source orientations to be estimated in 2D
or 3D space. Multiple methods based on the physical properties of thebrain and Maxwell's equations are available for the computation of each
Li (Sarvas, 1987). Finally, ɛ is a noise-plus-interference term where weassume, for simplicity, that the columns are drawn independently fromN(0,Σɛ). However, temporal correlations can easily be incorporated if
desired using a simple transformation outlined in Friston et al. (2008) or
using the spatio-temporal framework introduced in Bolstad et al. (2009).
To obtain reasonable spatial resolution, the number of candidate
source locations will necessarily be much larger than the number of
sensors (ds ≫db ). The salient inverse problem then becomes the illposed estimation of regions with signiﬁcant brain activity, which are
reﬂected by voxels i such that jjSi jj N 0; we refer to these as active
dipoles or sources. Because the inverse model is severely underdeterminedhi(the mapping from source activity conﬁguration
SJ ST1 ; N ; STds T to sensor measurement B is many to one), all efforts
at source reconstruction are heavily dependent on prior assumptions,which in a Bayesian framework are embedded in the distribution p(S).
Such a prior is often considered to be ﬁxed and known, as in the caseof minimum ℓ2-norm estimation (MNE) (Baillet et al., 2001),
minimum current estimation (MCE) (Uutela et al., 1999), minimum
variance adaptive beamforming (MVAB) (Sekihara and Nagarajan,
2008), FOCUSS (Gorodnitsky et al., 1995), and sLORETA (PascualMarqui, 2002). Alternatively, empirical Bayesian approaches have

642

D.P. Wipf et al. / NeuroImage 49 (2010) 641–655

been proposed that attempt a form of model selection by using thedata, whether implicitly or explicitly, to guide the search for an
appropriate prior. Examples include a rich variety of variational
Bayesian methods and hierarchical covariance component models(Friston et al., 2008; Mattout et al., 2006; Phillips et al., 2005; Sahani
and Nagarajan, 2004; Sato et al., 2004; Wipf et al., 2007; Zumer et al.,2007). These hierarchical models can also be handled by replacing
variational inference and related procedures with Markov-Chain
Monte Carlo sampling (Nummenmaa et al., 2007). In general, theseapproaches differ in the types of covariance components that are
assumed, the cost functions and approximations used to estimate
unknown parameters, and in how posterior distributions are invoked
to affect source localization. While advantageous in many respects,
none of these methods are explicitly designed to handle complex,
correlated source conﬁgurations with unknown orientation in thepresence of background interference (e.g., spontaneous brain activity,
sensor noise, etc.).
There are two types of correlations that can potentially disrupt thesource localization process. First, there are correlations within dipole
components (meaning the individual rows of Si are correlated), which
always exists to a high degree in real data (i.e., dc N 1). For example,
dipoles with a ﬁxed unknown orientation will have a correlation
coefﬁcient of 1.0; for rotating or wobbling dipoles this value will
typically be smaller. Secondly, there are correlations between different
dipoles that are simultaneously active (meaning rows of Si arecorrelated with rows of Sj for some voxels i ≠ j). These correlations
are more application speciﬁc and may or may not exist. The larger thenumber of active sources, the greater the chance that both types ofcorrelations can disrupt the estimation process. This issue can beproblematic for two reasons. First, and most importantly, failure toaccurately account for unknown orientations or correlations canseverely disrupt the localization process, leading to a very misleading
impression of which brain areas are active. Secondly, the orientations
and correlations themselves may have clinical signiﬁcance, although
this will not be our focus herein.
In this paper, we present an alternative empirical Bayesian schemethat attempts to improve upon existing methods in terms of sourcereconstruction accuracy and/or computational robustness and efﬁciency. The Modeling assumptions section presents the basic
generative model which underlies the proposed method anddescribes the associated inference problem. This model is designed
to estimate the number and location of a small (sparse) set of ﬂexible
dipoles that adequately explain the observed sensor data. The
Algorithm derivation section derives a robust algorithm, which wecall Champagne, for estimating the sources using this model andproves that each iteration is guaranteed to reduce the associated costfunction. It also describes how interference suppression is naturally
incorporated. The Conditions for perfect reconstruction section thenprovides a theoretical analysis of conditions under which perfectreconstruction of arbitrary, correlated sources with unknown orientation is possible, demonstrating that the proposed method has
substantial advantages over existing approaches. Finally, the Empiricalevaluation section contains experimental results using our algorithm
on both simulated and real data, followed by brief conclusions in theConclusion section.
Modeling assumptions
To begin we invoke the noise model from Eq. (1), which fully
deﬁnes the assumed likelihood
2 !

ds
X

1 
pðB jSÞ~ exp − B −Li Si 
;
2
Σ−1
i=1
e

ð2Þ

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
where jjXjjW denotes the weighted matrix norm trace½X T WX . The
unknown noise covariance Σɛ will be estimated from the data using avariational Bayesian factor analysis (VBFA) model as discussed in the
Learning the interference Σɛ section; for now we will consider that it is
ﬁxed and known. Next we adopt the following source prior for S:
" d
#!
s
X
1
T −1
Si Ci Si :pðS jCÞ~ exp − trace
2
i=1

ð3Þ

This is equivalent to applying independently, at each time point, azero-mean Gaussian distribution with covariance Γi to each source Si.
We deﬁne Γ to be the dsdc ×dsdc block-diagonal matrix formed byordering each Γi along the diagonal of an otherwise
zero-valuedhimatrix.
This implies, equivalently, that pðS jCÞ~ exp − 12 trace ST C − 1 S .
If Γ were somehow known, then the conditional distributionp(S|B,Γ) ∝ p(B|S)p(S|Γ) is a fully speciﬁed Gaussian distribution withmean and covariance given byT

T −1EpðS j B;CÞ ½S = CL ðΣe + LCL Þ B
h i


T
T −1
LC;
Covpðs j B;CÞ sj = C − CL Σe + LCL
j

8j;

ð4Þ



where LJ L1 ; N ; Lds and sj denotes the j-th column of S (i.e., thesources at the j-th time point) and individual columns are uncorrelated. However, since Γ is actually not known, a suitable approximation Ĉ≈C must ﬁrst be found. One principled way to accomplish this is
to integrate out the sources S and then maximizeZpðB jCÞ =pðBjSÞpðS jCÞdS~jΣb j

−

1
2


1 T −1
exp − B Σb B ;2TΣb JΣe + LCL :ð5ÞThis is equivalent to minimizing the cost functionhi−1LðCÞJ − 2 log pðBjCÞu trace Cb Σb + log jΣb j;

ð6Þ

T
1
where CbJd−
t BB is the empirical covariance. This process is sometimes referred to as type-II maximum likelihood, evidence maximization, or empirical Bayes (Berger, 1985).
The ﬁrst term of Eq. (6) is a measure of the dissimilarity between
the empirical data covariance Cb and the model data covariance Σb; ingeneral, this factor encourages Γ to be large because it is convex andnonincreasing in Γ (in a simpliﬁed scalar case, this is akin tominimizing 1/x with respect to x, which of course naturally favors x
being large). The second term provides a regularizing or sparsifying
effect, penalizing a measure of the volume formed by the model
covariance Σb.1 Since the volume of any high dimensional space ismore effectively reduced by collapsing individual dimensions as close
to zero as possible (as opposed to incrementally reducing all
dimensions isometrically), this penalty term promotes a model
covariance that is maximally degenerate (or non-spherical), which
pushes elements of Γ to exactly zero (hyperparameter sparsity). Thisintuition is supported theoretically by the results in the Conditions for
perfect reconstruction section.
Given some type-II ML estimate Ĉ computed by minimizing Eq.
(6), we obtain the attendant empirical prior pðS j ĈÞ. To the extent that
this ‘learned’ prior is realistic, the resulting posterior pðS jB; ĈÞquantiﬁes regions of signiﬁcant current density and point estimates
for the unknown source dipoles Si can be obtained by evaluating theposterior mean computed using Eq. (4). If a given Ĉi Y0 as describedabove, then the associated Ŝ i computed using Eq. (4) also becomeszero. It is this pruning mechanism that naturally chooses the numberof active dipoles and is consistent with the hypothesis that mostregions of the brain are approximately inactive for a given task.

1
The determinant of a matrix is equal to the product of its eigenvalues, a wellknown volumetric measure.

D.P. Wipf et al. / NeuroImage 49 (2010) 641–655

Algorithm derivation
Given Σɛ and Γ, computing the Gaussian posterior on S is
straightforward as outlined above. Consequently, determining these
unknown quantities is the primary estimation task. We will ﬁrst
derive an algorithm for computing Γ assuming Σɛ is known. Later inthe Learning the interference Σɛ section, we will describe a powerfulprocedure for learning Σɛ.
Learning the hyperparameters Γ

Since direct optimization of L(Γ) may be difﬁcult, we can instead
iteratively optimize L(Γ, X, Z) via coordinate descent over Γ, X, and Z.
In each case, when two are held ﬁxed, the third can be globally
minimized in closed form. This ensures that each cycle will reduceL(Γ, X, Z), but more importantly, will reduce L(Γ) (or leave it
unchanged if a ﬁxed-point or limit cycle is reached). The associatedupdate rules from this process are as follows.
The optimal X (with Γ and Z ﬁxed) is just the standard weighted
minimum-norm solution given bynewXi
The primary objective of this section is to minimize Eq. (6) withrespect to Γ. Of course one option is to treat the problem as a general
nonlinear optimization task and perform gradient descent or someother generic procedure. Alternatively, several methods in the MEG
literature rely, either directly or indirectly, on a form of the
expectation maximization (EM) algorithm (Friston et al., 2008; Satoet al., 2004). However, these algorithms are exceedingly slow when ds
is large and they have not been extended to handle arbitrary,
unknown orientations. Consequently, here we derive an optimization
procedure that expands upon ideas from Sato et al. (2004) and Wipf
et al. (2007), handles arbitrary/unknown dipole orientations, andconverges quickly.
To begin, we note that L(Γ) only depends on the data B through the
db × db sample correlation matrix Cb. Therefore, to reduce thefcomputational burden, we replace B with a matrix BaRdb × rankðBÞ
ffT
such that B B = Cb . This removes any per-iteration dependency on
dt, which can potentially be large, without altering the actual costfunction. It also implies that, for purposes of computing Γ, the numberof columns of S is reduced to match rank(B). We now re-express thecost function L(Γ) in an alternative form leading
to convenient  update
rules and, by construction, a proof that L Cðk + 1Þ VL CðkÞ at each
iteration.
The procedure we will use involves constructing auxiliary functions using sets of hyperplanes; an introduction to the basic ideas canbe found in Appendix A. First, the log-determinant term of L(Γ) is aconcave function of Γ and so it can be expressed as a minimum overupper-bounding hyperplanes via
log jΣb j = min
Zz0

" d
s
X

#


T

trace Zi Ci − h ðZ Þ ;

ð7Þ

i=1

X

"
#
2
ds
ds  2
X
X
f

 
 B −
X 
L
X
+
i
i
i

 −1
  −1 ;
i=1

Σe

i=1

ð8Þ

Ci

h
iT
is a matrix of auxiliary variables as before.
where XJ X1T ; N ; XdTs
Note that in this case, the implicit concave conjugate function exists inclosed form.
Dropping the minimizations and combining terms from Eqs. (7)and (8) lead to the modiﬁed cost function

2
ds
f X

LðC; X; Z Þ=  B −Li Xi 
i=1

Σ−1
e

+

" 
#dsX 2X  + trace Z T C− h ð Z Þ;i i
 i  −1

i=1

Ci

ð9Þ
where by construction L(Γ) = minXminZL(Γ, X, Z). It is straightforward
to show that if fĈ; X̂; Ẑg is a local (global) minimum to L(Γ, X, Z),
then C is a local (global) minimum to L(Γ).

T −1 ~YCi Li Σb B

ð10Þfor each i. The minimizing Z equals the slope at the current Γ of log|Σb|, which follows from simple geometric considerations (Boyd andVandenberghe, 2004). As such, we have
new

Zi

T

−1

YjCi log jΣb j = Li Σb Li :

ð11Þ

With Z and X ﬁxed, computing the minimizing Γ is a bit more difﬁcult
because of the constraint Γi ∈ H+ for all i, where H+ is the set of
positive-semideﬁnite, symmetric dc × dc covariance matrices. To
obtain each Γi, we must solve" # 2TX 
 i  −1 + trace Zi Ci :new
Ci Y arg minþ
Ci aH

ð12Þ

Ci

An unconstrained solution will satisfy
jCi LðCi ; Xi ; Zi Þ = 0;

ð13Þwhich, after computing the necessary derivatives and re-arranging
terms gives the equivalent condition
T

Xi Xi = Ci Zi Ci :

ð14ÞThere are multiple (unconstrained) solutions to this equation; we will
choose the unique one that satisﬁes the constraint Γi ∈ H+. This can befound usingT−1 = 2Zi −1 = 21=2
T 1=2 1=2
1=2
T 1 = 2 1 = 2 −1 = 2= Zi
Zi Xi Xi Zi
Zi Xi Xi Zi
Zi1 = 2−1 = 21=2
T 1=2−1 = 2= Zi
Zi Xi Xi Zi
Zi
Zi−1 = 2Xi Xi = Z ihiTwhere ZJ Z1T ; N ; ZdTs
is a matrix of auxiliary variables thatdifferentiates each hyperplane and h⁎(Z) is the concave conjugate oflog |Σb|. While h⁎(Z) is unavailable in closed form, for our purposes
below, we will never actually have to compute this function. Next, thedata ﬁt term is a concave function of Γ− 1 and so it can also be
expressed using similar methodology ashi−1
trace Cb Σb = min

643



1= 2

Zi

−1 = 2

 Zi

T 1= 2Xi Xi Zi

1=2
T 1 = 2 1 = 2 −1 = 2
Zi Xi Xi Zi
Zi
:

ð15Þ

This indicates (via simple pattern matching) the solution (or update
equation)newCi−1 = 2

YZi



1=2

Zi



T 1 = 2 1 = 2 −1 = 2
Zi;

Xi Xi Zi

ð16Þ

which satisﬁes the constraint. And since we are minimizing a convex
function of Γi (over the constraint set), we know that this is indeed aminimizing solution.
In summary then, to estimate Γ, we need simply iterate Eqs. (10),(11), and (16), and with each pass we are guaranteed to reduce (orleave unchanged) L(Γ); we refer to the resultant algorithm as Champagne. The per-iteration cost is linear in the number of voxels ds so thecomputational cost is relatively modest (it is quadratic in db, and cubicin dc, but these quantities are relatively small). The convergence rate isorders of magnitude faster than EM-based algorithms such as those in
Friston et al. (2008) and Sato et al. (2004) (Empirical evaluation
section contains a representative example).

644

D.P. Wipf et al. / NeuroImage 49 (2010) 641–655

Simpliﬁcations using constraints on Γ

Conditions for perfect reconstruction

Here we consider two constraints on the parameterization of Γ thatlead to much less complex updates and provide connections withexisting algorithms. First, we can require that off-diagonal terms ofeach Γi are equal to zero, i.e., the prior covariance of each source Si is
diagonal. We then restrict ourselves to learning the diagonal elementsof Γi via simpliﬁed versions of those presented above. A second
possibility is to further constrain each Γi to satisfy Γi = γiI, where γi is ascalar non-negative hyperparameter. In both cases, the resulting costfunctions and algorithms also fall out of the framework we discuss in
Wipf (2006) and Wipf et al. (2007). These variants are also similar tothe covariance component estimation model used in Mattout et al.
(2006) and Friston et al. (2008), albeit with a much larger number ofcomponents here.
As we will see in the Conditions for perfect reconstruction and
Empirical evaluation sections, the reduced parameterization associatedwith these models can potentially degrade performance in somesituations. Nonetheless, these approaches are very useful for comparison purposes. To distinguish all three cases, we use the designations
CHAMPS for the scalar version, CHAMPD, for the arbitrary diagonalcase, and CHAMPM when the matrices Γi are unrestricted. For the
special case where dc = 2, this gives the following parameterizations:

Whenever a new inverse algorithm is proposed, it is often insightful
to know what source conﬁgurations it can recover exactly and underwhat conditions. While in general this will require unrealistic
assumptions such as zero noise and a perfect forward model, theunderlying idea is that if only implausible source conﬁgurations can berecovered even given these strong simpliﬁcations, then perhaps thealgorithm may have serious difﬁculties.
For example, consider the classical MNE method (Baillet et al.,2001). Here exact recovery of the true sources requires that columns of
the source activity matrix S lie in the range space of the transposed
lead-ﬁeld, i.e., S = LTA, for some coefﬁcient matrix A. But this is a very
contrived requirement, because the lead-ﬁeld transpose is highly
overdetermined (many more rows than columns) meaning that thetrue S must be constrained to a small subspace of the total possible
source space, i.e., a subspace of Rds dc × dt . But this subspace has nothing
to do with any source activity or neurophysiological plausibility. The
lead-ﬁeld, which determines this subspace, only speciﬁes how a given
source conﬁguration maps to the sensors, it is unrelated to what theactual activity is. Therefore, MNE is effectively constraining thepossible source space in an ad hoc manner even before confounding
noise or forward modeling errors are introduced.
A related concept is the notion of localization bias, which isbasically a way of quantifying the ability of an algorithm to accurately(or sometimes perfectly) estimate a single dipolar source. Recently it
has been shown, both empirically and theoretically (Pascual-Marqui,
2002; Sekihara et al., 2005), that the MVAB and sLORETA algorithms
have zero location bias given noiseless data and an ideal forward
model, meaning their respective source estimates will peak exactly at
the true source. Note that this is a slightly different notion than perfectreconstruction, since both methods will still (incorrectly) produce
non-zero source estimates elsewhere. These ideas have been extended
to include certain empirical Bayesian methods (Sato et al., 2004; Wipf
et al., 2007). However, these results assume a single dipole with ﬁxed,
known orientation (i.e., dc = 1), and therefore do not formally handle
more complex issues that arise when multiple dipoles with unknown
orientations are present. The methods from Sahani and Nagarajan(2004) and Zumer et al. (2007) also purport to address these issues,
but no formal analyses are presented.
Despite its utilization of a complex, non-convex cost function L(Γ),
we now demonstrate relatively general conditions whereby Champagne can exhibit perfect reconstruction. We will assume that the fulllead-ﬁeld L represents a sufﬁciently high sampling of the source spacesuch that any active dipole component aligns with some lead-ﬁeld
columns (the ideal forward model assumption). Perfect reconstruction
can also be shown in the continuous case, but the discrete scenario is
more straightforward and of course more relevant to any practical task.
Some preliminary deﬁnitions are required to proceed. We deﬁne theempirical intra-dipole correlation matrix at the i-th voxel as Cii J d1t Si STi ;non-zero off-diagonal elements imply that correlations are present.
Except in highly contrived situations, this type of correlation will always
exist. The empirical inter-dipole correlation matrix between voxels i andj is Cij J d1t Si STj ; any non-zero element implies the existence of acorrelation. In practice, this form of correlation may or may not bepresent. With regard to the lead-ﬁeld L, spark is deﬁned as the smallest
number of linearly dependent columns (Donoho and Elad, 2003). By
deﬁnition then, 2 ≤ spark (L) ≤db + 1. Finally, da denotes the number ofactive sources, i.e., the number of voxels whereby jjSi jj N 0.

CHAMPS
γi 0
Ci =
0 γi

CHAMPDγ i1 0Ci =0 γ i2

CHAMPM
γi11 γ i12
Ci =:γi12 γ i22

ð17ÞLearning the interference ΣɛThe learning procedure described in the previous section boils
down to ﬁtting a structured maximum likelihood covariance estimate
Σb = Σɛ + LΓLT to the data covariance Cb. The idea here is that LΓLT will
reﬂect the brain signals of interest while Σɛ will capture all interferingfactors, e.g., spontaneous brain activity, sensor noise, muscle artifacts,
etc. Since Σɛ is unknown, it must somehow be estimated or otherwise
accounted for. Given access to pre-stimulus data (i.e., data assumed tohave no signal/sources of interest), stimulus evoked partitioned factor
analysis provides a powerful means of decomposing a data covariancematrix Cb into signal and interference components. While details canbe found in Nagarajan et al. (2007), the procedure computes theapproximation
TT

Cb ≈Λ + EE + FF ;

ð18Þ

where EaRdb × de represents a matrix of learned interference factors, Λ
is a diagonal noise matrix, and FaRdb × df represents signal factors.
Both the number of interference factors de and the number of signal
factors df are learned from the data via a variational Bayesian factor
analysis procedure. Using a generalized form of the expectation
maximization algorithm, the method attempts to ﬁnd a small numberof factors that adequately explains the observed sensor data
covariance during both the pre- and post-stimulus periods. The prestimulus is modeled with a covariance restricted to the terms Λ + EET,
while the post-stimulus covariance, which contains the signal
information FFT we wish to localize, is expressed additively as in Eq.
(18).
There are two ways to utilize the decomposition (Eq. (18)). First,
we can simply set Σɛ→Λ + EET and proceed as in the Learning thehyperparameters section. Alternatively, we can set Σɛ→0 (or set it to asmall diagonal matrix for robustness/stability) and then substitute FFT
for Cb, i.e., run the same algorithm on a de-noised signal covariance. In
practice, both methods have been successful; however, a full technicaldiscussion of the relative merits is beyond the scope of this paper.

Theorem 1. In the limit as Σɛ→0 (high SNR) and assuming dad-c b spark (L)−1, the cost function L(Γ) maintains the following two

properties:
1) For arbitrary Cii and Cij, the unique global minimum Γ ⁎ produces asource estimate S⁎ = Ep(S∣B,Γ⁎)[S] computed using Eq. (4) that

D.P. Wipf et al. / NeuroImage 49 (2010) 641–655

equals the generating source matrix S, i.e., it produces a perfectsource reconstruction.
2) If Cij = 0 for all active dipoles (although Cii is still arbitrary), thenthere are no local minima, i.e., the cost function is unimodal.
See the Appendix B for the proof. In words, this theorem says that
intra-dipole correlations do not disrupt the estimation process bycreating local minima, and that the global minimum always achieves aperfect source reconstruction. In contrast, inter-dipole correlations
can potentially create local minima, but they do not affect the global
minimum. Empirically, we will demonstrate that the algorithm
derived in the Algorithm derivation section is effective at avoiding
these local minima (see Empirical evaluation section). With added
assumptions these results can be extended somewhat to handle theinclusion of noise.
The cost functions from Sato et al. (2004) and Wipf et al. (2007)
bear the closest resemblance to L(Γ); however, neither possesses thesecond attribute from Theorem 1. This is a signiﬁcant failing because,
as mentioned previously, intra-dipole correlations are always present
in each active dipole. Consequently, reconstruction errors can occur
because of convergence to a local minimum. The iterative Bayesian
scheme from Zumer et al. (2007), while very different in structure,also directly attempts to estimate ﬂexible orientations and handle, tosome extent, source correlations. While details are omitted for brevity,
we can prove that the full model upon which this algorithm is based
fails to satisfy the ﬁrst property of the theorem, so the corresponding
global minimum can fail to reconstruct the sources perfectly. In
contrast, beamformers and sLORETA are basically linear methods with
no issue of global or local minima. However, the popular sLORETA and
MVAB solutions will in general display peaks that may be misaligned
from the true sources for multi-component dipoles (dc N 1) or when
multiple dipoles (da N 1) are present, regardless of correlations (theywill of course never produce perfect reconstructions of dipolar sources
because of spatial blur).
So in summary, our analysis provides one level of theoreticalcomparison. In some sense it is relevant to questions of the sort, do we
expect the underlying true current sources to lie in the subspace
formed by S = LTA, where LT is essentially an arbitrary overdeterminedmatrix in this context, or will they more likely be of the form S equals
the sum of da b sparkðdLcÞ − 1 arbitrary dipoles with unconstrained orientations and locations? This latter conﬁguration is a actually equivalent tothe union of a large number
 of dadc-dimensional subspaces (in fact the
number is combinatorial, dds , which represents all possible uniqueaconﬁgurations of da dipoles), a far more complex, and in some sense
richer, set of possible sources. Champagne is designed to recover sources
resembling the latter, an approximation to ‘real’ source conﬁgurations
that many neurophysiologists would say is both plausible and usefulclinically, especially in an event-related paradigm.
Empirical evaluation
In this section we test the performance of our algorithm on both
simulated and real data sets. We focus here on localization accuracy
assuming strong source correlations and unknown orientations. Notethat the primary purpose of the proposed algorithm is not really toestimate the actual orientations or correlations per se. It is toaccurately estimate the location and power of sources confoundedby the effects of unknown orientations and correlations. Consequently, accurate localization estimates implicitly indicate that these
confounds have been adequately handled, hence orientation (orcorrelation) estimates themselves are not stressed.
Simulated data
We ﬁrst conducted tests using simulated data with realistic source
conﬁgurations. The brain volume was segmented into 5 mm voxels

645

and a two orientation (dc = 2) forward lead-ﬁeld was calculated using
a single spherical-shell model (Sarvas, 1987). The data time courses
were partitioned into pre- and post-stimulus periods. In the prestimulus period (263 samples) there is only noise and interferingbrain activity, while in the post-stimulus period (437 samples) there is
the same (statistically) noise and interference factors plus source
activity of interest. The pre-stimulus activity consisted of the restingstate sensor recordings collected from a human subject presumed tohave spontaneous activity (i.e., non-stimulus evoked sources) andsensor noise; this activity was on-going and continued into the poststimulus period, where the simulated source signals were added.
Source time courses were seeded at locations in the brain withdamped-sinusoidal signals and this voxel activity was projected to thesensors through the lead-ﬁeld. The locations for the sources werechosen so that there was some minimum distance between sources
and a minimum distance from the center of the head. We could adjust
both the signal-to-noise-plus-interference ratio (SNIR), the correlations between the different voxel time courses (inter-dipole), and thecorrelations between the two orientations of the dipoles (intradipole) to examine the algorithm performance on unknown correlated sources and dipole orientations. For our purposes, SNIR isdeﬁned as
SNIRJ20 log

‖LS‖F
:
‖E‖F

ð19Þ

To obtain aggregate data on the performance of our method on manydifferent dipole conﬁgurations and noise levels, we ran 100 simulations of three randomly (located) seeded sources at SNIR levels of − 5,0, 5, 10 dB. The sources in these simulations always had an inter-dipole
correlation coefﬁcient and an intra-dipole correlation coefﬁcient of
0.95. We chose to test our algorithm against ﬁve representative source
localization algorithms from the literature: minimum varianceadaptive beamforming (MVAB) (Sekihara and Nagarajan, 2008), two
non-adaptive spatial ﬁltering methods, sLORETA (Pascual-Marqui,2002) and dSPM (Dale et al., 2000), and two variants of minimum
current estimation (MCE) specially tailored to handle multiple time
points and unconstrained dipoles (Wipf, 2006; Wipf and Nagarajan,2009). These two methods extend standard MCE by applying a ℓ2
norm penalty across time (an ℓ1 norm over space and an ℓ2 norm
over time, sometimes called an ℓ1,2 norm in signal processing). In one
case both dipole components are also included within the ℓ2 penalty(MCE1), in the other case they are treated individually (MCE2). Similar
to Champagne, both versions of MCE favor sparse/compact sourcereconstructions. Related algorithms are discussed in Bolstad et al.
(2009) and Huang et al. (2006).
We ran the simulations using a total of eight algorithms, the ﬁveabove plus the three variants of our Champagne method: CHAMPS,CHAMPD, and CHAMPM (see Simpliﬁcations using constraints on Γ
section for a description of these variants). In order to evaluate
performance, we used two features: source localization accuracy andtime-course estimation accuracy. To assess localization accuracy, weused the A′ metric (Snodgrass and Corwin, 1988) and to assess theestimation of the time courses, we used the correlation coefﬁcient
between the true and estimated time courses.
When assessing the localization accuracy, it is important to take
into account both the number of hits (sources that were correctly
localized) and the presence of false positives or spurious localizations.
A principled way to take these two features into account is the ROC(receiver-operator characteristic) method modiﬁed for brain imaging
results (Darvas et al., 2004), which is a measure of hit rate (HR) versus
false positive rate (FR). Speciﬁcally, we used the A′ metric which is a
way to approximate the area under the ROC for one HR/FR pair. We
determined the HR and FR at each SNIR level and each algorithm in thefollowing way. For each simulation, we calculated all the local peaks in
the image. A local peak is deﬁned as a voxel that is greater than its 30

646

D.P. Wipf et al. / NeuroImage 49 (2010) 641–655

Fig. 1. Performance evaluation. Left: Aggregate localization results (A′ metric) for MVAB, sLORETA, dSPM, two variants of MCE (MCE1 and MCE2), and the three variants of Champagne(CHAMPS, CHAMPD, and CHAMPM) for recovering three correlated sources with unknown orientations. Right: Estimated time-course correlation coefﬁcient results for the eightalgorithms. Error bars denote the standard error over 100 simulations.

three-dimensional nearest neighbors and is at least 10% of the
maximum activation of the image. (This thresholding at 10% is
designed to ﬁlter out any spurious peaks or ripples in the image thatare much weaker than the maximum peak.)After all the local peak locations were obtained, we tested whether
each true source location has at least one local peak within acentimeter of it. (We also tested the performance at a radius of half acentimeter and the trends were similar for all algorithms.) If there
were multiple local peaks within one centimeter of a particular true
source, we counted it as one hit. Each image has the possibility of
having 0, 1, 2, or 3 hits and we divided the number of hits by three toget the HR. Determining the false positive rate is more tenuous. There

Fig. 2. Localization of three highly correlated dipoles at SNIR = 5 dB, the green circles denote the true locations of the sources and the red-to-white colored surface plot shows the
maximum-intensity projection of the source power estimates produced by each algorithm. The color scale gradation goes from dark being the maximum to light being the minimum.

D.P. Wipf et al. / NeuroImage 49 (2010) 641–655

647

Fig. 3. Localization of 3 highly correlated dipoles at SNIR = − 5 dB, the green circles denote the true locations of the sources and the red-to-white colored surface plot shows the
maximum-intensity projection of the source power estimates produced by each algorithm. The color scale gradation goes from dark being the maximum to light being the minimum.

is not a clear maximum number of possible false positives, as there iswith hits. We decided that given the spatial smoothness of the images,
we could place a ceiling on the number of local peaks in an image at 5%of the total number of voxels or 100 voxels in the case of oursimulations. Thus we divided the number of false positives, those local
peaks not within one centimeter of a seeded source location, by 100 to
get the FR. Since each algorithm was treated the same for thedetermination of the FR, we do not think that this estimated ceiling on
the false positives caused any bias in our performance evaluation.

For evaluating the time-course reconstruction we used thecorrelation coefﬁcient between the true time course seeded for thesimulations and the estimated time series from each algorithm. Boththe A′ metric and correlation coefﬁcient range from 0 to 1, with 1
implying perfect localization and time-course estimation. For both thelocalization accuracy and the time-course estimation assessments 100
simulations were averaged. Fig. 1 displays comparative results for theeight algorithms at different SNIR levels with standard errors. Fig. 1 leftdemonstrates the A′ metric and Fig. 1 right shows the time-course

Fig. 4. Sensor data from all 275 MEG channels associated with Fig. 2 example. The redline denotes the pre- and post-stimulus periods.

Fig. 5. Sensor data from all 275 MEG channels associated with Fig. 3 example. The redline denotes the pre- and post-stimulus periods.

648

D.P. Wipf et al. / NeuroImage 49 (2010) 641–655

correlation coefﬁcients. All three variants of Champagne quite
signiﬁcantly outperform the others. For all further analyses, only
CHAMPM was used because of its superior performance.

Figs. 2 and 3 show sample reconstructions of three correlated
dipoles at 5 dB and − 5 dB SNIR respectively. The associated sensor
data is depicted in Figs. 4 and 5. Inter- and intra-dipole correlation

Fig. 6. Localization of 10 dipoles. Left column: Uncorrelated sources; right column: correlated sources. The green circles denote the true locations of the sources and the red-to-whitecolored surface plot shows the maximum-intensity projection of the source power estimates produced by each algorithm. The color scale gradation goes from dark being the
maximum to light being the minimum.

D.P. Wipf et al. / NeuroImage 49 (2010) 641–655

649

Fig. 7. Deep source example: 3 dipoles were seeded close to the center of the head to test the performance of the algorithm on deep sources. Champagne does the best at recovering
these three sources.

coefﬁcients were 0.90 in each case. All of the surface plots are
maximum-intensity projections of estimated power maps onto thecoronal plane. These power maps were obtained by calculating the
average power at every voxel from the reconstructed voxel timecourses. The green circles indicate the projection of the true locationsof the sources and the surface plot shows the maximum-intensity
projection of the power-map. For these images the black/dark red
regions are the maximum and the minimum is pale yellow/white. The
ﬁgures show reconstructions from CHAMPM against MVAB, MCE1, andsLORETA/dSPM. In practice, the localization maps for sLORETA anddSPM are nearly identical (e.g., see Fig. 1), so we chose to only show
the sLORETA images to represent both methods. CHAMPM performs

signiﬁcantly better at localizing the activity at both high and low SNIR
levels with highly correlated dipoles.
Fig. 6 shows a sample reconstruction of a much more complex
source conﬁguration involving 10 dipolar sources. We compared
CHAMPM versus MVAB, MCE1, and sLORETA/dSPM assuming (i) the
dipoles had zero inter- and intra-dipole correlations, and (ii) they had
a correlation coefﬁcient of 0.5 for both inter- and intra-dipole
interactions. The SNIR for these simulations was 10 dB. CHAMPM
performs signiﬁcantly better than the other algorithms on both thecorrelated and uncorrelated cases. Note that the reason superﬁcial
sources are easier to ﬁnd via Champagne (and with many othermethods) is simply because the effective SNIR of deep sources is

Fig. 8. Localization of diffuse sources. 10 clusters of 10 dipoles each were seeded in the brain and reconstructed using various algorithms. True and estimated sources were thenprojected to the rendered surface of the brain. Champagne comes the closest to matching the true sources.

650

D.P. Wipf et al. / NeuroImage 49 (2010) 641–655

Fig. 9. Convergence rate of proposed update rules relative to a conventional EM
implementation based on Friston et al. (2008); Sato et al. (2004) and Wipf (2006).

signiﬁcantly lower. This occurs because deep sources, when combinedwith superﬁcial ones, contribute very little to the numerator of Eq.
(19) because the norm of the associated lead-ﬁeld columns isrelatively small. It is not because there is an intrinsic bias favoring
superﬁcial dipoles; in fact, the global minimum of L(Γ) is invariant to
the lead-ﬁeld column norms in the sense that the activation/sparsity
pattern is completely unaffected.
As evidence that deeper sources need not necessarily be more
difﬁcult to ﬁnd for a ﬁxed SNIR, consider a scenario with three sourcesthat are all reasonably deep. In this situation, the relative SNIR of eachsource will be comparable, meaning the overall source power issomewhat evenly distributed. If we are assuming a ﬁxed SNIR forthree superﬁcial sources versus three deep ones, then the power
associated with the deep ones will necessarily be larger to compensate

for the small lead-ﬁeld values. As a result, these deep sources will be ofsimilar difﬁculty to ﬁnd. So it is really in the situation where the SNIR
is ﬁxed but both deep and superﬁcial sources are present that trouble
arises. The shallow sources will dominate the SNIR calculation andwill be readily located, while the deep ones contribute almost nothing
to the data covariance and will likely be overlooked.
Fig. 7 shows recovery results using three deep sources with 0.5
correlation (inter- and intra-dipole) and 2 dB SNIR, with interference
from real brain noise as in previous experiments. Clearly Champagne
has no difﬁculty in this scenario while the other algorithms struggle.
Of course in a practical situation, the SNIR for deep sources may bewell below 2 dB, rendering them more difﬁcult to ﬁnd.
We also tested Champagne on more distributed source conﬁgurations. Previously, we created sources constrained to one voxel; now
we consider sources formed from activity in 10 adjacent voxels toform clusters. We seeded 10 of these source clusters throughout thebrain volume of interest. The lead-ﬁeld used for this experiment is thesame as the previous experiments with a spatial resolution of 5 mm.
The SNIR was 10 dB and the intra-dipole correlation was 0.5 in allcases. The inter-dipole correlation was 0.8 for dipoles within the samecluster and 0.5 between dipoles in different clusters. We visualized theclusters by projecting the source power estimates to the surface of the
3-D rendered MRI.
From Fig. 8 we observe that Champagne was able to resolve muchof each cluster and 80 of the 100 voxels were found to be signiﬁcantly
active above any small spurious peaks. Comparing with the groundtruth plot, Champagne was able to reconstruct the clusters quite
accurately, while MCE, sLORETA, and MVAB are not able to do so.
Consequently, Champagne is found to be both applicable to focal
sources as well as small diffuse clusters.
Finally, Fig. 9 gives an example of the improvement in convergence
rate afforded by our method relative to an EM implementation
analogous to (Friston et al., 2008; Sato et al., 2004); the per-iteration
cost is the same for both methods. A simulated 2D dipolar source was
generated and projected to the sensors using the experimental

Fig. 10. Localization of SEF data. All four algorithms localize to somatosensory cortical areas, but Champagne and MCE are the most focal.

D.P. Wipf et al. / NeuroImage 49 (2010) 641–655

651

Fig. 11. Localization of AEF data (subject 1). Only Champagne is able to localize bilateral activity in primary auditory cortex.

paradigm described in Zumer et al. (2007) with ds = 1725 voxels. The
signal was corrupted by 10 dB additive Gaussian sensor noise. Fig. 9
displays the reduction in the Champagne cost function L(Γ) as afunction of the iteration number. Consistent with previous observations, the EM updates are considerably slower in reducing the cost.
While the detailed rationale for this performance discrepancy isFig. 12. Localization of AEF data (subject 2). Only Champagne is able to localize bilateral activity in primary auditory cortex.

652
D.P. Wipf et al. / NeuroImage 49 (2010) 641–655
Fig. 13. Top row: Sensor data from all 275 MEG channels associated with Figs. 10, 11, and 12. Bottom row: Recovered time courses at the maximum-intensity voxels as estimated via CHAMPM. Left: Left somatosensory cortex showing peak at
50 ms. Middle: Left auditory cortex from subject 1 (right auditory cortex, not shown, is similar) showing M50 and M100. Right: Right auditory cortex from subject 2 (left auditory cortex, not shown, is similar) showing M50 and M100.

D.P. Wipf et al. / NeuroImage 49 (2010) 641–655

beyond the scope of this paper, ultimately it is a consequence of thedifferent underlying bounds being used to form auxiliary functions.
EM leads to slower convergence because it is effectively using a muchlooser bound around zero than the bound described in the Learning
the hyperparameters section and therefore fails to fully penalize
redundant or superﬂuous components. This prevents the associatedhyperparameters from going to zero very quickly, drastically slowing
the convergence rate. More details on this topic can be found in Wipf
and Nagarajan (2009).
Real data
Three stimulus evoked data sets were collected from normal,
healthy research subjects on a 275-channel CTF System MEG device.
The ﬁrst data set was a sensory evoked ﬁeld (SEF) paradigm, where
the subject's right index ﬁnger was tapped for a total of 240 trials. A
peak is typically seen 50 ms after stimulation in the contralateral (inthis case, the left) somatosensory cortical area for the hand, i.e., dorsal
region of the postcentral gyrus. Champagne and MCE1 were able tolocalize this activation to the correct area of somatosensory cortex as
seen in Fig. 10 and the estimated time course shows the typical 50 ms
peak Fig. 13 (left). The other algorithms were also able to localize thisactivity, but the estimated activations are much more diffuse.
Champagne and MCE's sparsity characteristics result in focal activations that do not have spatial blur like MVAB and sLORETA/dSPM.
Note that for these examples, we do not perform a maximumintensity projection. Instead, we simply ﬁnd the voxels withmaximum power (obtained from the estimated time courses).
Champagne and MCE's activations were very focal, but for the othertwo algorithms we thresholded the image of estimated source power
to obtain a focal activation around the maxima.
The other two data sets analyzed were from an auditory evokedﬁeld (AEF) paradigm, where two subjects were presented tones
binaurally for a total of 120 trials. There are two typical peaks seen
after the presentation of an auditory stimulus, one at 50ms and one at
100 ms, called the M50 and M100 respectively. The auditory
processing of tones is bilateral at early auditory cortical areas andthe activations are correlated. Champagne was able to localize activity
in both primary auditory cortices for both subjects, Figs. 11 and 12, andthe time courses for these two activations reveal the M50 and M100
Fig. 13 (middle) and (left). In general, these ﬁgures demonstrate that
ChampM outperforms both MVAB, MCE, and sLORETA/dSPM. The
analysis of simple auditory paradigms is problematic for MVAB
because it cannot handle the bilateral correlated sources well.
sLORETA/dSPM does show some degree of bilateral activation but
there is a large amount of spatial blur, while Champagne is focal in its
localization. Note that the ﬁrst iteration of Champagne is very related
to sLORETA/dSPM (see Wipf et al., 2007) so this result is perhaps not
surprising. Finally, MCE like Champagne produces focal estimates, butwith additional spurious peaks and bilateral activations outside of
auditory cortical areas.
Conclusion
This paper derives a novel empirical Bayesian algorithm for MEG
source reconstruction that readily handles multiple correlated sources
with unknown orientations, a situation that commonly arises evenwith simple imaging tasks. Based on a principled cost function andfast, convergent update rules, this procedure displays signiﬁcanttheoretical and empirical advantages over many existing methods. We
have restricted most of our exposition and analyses to MEG; however,preliminary work with EEG is also promising. For example, on a realworld passive visual task where subjects viewed ﬂashing foreground/background textured images, our method correctly localizes activity to
the lateral occipital cortex while two state-of-the-art beamformers
fail. This remains an active area of research.

653

There are a number of directions for future research. First, while
not the focus of this paper, our model can readily be augmented to
produce model evidence approximations, which have been proposed
for Bayesian model comparison/selection purposes (Friston et al.,2008, 2006; Wipf and Nagarajan, 2009). While here such approximations emerge from the concavity bounds associated with thederivation of Champagne in the Learning the hyperparameters
section, in Friston et al. (2008, 2006), related bounds follow fromfactorial (or mean-ﬁeld) assumptions built into a variational freeenergy framework. The connection between these two different typesof bounds, as well as their relative performance in model selection
tasks, has not been explored to our knowledge.
Secondly, a full investigation of the effects of the temporalwindowing of sensor time courses, and algorithmic extensions forlearning optimal windows, is very important to empirical Bayesian
methods such as Champagne. There is an intrinsic trade-off here.
Large windows are most effective for localizing stationary sources thatremain active across the whole window length; however, if a source
has limited temporal extent or is moving, extending the window size
can drastically reduce performance. In contrast, small windows are
optimal for non-stationary, ephemeral sources, but the effective SNIR
will generally be worse. In preliminary studies with data from facial
processing experiments, adjusting the window size has been crucial
for localizing regions such as the fusiform face area. This issue isimplicitly considered in Bolstad et al. (2009), where temporal basis
functions are used in extended MCE framework (a notion that could
potentially be applied to Champagne as well). Spatial basis functions
can also be incorporated to allow more ﬂexible estimation of
distributed sources, a useful notion that has been applied in a wide
variety of settings (Bolstad et al., 2009; Limpiti et al., 2006; Phillips etal., 2002; Ramírez and Makeig, 2006; Wipf and Nagarajan, 2009).
Acknowledgments
This work was funded by NIH grants R01 DC006435 and DC004855
to SSN.
Appendix A. Introduction to conjugate duality
At the heart of the algorithm derived in the Algorithm derivation
section is the ability to represent a concave function in its dual form.
For example, given a concave function f ð yÞ : RYR, the dual form isgiven by 
f ð yÞ = inf υy − f ðυÞ ;υð20Þ

where f ⁎(υ) denotes what is called the conjugate function (Boyd andVandenberghe, 2004). Geometrically, this can be interpreted as representing f(y) as the lower envelope or inﬁmum of a set of lines
parameterized by υ. The selection of f ⁎(υ) as the intercept term
ensures that each line is tangent to f(y). If we drop the maximizationin Eq. (20), we obtain the boundf ð yÞVυy − f ðυÞ:

ð21Þ

Thus, for any given υ, we have a rigorous bound on f(y); we maythen optimize over υ to ﬁnd the optimal or tightest bound in a regionof interest. In higher dimensions, the strategy is the same, only now
the bounds generalize from lines to hyperplanes.
Appendix B. Proof of Theorem 1
We begin with Property 1. While it is possible to prove this result
using manipulations of L(Γ) in Γ-space, it is convenient to transform

654

D.P. Wipf et al. / NeuroImage 49 (2010) 641–655

the problem to an equivalent one in S-space. Speciﬁcally, the solution
S⁎ will be the global minimum of the dual optimization problem
S⁎ = lim

min g ðSÞð22Þ

Σe Y0 S:B = LS

Then we have
log jΣb j =DdcXlogðλi + eÞ + ðdb − Ddc Þ log e:

ð27Þ

i=1

This gives g(S) = O(1) + (db − Ddc) log ɛ.

■

where

g ðSÞJ min
C

" d
s
X



T −1
trace Si Ci Si + logjΣb j

Lemma 2. For any solution S such that B = LS, D will always satisfy
D ≥ da. The sources S that achieve equality are unique and satisfy S = Sgen.

#
ð23Þ

i=1

and Σb = Σɛ + LΓLT as deﬁned in the main text. The global and local
minima of Eq. (22) correspond with those of L(Γ) by straightforward
extension of results in Wipf and Nagarajan (2008). Therefore, we only
need to show that the global minimizer of Eq. (22) S⁎ satisﬁes
Property 1. For simplicity, we will assume that ∑ɛ = ɛI with ɛ→0 andthat spark (L) = db + 1 (i.e., every subset of db lead-ﬁeld columns arelinearly independent). The more general case is straightforward tohandle but complicates notation and exposition unnecessarily.
To begin we require two intermediate results for convenience andto avoid confusion, we deﬁne Sgen as the matrix of generating sourcesthat we would like to recover.

Proof. This result represents a simple extension of uncertainty
principles detailed in Cotter et al. (2005) and Donoho and Elad(2003). In particular, based on Lemma 1 in Cotter et al. (2005), if Sgensatisﬁes
da dc bðsparkðLÞ + rankðBÞ − 1Þ = 2VðsparkðLÞ + min ½dt ; da dc  − 1Þ = 2;ð28Þthen no other solution S can exist such that B = LS and D ≥ da.
Additionally, by directly applying results from Elad (2006), we ﬁndthat this condition will also hold for all S, except a set with measurezero, ifda dc b sparkðLÞ − 1:ð29Þ

Lemma 1. The function g(S) satisﬁes
g ðSÞ = Oð1Þ + ðdb − min ½db ; Ddc Þ log e;

ð24Þ

where D equals the cardinality of the set fSi : jjSi jj N 0g, i.e., D is thenumber of Si not equal to zero.2,3
In words this result states that g(S) will be O(1) unless Ddc b db, inwhich case g(S) will be dominated by the log ɛ term when ɛ is small.
Proof. Computing g(S) via Eq. (23) involves a minimization over two
terms. The ﬁrst term encourages each Γi to be large, the second
encourages each Γi to be small. Whenever a given Si = 0, the ﬁrst term
can be ignored and the associated Γi is driven to exactly zero by thesecond term. In contrast, for any Si ≠ 0, the minimizing Γi can never be
zero for any ɛ ≥ 0 or the ﬁrst term will be driven to inﬁnity. This amanifestation of the fact thatarg min 1 + logðx + eÞ N 0;
xz0 x

8ez0:

ð25Þ

Consequently, for any given S, the associated minimizing Γ will
necessarily have a matching sparsity proﬁle, meaning the indices ofzero-valued Si will align with zero-valued block-diagonal elementsin Γ.4
Whenever Ddc b db, the above analysis (and the assumption thatspark (L) = db + 1) ensures that the minimizing Σb will be full rankeven for ɛ = 0. This implies that g(S) = O(1) for essentially the same
reason that
min
xz0

1
+ logðx + eÞ = Oð1Þ:x

ð26Þ

In contrast, when Ddc b db, the minimizing Σb will become
degenerate when ɛ→0. Let λi denote the i-th non-zero eigenvalue ofLΓLT at the minimizing Γ. The spark assumption (coupled with theanalysis above) guarantees that there will be Ddc such eigenvalues.
2
Here we have adopted the notation f(x) = O(h(x)) to indicate that ∣f(x)∣ b C1∣h(x)∣for all x b C2, with C1 and C2 constants independent of x.
3
If S = Sgen, then by deﬁnition D = da.
4This point can be made more rigorous as shown in the ﬁrst author's PhD thesis, but
we omit lengthy details here.

■
Given these two results, the proof of Property 1 is simple. In thelimit as ɛ→0, Lemma 1 dictates that g(S) is minimized when D is
minimized. Lemma 2 then shows that D is uniquely minimized (withD ≥ da) when S = Sgen.
Property 2 is relatively easy to show by leveraging Theorem 2 from
Wipf and Nagarajan (2007). For our purposes, this result implies that
if the data covariance Cb satisﬁes Cb = Σb for some Γ, then the costfunction L(Γ) is unimodal and Cb = Σb at any minimizing solution.
When ɛ→0, the data covariance satisﬁes
Cb =


1 T
1
1
T
T
T
T
BB =
LSgen Sgen L = L
S S
L :
dt
dt
dt gen gen

ð30ÞGiven the conditions of Property 2, SgenSTgen will be zero except fordc × dc block-diagonal elements. Therefore, if Γ = 1 / dtSgenSTgen (whichis allowable given the speciﬁed block-diagonal structure of Γ), thenΣb = LΓLT = Cb and we have no local minima.
References
Baillet, S., Mosher, J., Leahy, R., Nov. 2001. Electromagnetic brain mapping. IEEE Signal
Process. Mag. 14–30.
Berger, J.O., 1985. Statistical Decision Theory and Bayesian Analysis, 2nd edition.
Springer-Verlag, New York.
Bolstad, A., Van Veen, B., Nowak, R., 2009. Space-time event sparse penalization formagneto-/electroencephalography. NeuroImage 46 (4), 1066–1081.
Boyd, S., Vandenberghe, L., 2004. Convex Optimization. Cambridge University Press.
Cotter, S.F., Rao, B.D., Engan, K., Kreutz-Delgado, K., April 2005. Sparse solutions to linear
inverse problems with multiple measurement vectors. IEEE Trans. Signal Process.
53 (5), 2477–2488.
Dale, A.M., Liu, A.K., Fischl, B.R., Buckner, R.L., Belliveau, J.W., Lewine, J.D., Halgren, E.,
April 2000. Dynamic statistical parametric mapping: combining fMRI and MEG for
high-resolution imaging of cortical activity. Neuron 26 (1), 55–67.
Darvas, F., Pantazis, D., Kucukaltun-Yildirim, E., Leahy, R.M., Sept. 2004. Mapping humanbrain function with MEG and EEG: methods and validation. NeuroImage 23 (Suppl. 1),
S289–S299.
Donoho, D.L., Elad, M., March 2003. Optimally sparse representation in general(nonorthogonal) dictionaries via ℓ1 minimization. Proc. Natl. Acad. Sci. 100 (5),2197–2202.
Elad, M., 2006. Sparse representations are most likely to be the sparsest possible.
EUROSIP J. Appl. Signal Process. 1–12.
Friston, K., Harrison, L., Daunizeau, J., Kiebel, S., Phillips, C., Trujillo-Barreto, N., Henson,
R., Flandin, G., Mattout, J., 2008. Multiple sparse priors for the MEG/EEG inverse
problem. NeuroImage 39 (1), 1104–1120.
Friston, K., Mattout, J., Trujillo-Barreto, N., Ashburner, J., Penny, W., 2006. Variationl freeenergy and the Laplace approximation. NeuroImage 34, 220–234.

D.P. Wipf et al. / NeuroImage 49 (2010) 641–655
Gorodnitsky, I., George, J., Rao, B., Oct. 1995. Neuromagnetic source imaging with
FOCUSS: a recursive weighted minimum norm algorithm. J. Electroencephalogr.
Clin. Neurophysiol. 95 (4), 231–251.
Huang, M., Dale, A., Song, T., Halgren, E., Harrington, D., Podgorny, I., Canive, J., Lewis, S.,
Lee, R., 2006. Vector-based spatial–temporal minimum ℓ1-norm solution for MEG.
NeuroImage 31 (3), 1025–1037.
Limpiti, T., Veen, B.V., Wakai, R., Sept. 2006. Cortical patch basis model for spatiallyextended neural activity. IEEE Trans. Biomed. Eng. 53 (9), 1740–1754.
Mattout, J., Phillips, C., Penny, W., Rugg, M., Friston, K., 2006. MEG source localization undermultiple constraints: an extended Bayesian framework. NeuroImage 30, 753–767.
Nagarajan, S.S., Attias, H.T., Hild, K.E., Sekihara, K., Sept. 2007. A probabilistic algorithm
for robust interference suppression in bioelectromagnetic sensor data. Stat. Med. 26(21), 3886–3910.
Nummenmaa, A., Auranen, T., Hämäläinen, M., Jääskeläinen, I., Lampinen, J., Sams, M.,
Vehtari, A., 2007. Hierarchical Bayesian estimates of distributed MEG sources:
theoretical aspects and comparison of variational and MCMC methods. NeuroImage
35 (2), 669–685.
Pascual-Marqui, R.D., 2002. Standardized low resolution brain electromagnetic tomography(sloreta): technical details. Methods Find. Exp. Clin. Pharmacol. 24 (Suppl. D), 5–12.
Phillips, C., Mattout, J., Rugg, M., Maquet, P., Friston, K., January 2005. An empirical Bayesiansolution to the source reconstruction problem in EEG. NeuroImage 24, 997–1011.
Phillips, C., Rugg, M., Friston, K., July 2002. Anatomically informed basis functions for
EEG source localization: combining functional and anatomical constraints. NeuroImage 13 (3), 678–695 pt. 1.
Ramírez R. and Makeig S., June 2006. “Neuroelectromagnetic source imaging using
multiscale geodesic neural bases and sparse Bayesian learning,” Human Brain
Mapping, 12th Annual Meeting, Florence, Italy.
Sahani, M., Nagarajan, S.S., 2004. Reconstructing MEG sources with unknown
correlations. Adv. Neural Inf. Process. Syst. 16.

655

Sarvas, J., 1987. Basic mathematical and electromagnetic concepts of the biomagnetic
inverse problem. Phys. Med. Biol. 32, 11–22.
Sato, M., Yoshioka, T., Kajihara, S., Toyama, K., Goda, N., Doya, K., Kawato, M., 2004.
Hierarchical Bayesian estimation for MEG inverse problem. NeuroImage 23, 806–826.
Sekihara, K., Nagarajan, S.S., 2008. Adaptive Spatial Filters for Electromagnetic Brain
Imaging. Springer.
Sekihara, K., Sahani, M., Nagarajan, S.S., 2005. Localization bias and spatial resolution ofadaptive and non-adaptive spatial ﬁlters for MEG source reconstruction. NeuroImage 25, 1056–1067.
Snodgrass, J.G., Corwin, J., 1988. Pragmatics of measuring recognition memory:
applications to dementia and amnesia. J. Exp. Psychol.: Gen. 17 (1), 34–50.
Uutela, K., Hamalainen, M., Somersalo, E., 1999. Visualization of magnetoencephalographic data using minimum current estimates. NeuroImage 10, 173–180.
Wipf D.P., 2006. “Bayesian methods for ﬁnding sparse representations,” PhD Thesis,
University of California, San Diego.
Wipf D.P. and Nagarajan S., June 2007. “Beamforming using the relevance vector
machine,” International Conference on Machine Learning.
Wipf, D.P., Nagarajan, S., 2008. A new view of automatic relevance determination.
Advances in Neural Information Processing Systems 20. MIT Press.
Wipf, D.P., Nagarajan, S., February 2009. A Uniﬁed Bayesian Framework for MEG/EEG
Source Imaging. NeuroImage 44 (3).
Wipf D.P. and Nagarajan S., 2009. “Iterative reweighted ℓ1 and ℓ2 methods for ﬁnding
sparse solutions,” Submitted.
Wipf, D.P., Ramírez, R.R., Palmer, J.A., Makeig, S., Rao, B.D., 2007. Analysis of empirical
Bayesian methods for neuroelectromagnetic source localization. Adv. Neural Inf.
Process. Syst. 19.
Zumer, J.M., Attias, H.T., Sekihara, K., Nagarajan, S.S., 2007. A probabilistic algorithm
integrating source locations and noise suppression for MEG and EEG data.
Neuroimage 37, 102–115.