Efficient hierarchical Bayesian inference for spatio-temporal regression models in neuroimaging
Authors: Ali Hashemi
Yijing Gao
Chang Cai
Sanjay Ghosh
Klaus Robert Muller
Srikantan Nagarajan
Stefan Haufe
Venue: NeurIPS 2021 Poster (2021)
Keywords: Multi-task Linear Regression
Hierarchical Bayesian Learning
Type-II Maximum-Likelihood Inference
Nonconvex Optimization
Majorization Minimization
Riemannian Geometry
Toeplitz
Circulant Embedding
Sparse Bayesian Learning
Sparse Models
Sparse Regression
Compressed Sensing
Sparse Denoising
Neuroimaging
EEG/MEG
Brain Source Imaging

Abstract: Several problems in neuroimaging and beyond require inference on the parameters of multi-task sparse hierarchical regression models. Examples include M/EEG inverse problems, neural encoding models for task-based fMRI analyses, and climate science. In these domains, both the model parameters to be inferred and the measurement noise may exhibit a complex spatio-temporal structure. Existing work either neglects the temporal structure or leads to computationally demanding inference schemes. Overcoming these limitations, we devise a novel flexible hierarchical Bayesian framework within which the spatio-temporal dynamics of model parameters and noise are modeled to have Kronecker product covariance structure. Inference in our framework is based on majorization-minimization optimization and has guaranteed convergence properties. Our highly efficient algorithms exploit the intrinsic Riemannian geometry of temporal autocovariance matrices. For stationary dynamics described by Toeplitz matrices, the theory of circulant embeddings is employed. We prove convex bounding properties and derive update rules of the resulting algorithms. On both synthetic and real neural data from M/EEG, we demonstrate that our methods lead to improved performance.

Efficient Hierarchical Bayesian Inference for Spatio-temporal Regression Models in Neuroimaging
Hashemi' Yijing Gao", Chang Cai Sanjay Ghosh
Roburt Müller2.5.677 Srikantan S. Nagarajan" and Stefan Haufe 1,89,,10
Unnertainty, Inverse Modeling and Machine Learning Group, Technische Universität Berlin, Germany. "aahhnee Learning Group, Technischc Universität Berlin, Germany.
"Department of Radiology and Biomedical Imaging, University California, San Francisco, USA.
National Engineering Research Center for E-Learning, Central China Normal University, China. BIFOLD Berlin Institute for the Foundations of Learning and Data, Berlin, Germany.
'Department of Artificial Intelligence, Korea University, South Korea.
Max Planck Institute for Informatics, Saarbrücken, Germany. 8 PPhysikalisch-Techniice Bundesanstalt, Berlin, Germany. 10 Bemstein Center for Computational Neuroscience, Berlin, Germany.
Charité Universilätsmediizn Berlin, Germany.
Abstract
Several problems in neuroimaging and beyond require inference the parameters multi-task sparse hierarchical regression models. Examples include M/EEG inverse problems, neural encoding models for task-based fMRI analyses, and climate science. In these domains, both the model parameters to be inferred and the measurement noise may exhibit a complex spatio-temporal structure. Existing work either neglects the temporal structure or leads to computationally demanding inference schemes. Overcoming these limitations, we devise novel flexible hierarchical Bayesian framework within which the spatio-temporal dynamics of model parameters and noise are modeled to have Kronecker product covariance structure. Inference in our framework is based on optimization and has guaranteed convergence properties. Our highly efficient algorithms exploit the intrinsic Riemannian geometry of temporal autocovariance matrices. For stationary dynamics described by Toeplitz matrices, the theory of circulant embeddings is employed. We prove convex bounding properties and derive update rules of the resulting algorithms. both synthetic and real neural data from M/EEG, demonstrate that methods lead to improved performance.
minimizzationn
Introduction
Probabilistic graphical models for regression problems, where both the model parameters and observation noise can have complex spatio-temporal correlation structure, are prevalent in neuroscience, neuroimaging and beyond 11 Example domains where these models are being
applied include Gaussian process inference 22] sparse signal recovery and multi-task learning (3]6),
array signal processing (7]. fMRI analysis
8 climate science [9] computer hardware monitoring
The focus in this paper will be on brain source imaging (BSI) [13), i.e. the reconstruction of brain source activity from non-invasive magneto- electroetoencephallggaaa. (M/EEG) sensors [14. Although a number of generative models that encompass spatio-temporal correlations have been proposed within the BSI regression problem, inference solutions have either imposed specific simplifications and model constraints [15] have ignored temporal correlations overall (5), and therefore not fully addressed the inherent spatio-temporal problem structure.
[0), and crain-computer interfaces (c.g.
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

We propose a novel and efficient algorithmic framework using Type-Il Bayesian regression that ex- plicitly considers the spatio-temporal covariance structure in both model coefficients and observation noise. To this end, we focus on the M/EEG inverse problem of brain source reconstruction, for which we adopt a multi-task regression approach and formulate the source reconstruction problem as a probabilistic generative model with Gaussian scale mixture priors for sources with spatio-temporal prior covariance matrices that are expressed as a Kronecker product of separable spatial and temporal covariances. Their solutions are constrained to the Riemannian manifold of positive definite (P.D.) temporal covariance matrices. Exploiting the concept of geodesic convexity on the Riemannian manifold of covariance matrices with Kronecker structure, and efficient majorization iiniiiinnniiiiiaiiionn (MM) optimization algorithms for model inference with provable convergence guarantees. In addition to deriving update rules for full-structural P.D. temporal covariance matrices, we also show that assuming a Toeplitz structure for the temporal covariance matrix results in computationally efficient update rules within the proposed MM framework.
then able to derive robust, fast
2 Spatio-temporal generative model
Let us consider a multi-task linear regression problem, mathematically formulated as
Y, g = LXg g Eg for g I,
(1)
in which forward matrix, L M maps set of coefficients source components X, E RXT to measurements Y, E RMXT, with independent Gaussian noise E, IR MxT. G denotes the number of sample blocks in the multi-task problem, while N, M, and T denote the number of sources or coefficients, the number of sensors or observations, and the number of time instants, respectively. The problem of estimating (Xg]g= given L and (Y,]g=) can represent an inverse problem in physics, a multi-task regression problem in machine learning, or multiple measurement vector (MMV) recovery problem in signal processing [8]
G
In the context of BSI, {Y, JG refers to M/EEG sensor measurement data, and (X,) J refers to brain source activity contributing to the observed sensor data. Here, G can be defined as the number of epochs, trials or experimental tasks. The goal of BSI is to infer the underlying brain activity from the M/EEG measurements given the lead field matrix L. In practice, L can be computed using discretization methods such the finite element method for given head geometry and known electrical conductivities using the stasi- approximation of Maxwell's s equations [19]. As the number of locations of potential brain sources is dramatically larger than the number of sensors, the M/EEG inverse problem is highly posep, which can be dealt with by incorporating prior assumptions. Adopting Bayesian treatment is useful in this capacity because it allows these assumptions to be made explicitly by specifying prior distributions for the model parameters. Inference can be performed either through Maximum-a-PPsteriorr (MAP) estimation (Type-I Bayesian learning) [20-22] or, when the model has unknown hyperparameters, through Type-1l Maximum- Likelihood (ML) estimation (Type-1l Bayesian learning) 23
Here we focus on full spatio-tempore Type-II Bayesian learning, which assumes a family of prior distributions p( pXX,, g ()) parameterized by set of hyperparameters +. We further assume that the spatio-tempordl correlation structure of the brain sources can be modeled by a Gaussian p(x,, T, B) ~ (0, Eo). where X, = vec(X.) E NT and where the covariance matrix Eg = T & B has Kronecker structure, meaning that the temporal and spatial correlations of the sources modeled independently through matrices and B, where @ stands for the Kronecker product. Note that the decoupling of spatial and temporal covariance is selooclysi plausible assumption as neuronal sources originating from different parts of the brain can be assumed to share a similar autocorrelation spectrum. In this paper, we set r diagg()) (7), 7 = (-,....77 7N] T, and B E 844,. where Sy_ denotes the set of positive definite matrices with size T T. Note that using a diagonal matrix I amounts to assuming that sources in different locations are independent. This assumption may be relaxed in our future work. On the other hand, modeling full temporal covariance matrix B amounts to assuming non- -tatiooio dynamics of the sources, which is appropriate in event-related experimental designs. To deal with the more stationary sources (ongoing brain activity), we will also constrain B to have Toeplitz structure later on. Based on these specifications, the prior distribution of the i-th brain source is modeled as p( (X Jo.11%: B) ~V(0, BB) for 1,..., N, where (X,)). denotes the i-th row of source matrix X,
2

---

Analogous to the definition of the sources, we here also model the noise to have temtoral structure. The matrix 3 E
IR MxT represents T time instances of zero-mean Gaussian noise with full covariance, Σ. = & Y, where cg vec(E, , E RT (0, Σe), and where A and T denote spatial and temporal noise covariance matrices, respectively. Here, we assume that the noise and sources share the same temporal structure, i.e., = B. We later investigate violations of this assumption in the simulation section. Analogous to the sources, we consider spatially independent noise characterized by diagonal spatial covariance A = diag(A), À = (AI....;>uu AM] We consider both the general, heleroscedastic, case in which the noise level may be different for each sensor as well as the homoscedastic case in which the noise level uniform.
For later use, we also define augmented versions of the source and noise covariance matrices as well as the lead field. Specifically, we set H and & so that H := [r, O: O, A], and definitions allow us to unify source and noise covariance parameters within the single variable H, which facilitates the concurrent estimation of both quantities. & now plays the same role as the lead field L, in the sense that it maps ny == [x Cg to the measurements; i.e., V, = @yg. Figure | (a) illustrates probabilistic graphical model summarizing the spatio-temporal generative model of our multi-task linear regression problem.
[L. I]. These
The MMV model Eq. ① be formulated equivalently in single measurement vector (SMV) form by vectorizing the spatio-temporal data matrices and using the Kronecker product as follows: Y, = Dx, + eg. where y, vec (Y, RMTx1 and D L Ir. Type-II Bayesian learning framework, the hyper-parameters of the spatio-temporal source model are optimized jointly with our case, these paper-parameters comprise the unknown source,
the model parameters (X, |
noise and temporal covariance matrices, i.e., e A....
3 Proposed method - full Dugh
The unknown parameters r, A, and B can be optimized in an alternating iterative process. Given ini- tial estimates, the posterior distribution of the sources is Gaussian of the form p(x (y T, A, B) (8,,, Ex), whose mean and covariance are defined as follows:
x, D (A B DT, ly, EnD
Σx = Σ - D D Σ, DEO,
(3)
where Ey = LLL + A, and where Σy = Ey @ B denotes the spatio-tempores statistical model
covariance matrix. The estimated posterior parameters Xg and Ex are then in turn used update r, A, and B as the minimizers of the negative log of the marginal likelihood p( (|T,, A, B), which is given by:
L (T, A. B) Tlog Σyl + M log B EirE, Y Y -Y,)
(4)
where |- denotes the determinant of matrix. A detailed derivation is provided in Appendix A Given the solution of the hyperparameters r, A, and B, the posterior source distribution obtained by plugging these estimates into Eqs. ( This process is repeated until convergence.
The challenge in this high -igh-iimensional inference problem is to find (locally) optimal solutions to Eq. 4 which is a non-convex cost function, in adequate time. Here we propose a novel efficient algorithm, which is able to do so. Our algorithm thus learns the full spatio-temporal correlation structure of sources and noise. Hashemi et al. (26] have previously demonstrated that majorization- minimization (MM) (27] is a powerful non-lineat optimization framework that can be leveraged to solve similar Bayesian Type-II inference problems. Here extend this work to spatio- temporal setting. Building on the idea of majorization-mininn minimization,, we construct convex surrogate functions that majorize Limm (T, A, B) in each iteration of the proposed algorithm. Then, we show the minimization equivalence between the constructed majorizing functions and Eq. (4). These results presented in theorems [ and 3 Theorems 2 and propose
efficient alternating optimization algorithm for solving Lkron (T. A, B), which leads to update rules for the spatial and temporal covariance matrices I' and B as well as the noise covariance matrix A.
Starting with the estimation of the temporal covariance based on the current source estimate, we can state the following theorem:
3

---

B
Spatial Covariance of Noise
Temporal Covariance 8'+
Spatial Covariance
Tangrst Space
Geometric Mean
S
E
P. M.naifoodd
M
(b)
Figure 1: (a) Probabilistic graphical model for the multi-task linear regression problem. (b) Geometric representation of the geodesic path between the pair of matrices {B², M Mime) on the P.D. manifold and the geometric mean between them, which is used to update B
Theorem 1. Optimizing the non-convex Type-Hl cost function in Eq. (4), Li (T, A, B), with respect to the temporal covariance matrix is equivalent to optimizing the following convex surrogate function, which majorizes Eq. (4)
Clime
(T A B, = tr ((B+)-1B) + tr time B-)),
(5)
where M is defined as
M time
1 MG
Σ (2) Y 9
(6)
and where Ik A and Σy denote the source, noise and statistical model covariance matrices at the k-th iteration, which are treated as constants when optimizing over B.
Proof. A detailed proof is provided in Appendix B
The solution of Eq. (5) can be obtained by taking advantage of the Riemannian geometry of the search space. The following theorem holds.
Theorem 2. The cost function Ctime (T* A A*,B B) is strictly geodesically convex with respect the P.D. manifold and its minimum with respect to B can be attained by iterating the following update rule until convergence:
B B--1 -- (B*)1/2
( (BB M
time (B*) (B²)- 1/2 (B*)
Proof. A detailed proof is provided in Appendix C
Remark 1. A geometric interpretation of the update rule 7 is that it finds the geometric mean between a spatially whitened version of the empirical temporal covariance matrix, M and the model temporal covariance matrix from the previous iteration, B A geometric representation of the geodesic path between the pair of matrices {Bk Mk (ine) ) on the P.D. manifold and the geometric mean between them, representing the update for B Bk+1 is provided in Figure (b).
We now derive an update rule for the spatial source and noise covariance matrices using an analogous approach. To this end, we first construct convex surrogate function that majorizes Lxvon (T, A. B) in each iteration of the optimization algorithm by considering augmented variables comprising both the source and noise covariances, i.e., H [T 0; A), and & LL. IJ. The following theorem holds.
Theorem 3. Minimizing the non-convex Type-Il cost function in Eq. (4), Liron (I, A, B), with respect the spatial covariance matrix H equivalent to minimizing the following convex surrogate function. which majorizes Eq. 4
espace cony
(T, A, B²)
(H, B.)
tr (E2)) -E
(H + tr (M "
(8)
conv

---

where MS is defined as
S
H T
*H Y (B*
MM
Maace (25) with
1 G
M space
Σ
TG 9=1
9
and where B denotes the temporal model covariance matrix estimated in the k-th iteration.
Proof. A detailed proof is provided in Appendix D
While, in principle, update rules for full-structured spatial source and noise covariances may be conceived, in analogy to the ones presented for the temporal covariance, we here restrict ourselves to the discussion of diagonal spatial covariances. Note, though, that such a choice of prior does not prohibit the reconstruction of parameters with more complex correlation structure. We have H diag (h),
Mi and Σy 4H T_ The update rule for H then takes simple form, as stated in the following theorem:
[Ti..... YN, .i..... ..........
Theorem 4. The cost function L H, B²) ) is convex in h, and its minimum with respect to h can be obtained according to the following flosed-form update rule, which concurrently estimates the scalar source variances and heteroscedastic noise variances:
H diagghh*))) 1841
g
for i = I...., N M where
(10)
g = diag M SN
(11)
z = diag( + ((=)))))0
(12)
Proof. A detailed proof can be found in Appendix
The final estimate of the posterior distribution (2 can be obtained by starting from any initial guess
for the model parameters, H = (To; A and B and iterating between update rules 7 (TO).
and ( until convergence. We call the resulting algorithm full Dugh. Convergence is shown in the following theorem:
Theorem 5. Optimizing the non-convex ML cost function Eq. 44 with alternating update rules for B and H = AT in Eq. 7 and Eq. (10) defines MM algorithm, which is guaranteed to converge to a stationary point.
Proof. A detailed proof can be found in Appendix G
Efficient inference for stationary temporal dynamics - - -thin Dugh
With full Dugh, could present general inference method for data with full temporal covariance structure. This algorithm may be used for non-stationary data, for which the covariance between two time samples depends the absolute position of the two samples within the studied time window. An example is the reconstruction of brain responses to repeated external stimuli from M/EEG (event- related electrical potential or magnetic field) data, where data blocks are aligned to the stimulus onset. In other realistic settings, however, stationarity can be assumed in the sense that the covariance between two samples only depends on their relative position (distance) each other but not to any external trigger. An example is the reconstruction of ongoing brain dynamics in absence of known external triggers. This situation can be adequately modeled using temporal covariance matrices with Toeplitz structure. In the following, we devise an efficient extension of full Dugh for that setting.
Temporal correlation has been incorporated in various brain source imaging models through different priors [see 13 28 and references therein]. According to [29] the temporal dynamics of neuronal populations can be approximated by a first order auto-regressive (AR(1)) model of the form xr, 1,... N; 1, 1,..., T with AR coefficient 3 E (-1,1) and innovation noise &n (t). It can be shown that the temporal correlation matrix corresponding to this AR
1) frr (1) + /1 (8ggn())
5

---

model has Toeplitz structure, Br,, 318 Eq. ( the set of Toeplitz matrices:
30]. Consequently, we now constrain the cost function
B R++
arg min (r((BB)BB))) + tr(M time B B)
(13)
BEB, H=
where set denotes the set vallev positive-definite Toeplitz matrices of size T x T.
In order to be able to derive efficient update rules for this setting, we bring B into diagonal form. The following proposition holds.
Proposition Let I& with T, be the circulant embedding of matrix BL, where BL is assumed to be the subset of B that guarantees that P is a real-valued circulant P.D. matrix. Then the Toeplitz matrix can be diagonalized as follows:
B = QPQ with Q=
and
(14)
Fi]m,
126t-111 (m-1)
(15)
VL
= diag(po- p1,.... -PL-11 is a diagonal matrix. The main diagonal coefficients Pf given the normalized discrete Fourier transform (DFT), represented by the linear operator F i. of the first row of B: F₁, (Note that a Toeplitz matrix can be represented by its first row or column, B
where == d
Proof.
This is a direct implication of the fact that Toeplitz matrices can be embedded into circulant matrices of larger size [3]] and the result that circulant matrices can be approximately diagonalized using the Fourier transform
32 Further details can be found in Appendix
The solution of Eq. 13 be obtained by direct application of Proposition The results is summarized in the following theorem:
Theorem 6. The cost function Eq. (13) is convex in p. and its minimum with respect to can be obtained by iterating the following closed-form update rule until convergence:
g
P1
for = 1,... L. where
(16)
%
g = diag (pag
pl Q (B*) M
time QP))
(17)
2 = diag Q4d (B*)-1Q)
(18)
Proof. A detailed proof is provided in Appendix
4.1 Efficient computation of the posterior
Given estimates of the spatial and temporal covariance matrices, can efficiently compute the posterior mean by exploiting their intrinsic diagonal structure.
Theorem 7, Given the diagonalization = QPQ" of the temporal correlation matrix and the eigenvalue decomposition LILT D. U of LrLT, where Dx diag(di,.... ......d...... the posterior mean is efficiently computed as
8 (- tPP Y U.) (U)
(19)
where 3 denotes the Hadamard product between the corresponding elements two matrices of equal size. In addition, II is defined as follows: [[] tpid for -,...... L and m = I,..., M.
Proof. A detailed proof is provided in Appendix
The resulting algorithm, obtained by iterating between update rules (16). (T and 19 is called thin
Dugh (as opposed to full Dugh introduced above).

---

5 Simulations
We present two sets of experiments to assess the performance of our proposed methods. In the first experiment, we compare the reconstruction performance of the proposed Dugh algorithm variants to that of Champagne [3 and two other competitive methods - cLORETA 34 and S-FLEX [35 for a range of SNR levels, numbers of time samples, and orders of AR coefficients. In the second experiment, we test the impact of model violations on the temporal covariance estimation. All experiments are performed using Matlab on a machine with a 2.50 GHz Intel(R) Xeon(R) Platinum 8160 CPU. The computational complexity of each method in terms of the average running time in units of seconds for 1000 iterations is as follows: Full Dugh: 67.094s, Thin Dugh: 62.289s, Champagne: 1.533s, eLORETA: 2. 2.653s, and S-FLEX: 20.963s. The codes are publicly available at https://giubb.com/AHHaaaasss -i/ Newrrll--u2P2
5.1 Pseudo-EEG signal generation and benchmark comparison
We simulate sparse set of No = 3 active sources placed random locations on the cortex. To simulate the electrical neural activity, we sample time series of length T E {10, 20, 50, 100) from a univariate linear autoregressive AR(P) process. We use stable AR systems of order P E 21, 7}. The resulting source distribution is then projected to the EEG sensors, denoted by Ysignal, using a realistic lead field matrix, L E I258 2004, We generate L using the New York Head model (36 taking into account the realistic anatomy and electrical tissue conductivities of an average human head. Finally, we add Gaussian white noise to the sensor space signal. Note that the simulated noise and source time courses do not share similar temporal structure here sources are modeled with a univariate autoregressive AR(P process while a temporally white Gaussian distribution is used for modeling noise. Thus, we could assess the robustness of our proposed method under violation of the model assumption that the temporal structure of sources and noise is similar. The resulting
noise matrix E- [e(1),...,. e(T)) is first normalized and then added to the signal matrix ysigual as follows: Ysignal + (1-a) Ysisnal
E. where a determines the signal-to-noise ratio (SNR) in
E| F
sensor space. Precisely, SNR is defined as SNR 20log10 (+/1-a). In this experiment the following values of a are used: a E [0.55, 0.65, 0.7, 0.8 which correspond to the following SNRs: SNR E 117,,,4, 7.4, 12) (dB). Interested readers can refer to Appendix K and (37) for further details on the simulation framework. We quantify the performance of all algorithms using the earth mover's 's distance (EMD) 21 38 and the maximal correlation between the time courses of the simulated and the reconstructed sources (TCE). Each simulation is carried out 100 times using different instances of X and E, and the mean and standard error of the mean (SEM) of each performance measure across repetitions is calculated.
In order investigate the impact of model violations on the temporal covariance estimation, we generate a random Gaussian source matrix, X E 3000000000 representing the brain activity of 2004 brain sources at 30 time instances for different numbers of trials G E 10, 20, 30, 40, 50. In all trials, sources are randomly sampled from a zero-mean normal distribution with spatio-temppral covariance matrix I @ B, where B E 30300 is either a full-structural random PSD matrix or a Toeplitz. matrix with 3 = 0.8. Gaussian noise E sharing the same temporal covariance with the sources is added to the measurements Y LX, so that the overall SNR is dB. We evaluate the accuracy of the temporal covariance reconstruction using Thin and full Dugh. The performance is evaluated using two measures: Pearson correlation between the original and reconstructed temporal covariance matrices, B and B, denoted by r(B, B). and the normalized mean squared error (NMSE) defined as: NMSE = | B BI] The similarity error is defined as: 1 - r(B, B). Note that NMSE measures the reconstruction at the true scale of the temporal covariance; while r(B, B) is inver-iivariant and hence only quantifies the overall structural similarity between simulated and estimated noise covariance matrices.
In Figure we show the source reconstruction performance (mean + SEM) of the four different brain source imaging schemes, namely thin Dugh, Champagne, eLORETA and S-FLEX. We notice that Dugh achieves superior performance in terms of EMD metric, whereas it competitive in terms of TCE. Note that since thin Dugh incorporates the temporal structure of the sources into the inference scheme, its performance with respect to EMD and TCE can be significantly improved by increasing the number of time samples. diqureedemoostratee from our two proposed spatio-temporal learning schemes, namely thin and full Dugh, indicated by cyan and magenta colors, respectively. The upper panel illustrates the reconstruction results for
the estimated temporal covariance matrix obtained
7

---

eLORETA -S-FEX
Champagne
Thin Dugh
0,4
0.4
0.4
EMD
EMD
EMD
0.2
0.2
0.2
0
2
6
10
50
100
2
4
6
SNB(dB)
# Time samples
AR order
0.2
0.2
0.2
TCE
TCE
TCE
0.1
0.1
0.1
0
0
0
2
6 10
50
100
2
6
SNR(dB)
#Time samples
AR order
Figure 2: Source reconstruction performance (mean + SEM) the four different brain source imaging schemes, namely Thin Dugh, Champagne, eLORETA, and S-FLEX, for data generated by realistic lead field matrix. Performance assessed for different settings including a wide range of SNRs, different numbers of time samples, and different AR model orders. Performance is evaluated in terms of the earth mover's distance (EMD) and time-course correlation error (TCE) between each simulated source and the reconstructed source with highest maximum absolute correlation.
setting where the ground-trun truth temporal covariance matrix has full structure, while the lower panel shows a case with Toeplitz temporal covariance matrix structure. It can be seen that full Dugh can better capture the overall structure of ground truth full-structure temporal covariance as evidenced by lower NMSE and similarity errors compared thin Dugh that is only able to recover Toeplitz matrix. As can be expected, the behavior is reversed in the lower-panel of Figure 3 where the ground truth temporal covariance matrix is indeed Toeplitz. Full Dugh, however, still provides reasonable performance in terms of NMSE and similarity error, even though it estimates full-structure temporal covariance matrix. Finally, the performance of both learning schemes are significantly improved by increasing the number of trials.
Analysis of real MEG and EEG data recordings
We now demonstrate the performance of the novel algorithms on two MEG datasets (see Figure and 5) and one EEG dataset (see Appendix D. All participants provided informed written consent prior to study participation and received monetary compensation for their participation. The studies were approved by the University of California, San Francisco Committee on Human Research. The MEG datasets from two female subjects comprised five trials of visual or auditory stimulus presentation, where the goal was to reconstruct cortical activity reflecting auditory and visual processing. prior work has shown success in reconstruction in such extreme low SNR data.
Figure 4 shows the reconstructed sources for auditory evoked fields (AEF) from representative subject using eLORETA, MCE, thin and full Dugh. In this case, we tested the reconstruction performance of all algorithms with the number of trials limited 5. As Figure 4 demonstrates, the reconstructions of thin and full Dugh both show focal sources at the expected locations of the auditory cortex. Limiting the number of trials as few as 5 does not negatively influence the reconstruction result of Dugh methods, while it severely affects the reconstruction performance of competing methods. For the visual evoked field (VEF) data in Figure 5) thin Dugh was able to reconstruct two sources in the visual cortex with their corresponding time-courses demonstrating an early peak 100 ms and a later peak around 200 ms (similar performance was observed for full Dugh).

---

Ground Truth
Thin Dugh
Full Dugh
-1
Fi-Full Dugh Thin Dugh
0.2
0.2
Full Structure
0.15
10
10
1-7(R,
NMSE
0.1
20
20
20
0.05
30
30
10
20
30
10
30
10
20
#Triils
# Trials
0.02
0.02
Toeplitz
0.8
1- .8)
10
10
NMSE
0.01
0.01
20
20
20
0.2
0
50
10
20
30
10
30
20
# Trials
Figure 3: Accuracy of the temporal covariance matrix reconstruction incurred by two different temporal learning approaches, namely thin (cyan) and full (magenta) Dugh, assuming Toeplitz and full temporal covariance structure, respectively. The ground-trute temporal covariance matrix (first column) has full covariance structure in the first row, and Toeplitz structure in the second row. Performance is assessed in terms of the Pearson correlation between the entries of the original and reconstructed temporal covariance matrices, B and B, denoted by r rBB, B). Shown is the similarity error r(B, B) (forth column). Further, the normalized mean squared error (NMSE) between B and B, defined as NMSE [B B / is reported (fifth column).
eLORETA
MCE
Thin Dugh
Full Dugh
trials
Figure 4: Analysis of auditory evoked fields (AEF) of one representative subject using eLORETA, MCE, thin and full Dugh. As can be seen, thin and full Dugh can correctly localize bilateral auditory activity Heschl's gyrus, which is the characteristic location of the primary auditory cortex, with as few as 5 trials. In this challenging setting, all competing methods show inferior performance.
Dugh
MCE
1st peak
0.
100
Timel1s)
200
250
50
100
160
200
209
2 peak
Figure 5: Spatial maps and corresponding time-courses reconstructed from five trials of visual stimulus presentation using Dugh and MCE, Distinct activity in visual cortex was reconstructed.

---

7 Discussion
Inverse modeling is a challenging problem that can be cast as a multi-task regression problem. Machine carning-based methods have contributed by systematically incorporating prior information into such models (e.g. [33 (T)))) While, so far, explicit models of spatio-temporal dynamies are rare, we contribute in this paper by deriving efficient inference algorithms for regression models with spatio-temporal dynamics in model parameters and noise. Specifically, we employ separable Gaussian distributions using Kronecker products of temporal and spatial covariance matrices. We assume sparsity and independent activations in the spatial domain, while the temporal covariance is modeled to have either full or Toeplitz structure. The proposed Dugh framework is encompassing efficient optimization algorithms for jointly estimating these distributions within a hierarchical Bayesian inference and MM optimization framework. Interestingly, we could theoretically prove convergence for the proposed update rules yielding estimates of the spatial and temporal covariances. In careful simulation studies, parameters in the model indeed leads to significantly improved reconstruction performance. Finally, the utility of our algorithms showcased in challenging real-world applications by reconstructing expected sources from real M/EEG data based on very few experimental trials.
have demonstrated that the inclusion of both spatial and temporal model
The M/EEG BSI problem has been previously framed as multi-task regression problem in the context of Type-I learning methods [5] 42). Bertrand et al. [5] proposed a method to extend the group Lasso class of algorithms to the multi- learning case. Their method shows promising results for large numbers of trials but is not competitive for smaller sample sizes, presumably due to not modeling the temporal structure of the sources. In contrast, Dugh is Type-Il method that learns the full spatio-temporal prior source distribution as part of the model fitting. Note that the assumption of spatially independent noise made here can be readily dropped, in ( A number of alternative approaches have been proposed to estimate the spatio-tempore correlation structure of the noise 44.. These works, however, do not estimate the noise characteristies as part of the source reconstruction problem on the same data but require separate noise recordings. Our proposed algorithm substantially differs in this respect, as it learns the full noise covariance jointly with the brain source distribution. This joint estimation perspective is very much in line with the end-to-en to- learning philosophy as opposed to a step-wise independent estimation process that can give rise to стror accumulation.
With respect to the limits of our Dugh framework, we would like to note that using the same temporal correlation prior for noise and sources is a potentially restricting assumption in our modeling. It is made here achieve tractable inference due the fact that the spatio-temporal statistical model covariance matrix can be formulated in a separable form, i.e., Σy Ey @ B. Although we have demonstrated empirically that the reconstruction results of our proposed learning schemes are fairly robust against violating this assumption. this constraint may be further relaxed by exploiting eigenvalue decomposition techniques presented in [2] 49]. Another potentially limiting assumption in our model is to assume Gaussian distributions for the sources and noise. Although Gaussian priors are commonly justified, the extension of our framework to heavy-tailed noise distributions, which are more robust to outliers, is a direction of future work. Regarding the societal impact of this work, we note that we intend to solve inverse problems that have non-unique solutions, which depend on explicit or implicit assumptions of priors. These have applications in basic, clinical, and translational neuroscience imaging studies. Our use of hierarchical empirical Bayesian framework allows for explicit specifications of the priors in our model that are learned from data, and a more clear interpretation of the reconstructions with reduced bias. Nevertheless, in our algorithms, we assume sparse spatial priors, and in scenarios where this assumption embodied in the priors is incorrect, the resulting reconstructions will be inaccurate. Users and neuroscientists will need to be cognizant of these issues.
In conclusion, we could derive novel flexible hierarchical Bayesian algorithms for multi-task regres- sion problems with spatio-temporal covariances. Incorporating prior knowledge, namely, constraining the solutions to Riemannian manifolds, and using ideas of geodesic convexity, circulant embeddings, we derive inference update rules and prove convergence guarantees. The proposed Dugh algorithms show robust and competitive performance both on synthetic and real neural data from M/EEG recordings and thus contribute to a well-founded solution to complex inverse problems in neuroscience.
and mngorization-minnnizzze
10

---

Acknowledgments and Disclosure of Funding
This result is part of project that has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 758985). AH acknowledges scholarship support from the Machine Learning/Intellieent Data Analysis research group at Technische Universität Berlin. He further wishes to thank the Charité Universitätsmedzzin Berlin, the Berlin Mathematical School (BMS), and the Berlin Mathematics Research Center MATH+ for partial support. CC was supported by the National Natural Science Foundation of China under Grant 62007013. KRM was funded by the German Ministry for Edu- cation and Research as BIFOLD Berlin Institute for the Foundations of Learning and Data (ref. 011S18025A and ref. 011S18037A) A), and the German Research Foundation (DFG) Math+: Berlin Mathematics Research Center (EXC 2046/1. project-ID 390685689), Institute of Information & Communications Technology Planning & Evaluation (ITTP) grants funded by the Korea Government (No. 2019-0-00072 Artificial Intelligence Graduate School Program, Korea University). SSN was funded in part by National Institutes Health grants (R0IDC004855, ROIEB0227 17, R01DCI76960, ROI ROINS 100440, R01AG062196 and R01DC013979), University of California MRPI MRP-17-454755, the US Department Defense grant (W8IXWH--1311-099))
References
[1] Richard M Shiffrin, Michael D Lee, Woojae Kim, and Eric-Jan Wagenmakers. A survey of model evaluation approaches with a tutorial on hierarchical Bayesian methods. Cognitive Science, 32(8):1248-12) 2008.
[2] Barbara Rakitsch, Christoph Lippert, Karsten Borgwardt, and Oliver Stegle. It is all in the noise: Efficient multi-task Gaussian process inference with structured residuals. In Proceedings of the 26th International Conference on Neural Information Processing Systems- Volume 1, pages 1466 1674, 2013.
[3] Emmanuel J Candès, Justin K Romberg, and Terence Tao. Stable signal recovery from incomplete and inaccurate measurements. Communications on pure and applied mathematics, 59(8):1207-1223 2006.
[4] David L Donoho. Compressed sensing. IEEE Transactions on information theory, 52(4): 1289-1306, 2006.
[5] Quentin Bertrand, Mathurin Massias, Alexandre Gramfort, and Joseph Salmon. Handling correlated and repeated measurements with the smoothed multivariate square-root Lasso. In Advances in Neural Information Processing Systems, pages 3959-3970, 2019.
[6] Jérôme-Alext Chevalier, Alexandre Gramfort, Joseph Salmon, and Bertrand Thirion. Statisti- cal control for spatio-temporal MEG/EEG source imaging with desparsified multi-task lasso. arXiv preprint 20xii:2009.143100 2020.
[7] Dmitry Malioutov, Müjdat Cetin, and Alan S Willssky A sparse signal reconstruction perspec- tive for source localization with sensor arrays. IEEE Transactions on Signal Processing, 53(8): 2005.
[8] Ming Bo Cai, Michael Shvartsman, Anqi Wu, Hejia Zhang, and Xia Zhu. Incorporating structured assumptions with probabilistic graphical models in fMRI data analysis. Neuropsy- chologia, page 107500, 2020.
[9] S Beirle, U Platt, Mr Wenig, and T Wagner. Weekly cycle of NO2 by GOME measurements: a signature of anthropogenic sources. Atmospherie Chemistry and Physics, 3(6):2225 2232, 2003.
[10] Juri Ranieri, Alessandro Vincenzi, Amina Chebira, David Atienza, and Martin Vetterli. Near- optimal thermal monitoring framework for many-core systems-on-chip. IEEE Transactions on Computers, 64(1 (4(11):31773309 2015.
[11] Steven Lemm, Benjamin Blankertz, Gabriel Curio, and K-R Muller. Spatio-spectral filters for improving the classification of single trial eeg. IEEE transactions on biomedical engineering, 1541-1554) 2005.

---

[12] Daniel Bartz and Rlaur-Robert Müller. Covariance shrinkage for autocorrelated data. In Z. Ghahramani, M. Welling. C. Cortes, N. Lawrence, and K. Q. Weinberger, editors, Advances in Neural Information Processing Systems, volume 27. Curran Associates, Inc., 2014.
[13] Wei Wu, Srikantan Nagarajan, and Zhe Chen. Bayesian machine learning: EEG MEG signal processing measurements. IEEE Signal Processing Magazine, 33(1): 14-36, 2016.
[14] Sylvain Baillet, John C Mosher, and Richard M Leahy, Electromagnetic brain mapping. IEEE Signal Processing Magazine, 18(6):14-3)) 30, 2001.
[15] Mingho Cai, Nicolas W Schuck, Jonathan W Pillow, and Yael Niv. A Bayesian method for reducing bias in neural representational similarity analysis. In Advances in Neural Information Processing Systems, pages 4951 1951,4959, 2016.
[16] Michael Shvartsman, Narayanan Sundaram, Mikio Aoi, Adam Charles, Theodore Willke, and Jonathan Cohen. Matrix-normal models for fMRI analysis. In International Conference on Artificial Intelligence and Statistics, pages 1914-1923. PMLR, 2018.
[17] Ming Bo Cai, Nicolas W Schuck, Jonathan W Pillow, and Yael Niv. Representational structure or task structure? Bias in neural representational similarity analysis and Bayesian method for reducing bias. PLoS computational biology. 15(5)ee100629) 1006299. 2019.
[18] Shane Cotter, Bhaskar Rao, Kjersti Engan, and Kenneth Kreutz-Delgado. Sparse solutions to linear inverse problems with multiple measurement vectors. IEEE Transactions on signal processing, 53(7):2477-2488 2005.
[19] Matti Hämäläinen, Ritta Hari, Risto J Ilmoniemi, Jukka Knuutila, and Olli V Lounasmaa. Magmetoonccephalographyy theory, instrumentation, and applications noninvasive studies of the working human brain. Reviews of modern Physics, 65(2):413, 1993.
[20] Roberto D Pascual- Pascual-Marqui, Christoph M Michel, and Dietrich Lehmann. Low resolution electromagnetic tomography: new method for localizing electrical activity in the brain. International Journal of psychophysiology, 18(1):49-66 65, 1994.
[21] Stefan Haufe, Vadim V Nikulin, Andreas Ziche, Klaus-Robert Müller, and Guido Nolte. Combining sparsity and rotational invariance in EEG/MEG source reconstruction. Neurolmage. 42(2):726-738, 2008.
[22] Alexandre Gramfort. Matthicu Kowalski, and Matti Hämäläinen. Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods. Physics medicine and biology, 57(7): 2012.
[23] Sebastian Mika, Gunnar Rätsch, and Klaus-Robert Müller. A mathematical programming approach to the kernel fisher algorithm. Advances in neural information processing systems, 13:591-597, 2000.
[24] Michael E Tipping. Sparse bayesian learning and the relevance vector machine. Journal of Machine Learning Research, ((June):211-2 1-244, 2001.
[25] David Wipf and Srikantan Nagarajan. unified Bayesian framework for MEG/EEG source imaging. Neurolmage, 44(3):947 2009.
[26] Ali Hashemi, Chang Cai, Gitta Kutyniok, Klaus-Robert Müller, Srikantan Nagarajan, and Stefan Haufe. Unification of sparse Bayesian learning algorithms for electromagnetic brain imaging with the majorization minimization framework. Neurolmage, 239:118309, 2021.
[27] Ying Sun, Prabhu Babu, and Daniel P Palomar. aisiorzzttio-mmmimizzzzn algorithms in signal processing, communications, and machine learning. IEEE Transactions on Signal Processing, 65(3):794-8) 816,, 2017.
[28] Elvira Pirondini, Behtash Babadi, Gabriel Obregon-Henno, Camilo Lamus, Wasim Q Malik, Matti S Hamalainen, and Patrick I. Purdon. Computationall--effiiees algorithms for sparse, dynamic solutions to the eeg source localization problem. IEEE Trunsactions on Biomedical Engineering, 2017.
12

---

[29] Camilo Lamus, Matti S Hämäläinen, Simona Temereanca, Emery N Brown, and Patrick L Pur- don. spatiotemporal dynamic distributed solution to the MEG inverse problem. Neurolmage. 63(2):894-909) 909, 2012.
[30] Zhilin Zhang and Bhaskar D Rao. Sparse signal recovery with temporally correlated source vectors using sparse Bayesian learning. IEEE Journal of Selected Topics in Signal Processing. 5(5):912- 9269 2011.
[31] Amir Dembo, Colin I. Mallows, and Lawrence A Shepp. Embedding nonnegative definite Toeplitz matrices in nonnegative definite circulant matrices, with application to covariance estimation. IEEE Transactions on Information Theory, 1566):1206-1212, 1989.
[32] Ulf Grenander and Gabor Szegö. feeplitz,oorm and their applications. Univ of California Press, 1958.
[33] David P Wipf, Julia Owen, Hagai T Attias, Kensuke Sckihara, and Srikantan S Nagarajan. Robust Bayesian estimation of the location, orientation, and time course of multiple correlated neural sources using MEG. Neurolmage, 49(1):641-655, 2010.
[34] Roberto D Marquu.. Discrete, 3D distributed, linear imaging methods of electric neuronal activity, part I: exact, zero error localization. arXiv preprint arXiv:0710.3341, 2007.
[35] Stefan Haufe, Ryota Tomioka, Thorsten Dickhaus, Claudia Sannelli, Benjamin Blankertz, Guido Nolte, and -lbbrr Müller. Large-scale EEG/MEG source localization with spatial flexibility.. Neurolmage, 54(2):851-855) 2011.
[36] Yu Huang, Lucas C Parra, and Stefan Haufe. The New York head-a precise standardized volume conductor model for EEG source localization and tes targeting. Neurolmage, 140: 150-162, 2016.
[37] Stefan Haufe and Arne Ewald. A simulation framework for benchmarking EEG- based brain connectivity estimation methodologies. Brain topography, pages 1-18, 2016.
[38] Y. Rubner, C. Tomasi, and L. J. Guibas. The earth mover's distance metric for image retrieval. Int J Comput Vision, 40:99-121, 2000.
[39] Michael T McCann, Kyong Hwan Jin, and Michael Unser. Convolutional neural networks for inverse problems in imaging: A review. IEEE Signal Processing Magazine, 34(6):85-95, 2017.
[40] Simon Arridge, Peter Maass, Ozan Öktem, and Carola-Bibiane Schönlieb, Solving inverse problems using data-driven models. Acta Numerica, 28:1-174, 2019.
[41] Tatiana A Bubba, Gitta Kutyniok, Matti Lassas, Maximilian Maerz, Wojciech Samek, Samuli Siltanen, and Vignesh Srinivasan. Learning the invisible: hybrid deep learning-shearlet framework for limited angle computed tomography. Inverse Problems, 35(6):064002, 2019.
[42] Mathurin Massias, Olivier Fercoq, Alexandre Gramfort, and Joseph Salmon. Generalized concomitant multi-task lasso for sparse multimodal regression. In International Conference on Artificial Intelligence and Statistics, pages 998- 10081 2018.
[43] Ali Hashemi, Chang Cai, Klaus Robert Müller, Srikantan Nagarajan, and Stefan Haufe. Joint learning of full-structure noise hierarchical Bayesian regression models. 2021. URL
net:forum?id=y/aae _Bxddu
[44] Hilde M Huizenga, Jan C De Munck, Lourens Waldorp, and Raoul PPP Grasman. Spa- tiotemporal EEG/MEG source analysis based on a parametric noise covariance model. IEEE Transactions on Biomedical Engineering, 49(6):533-533, 2002.
[45] Jan Casper De Munck, Hilde M Huizenga, Lourens J Waldorp, and RA Heethaar. Estimating stationary dipoles from MEG/EEG data contaminated with spatially and temporally correlated background noise. IEEE Transactions Signal Processing, 1557 2002.
13

---

[46] Fetsje Bijma, Jan C De Munck, Hilde M Huizenga, and Rob M Heethaar. A mathemati- cal approach to the temporal stationarity of background noise in MEG/EEG measurements. Neurolmage, (0(1):233-243 2003.
[47] Jan Casper De Munck, Fetsje Bijma, Pawel Gaura, Cezary Andrzej Sieluzycki, Maria Inês Branco, and Rob M Heethaar. maximum-likelihodd estimator for trial-to-trial variations in noisy MEG/EEG data sets. IEEE Transactions on Biomedical Engineering, 51 51(12):2123-2128 2004.
[48] Sung C Jun, Sergey M Plis, Doug M Ranken, and David M Schmidt. Spatiotemporal noise covariance estimation from limited empirical magnetoencephallgraphii data. Physics in Medicine & Biology, 51(21):5549, 2006.
[49] Anqi Wu, Stan Pashkovski, Sandeep R Datta, and Jonathan W Pillow. Learning a latent manifold of odor representations from neural responses in piriform cortex. Advances in Neural Information Processing Systems, 31:5378-5533 5388,, 2018.
[50] Yonina C Eldar and Moshe Mishali. Robust recovery of signals from a structured union of subspaces. IEEE Transactions information theory, 55(11):5302-55) 5316, 2009.
[51] Yuzhe Jin and Bhaskar D) Rao. Support recovery of sparse signals in the presence of multiple measurement vectors. IEEE Transactions on Information Theory, 59(5):3139-3157, 2013,
[52] Piya Pal and PP Vaidyanathan. Pushing the limits of sparse support recovery using correlation information. IEEE Transactions on Signal Processing, 63(3):711-726, 2014.
[53] Rajendra Bhatia. Positive definite matrices, volume 24. Princeton University Press, 2009.
[54] Pourya Zadeh, Reshad Hosseini, and Suvrit Sra. Geometric mean metric learning. In Interna- tional conference machine learning, pages 2464-2477 2016.
[55] Ami Wiesel, Teng Zhang, et al. Structured robust covariance estimation. Foundations and Trendss) in Signal Processing, 8(3): 2015.
[56] Athanase Papadopoulos. Metric spaces, convexity and nonpositive curvature, volume European Mathematical Society, 2005.
[57] Tamas Rapesak. Geodesic convexity in nonlinear optimization. Journal of Optimization Theory and Applications, 69(1):169-183, 1991.
[58] Aharon Ben-Tal. On generalized means and generalized convex functions. Journal of Opti- mization Theory and Applications, 21(1):1-13, 1977.
[59] Leo Liberti. On class of nonconvex problems where all local minima are global. Publications de I' I'Institut Mathémathique, 76(90))101-109 2004.
[60] Diethard Ernst Pallaschke and Stefan Rolewicz. Foundations of mathematical optimization: convex analysis without linearity, volume 388. Springer Science & Business Media, 2013.
[61] Silvere Bonnabel and Rodolphe Sepulchre. Riemannian metric and geometric mean for positive semidefinite matrices of fixed rank. SIAM Journal on Matrix Analysis and Applications, 31))
1055-1070, 2009.
[62] Maher Moakher. differential geometric approach to the geometric mean of symmetric positive-definit matrices. SIAM Journal on Matrix Analysis and Applications, 26(3):735-747, 2005.
[63] Suvrit Sra and Reshad Hosseini. Conic geometric optimization on the manifold of positive definite matrices. SIAM Journal on Optimization, 25(1):713-739, 2015.
[64] Nisheeth K Vishnoi. Geodesic convex optimization: Differentiation manifolds, geodesics, and convexity, arXiv preprint arXiv: 06666.6663,. 2018.
[65] Constantin Niculescu and Lars-Erik Persson. Convex functionn and their applications. Springer, 2006.
14

---

[66] Jason V Davis, Brian Kulis, Prateek Jain, Suvrit Sra, and Inderjit S Dhillon. Information- theoretic metric learning. In Proceedings of the 24th international conference on Machine learning, pages 209-216. ACM, 2007.
[67] Stephen P Boyd and Lieven Vandenberghe. Convex optimization. Cambridge university press, 2004.
[68] Chang Cai, Ali Hashemi, Mithun Diwakar, Stefan Haufe, Kensuke Sekihara, and Srikantan S Nagarajan. Robust estimation of noise for electromagnetic brain imaging with the Champagne algorithm. Neurolmage, 225:1172:0 117411 2021.
[69] David Hunter and Kenneth Lange. tutorial on MM algorithms. The American Statistician, 58(1):30-37, 2004.
[70] Meisam Razaviyayn, Mingyi Hong, and Zhi-Quan Luo. A unified convergence analysis of block successive minimization methods for nonsmooth optimization. SIAM Journal on Optimization, 23(2): 126- 1156,,113,0000 2013.
[71] Matthew W Jacobson and Jeffrey A Fessler. An expanded theoretical treatment of iteration- dependent minimize algorithms. IEEE Transactions on Image Processing, 16(10); 2411-2422, 2007.
[72] Tong Tong Wu, Kenneth Lange, al. The MM alternative to EM. Statistical Science, 25(4): 492-505, 2010.
[73] Eric Mjolsness and Charles Garrett. Algebraic transformations of objective functions. Neural Networks, ((6):651-66)) 1990.
[74] Alan L Yuille and Anand Rangarajan. The concave-convex procedure. Neural computation, 15(4):9159936, 2003.
[75] Thomas Lipp and Stephen Boyd. Variations and extension the convex-concave procedure. Optimization and Engineering, 17(2):663-877, 2016.
[76] Dylan Fagot, Herwig Wendt, Cédric Févotte, and Paris Smaragdis, Majjorzatioon-ninnnizaaooo algorithms for convolutive NMF with the beta-divergence. In ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 8202- 8206. IEEE, 2019.
[77] Sandeep Kumar, Jiaxi Ying. José Vinícius de Miranda Cardoso, and Daniel P Palomar. A unified framework for structured graph learning via spectral constraints. Journal of Machine Learning Research, 21(2)):1666 2020.
[78] Yiyong Feng, Daniel P Palomar, et al. signal processing perspective on financial engineering. Foundations and Trends in Signal Processing, 9(1-2):1-23) 2016.
[79] Ranjitha Prasad, Chandra R Murthy, and Bhaskar Rao. Joint channel estimation and data detection in MIMO-OFDM systems: A sparse Bayesian learning approach. IEEE Transactions on Signal Processing, 63(20):5369-5338 2015.
[80] Peter Gerstoft, Christoph F Mecklenbräuker, Angeliki Xenaki, and Santosh Nannuru. Mul- tisnapshot sparse Bayesian learning for DOA. IEEE Signal Processing Letters, 23(10): 1469-1473, 2016.
[81] Sacid Haghighatshoar and Giuseppe Caire. Massive MIMO channel subspace estimation from low-dimensional projections. IEEE Transactions on Signal Processing, 65(2):303-318, 2017.
[82] Mahdi Barzegar Khalilsarai, Tianyu Yang, Saeid Haghighatshoar, and Giuseppe Caire. Struc- tured channel covariance estimation from limited samples in Massive MIMO. In ICC 2020- 2020 IEEE International Conference on Communications (ICC), pages 1-7. IEEE, 2020.
[83] Alexander Fengler, Giuseppe Caire, Peter Jung, and Saeid Haghighatshoar, Massive MIMO unsourced random access. arXiv preprint arXiv: 0000808 2019.
15

---

[84] Alexander Fengler, Saeid Haghighatshoar, Peter Jung, and Giuseppe Caire. Non- Non-Bayesian activity detection, large-scale fading coefficient estimation, and unsoureed random access with a massive MIMO receiver. arXiv preprini arXiv: 1266, 2019.
[85] Martin Luessi Matti S Hämäläinen, and Victor Solo. Sparse component selection with application to MEG source localization. In 2013 IEEE 10th International Symposium on Biomedical Imaging, pages 556-559. IEEE, 2013.
[86] Ali Hashemi and Stefan Haufe. Improving EEG source localization through apatio-temports sparse Bayesian learning. In 2018 26th European Signal Processing Conference (EUSIPCO), pages 1935- 1939. IEEE, 2018.
[87] Yousra Bekhti, Felix Lucka, Joseph Salmon, and Alexandre Gramfort. A hierarchical Bayesian perspective on M/EEG source imaging. Inverse Problems, 3488)008500,, 2018.
for non-convex sparse regression: application to
[88] Prabhu Babu, MELT Sag---mmmmmmm--iieiiiiooooooo estimation of low-rank Toeplitz covariance matrix. IEEE Signal Processing Letters, 23(11): 2016.
[89] Petre Stoica, Randolph L Moses, et al. Spectral analysis of signals, volume 1. Pearson Prentice Hall Upper Saddle River, NJ, 2005.
[90] Robert M Gray et al. Toeplitz and circulant matrices: A review. Foundations and Trends in Communications and Information Theory, 2(3):155-239, 2006.
[91] Ying Sun, Prabhu Babu, and Daniel P Palomar. Robust estimation of structured covariance matrix for heavy-tailed elliptical distributions. IEEE Transactions on Signal Processing, 64 (14):3576-359) 2016.
[92] Yunus Saatçi. Scalable inference for structured gaussian process models. Ph.D. thesis, University of Cambridge Cambridge, UK., 2012.
[93] Arno Solin, Pasi Jylänki, Jaakko Kauramäki, Tom Heskes, Marcel AJ van Gerven, and Simo Särkkä. Regularizing solutions the MEG inverse problem using space-time separable covariance functions. arXiv preprint (004:04600004311,, 2016.
[94] R. Oostenveld and P. Praamstra. The five percent electrode system for high-resolution EEG and ERP measurements. Clinical Neurophysiology: Official Journal of the International Federation of Clinical Neurophysiology, 112(4):713-719, 2001. ISSN 1 1388-2457. doi: 10. (0016//138--245000000 2477500)
[95] Alexandre Gramfort. Mapping, timing and tracking cortical activations with MEG and EEG: Methods and application to human vision. PhD thesis, Ecole nationale supérieure des ielecommunicationn-Pliee ENST, 2009.
[96] Kanta Matsuura and Yoichi Okabe. Selective minimum-norm solution of the biomagnetic inverse problem. IEEE Transactions Biomedical Engineering, 42(6))608-611 1995.
[97] Sarang Dalal, JM Zumer, V Agrawal, KE Hild, Kensuke Sekihara, and SS Nagarajan. NUT- MEG: a neuromagnetic source reconstruction toolbox. Neurology & Clinical Neurophysiology: NCN, 2004:52, 2004.
[98] Sarang Dalal, Johanna M Zumer, Adrian G Guggisherg, Michael Trumpis, Daniel DE Wong, Kensuke Sckihara, and Srikantan S Nagarajan. MEG/EEG source reconstruction. statistical evaluation, and visualization with NUTMEG. Computational Intelligence and Neuroscience, 2011, 2011.
[99] Julia P Owen, David P Wipf, Hagai T Attias, Kensuke Sekihara, and Srikantan S Nagarajan. Performance evaluation of the Champagne source reconstruction algorithm on simulated and real M/EEG data. Neuroimage, (0(1):305-325 2012.
[100] Chang Cai, Mithun Diwakar, Dan Chen, Kensuke Sekihara, and Srikantan S Nagarajan. Robust empirical Bayesian reconstruction of distributed sources for electromagnetic brain imaging. IEEE Transactions on Medical Imaging, 39(3):567-577, 2019. 2019.
16


Reviews:
[{"review_id": "z799wk0sCT5", "reviewer": "Reviewer_exED", "summary": "Inverse modeling is at the heart of medical imaging techniques, especially in neuroimaging. Many techniques have been developed in recent years, first based on PDE techniques and, more recently, by machine learning approaches. The addition of a priori knowledge to these techniques allows a better understanding of the studied phenomenon. Therefore, it is necessary to develop robust methods combining modeling, a priori knowledge, while controlling the computational cost. This article proposes to derive an efficient inference algorithm with spatio-temporal dynamics in model parameters and noise. The spatial and temporal duality is realized by using the product of \nKronecker product. They assume sparse spatial covariance matrices, while the temporal covariance is modeled to have either full or Toeplitz structure. They then define a Bayesian model whose parametric estimation is performed by a Maximization-Minimization algorithm. The proposed results are motivated by demonstrations (detailed in appendix) and numerical experiments.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The article is well constructed, as are the appendices, which provide many insights.\n\nThe demonstrations, postponed to the appendix, are pleasant to follow. Moreover, for those I have reread, at least, they seem correct unless I am mistaken. However, they are essentially based on Riemannian geometry, which may reduce their scope. The authors have nevertheless made an effort to give the essential definitions in appendix C (even if I doubt that this is sufficient for a \"neophyte\" in this field).\nI particularly appreciated that each theorem stated in the article is precisely referred to its proof.\n\nThe chosen approach seems innovative. However, I wonder what sense it makes to assume the same temporal correlation for noise and sources in terms of modeling. And to what extent is this assumption realistic, encountered in practice, reasonable? I am aware that this makes the computation easier, as claimed by the authors in the discussion. Are the authors convinced that this assumption can also be relaxed theoretically and not just empirically as is currently the case?\n\nSmall remarks: prefer to use a git repository rather than a folder on a dropbox, which seems to me less perennial.\nAlso, Figure 1 is not readable when printed in black and white: I suggest that authors use different symbols for each node of the curves and include them in the legend.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The article discusses methods to improve medical imaging techniques. Unless medicine is misused, it seems to me that it cannot have a negative social impact (but I could be wrong).", "review_text": "", "extra_scores": {"ethical_concerns": "I did not identify any ethical issues in this document. A call for volunteers was required for this work. Participants had to provide written informed consent and received monetary compensation. ", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "KlKKTXQ4p6X", "reviewer": "Reviewer_KaNL", "summary": "The paper proposes a Bayesian framework for source estimation in M/EEG. The Bayesian framework explicitly estimates spatial and temporal covariance matrices for both signal and noise, before estimating the posterior mean. A majorization-minimization algorithm is proposed to estimate the model, and detailed proofs are provided, including a convergence proof. The model is benchmarked against three source-estimation methods on simulated signals. Finally, the model is used in empirical recordings with a challenging low number of trials.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "**Originality**: The problem of M/EEG source estimation is a well-known and challenging problem. This paper propose a new Bayesian model to solve this problem. The novelty resides in the explicit spatial and temporal priors, where the covariances are estimated jointly with the model. Related work seems to be well cited and discussed.\n\n**Quality**: The model is sound and technically advanced, and the methods used are appropriate. The technical derivations (in supplementary) are thorough and well described. The proofs look reasonable, though I am not familiar with this sort of work, and I have only skimmed through them. The simulation benchmark is short, and only consider a subset of the competing methods, but results are convincing. The empirical result seems impressive. The discussion of the pros and cons of this method is appreciated, although the discussion against type-I methods could have been supported by a benchmark.\n\n**Clarity**: The submission is clearly written and well organized. A large part of the work is located in the long supplementary document, with detailed proofs and discussions. The experimental part is well described, although less well polished than the technical derivations. The figure quality is mediocre. Importantly, the empirical results (Figure 3) are not interpretable without the help of Figure 4-5-6 (all in supplementary). Improving the organization of these figures could help the reader understand the interest of these results.\n\n**Significance**: The model is sounds and seems to be more robust than other models. I am not able to judge if it is a game changer, but it does advance the state of the art in this domain, although the M/EEG practitioners often ignore new advanced models and only use simple ones. To make this work actually used by the community, an effort should be made to release the code with extensive documentation and examples, ideally in Python (the shared code is in a proprietary language (matlab)).\n\nminor comments:\n- hyperlinks to references do not work\n- L44: M, N and T are not defined\n- Fig1: really low resolution", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Some discussion is provided about the limitations of the model assumptions.\n\nReporting computational cost of each method in the experiments would be appreciated.\n\nSocietal impact is not discussed, although it is probably not necessary.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "8uWPiuIOrjN", "reviewer": "Reviewer_bAFS", "summary": "This manuscript presents an efficient majorization-minimization optimization framework to infer the parameters of the spatio-temporal covariance matrix in the MEEG signal source reconstruction problem. The main contribution lies in simplifying a non-convex optimization problem using a convex surrogate function with a convergence guarantee. The method has been tested in few limited scenarios on simulated and real MEEG data and the results are compared against three alternative approaches. \n", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "*** Main Comments ***\n- This paper is addressing an interesting and significant problem in the neuroimaging context.\n- The idea of using the specifications of the Kronecker product in disentangling the spatio-temporal structure of the covariance matrix is interesting, however, this idea is not original and was presented in previous studies for the purpose other than MEEG source reconstruction.\n- The text is clear but dense in mathematical notations and representations.\n- This work's significance is limited by the strong (restrictive and wrong) assumptions that are made on 1) the independence of the sources across different brain locations, 2) identical temporal distribution for signal and noise. The authors are aware of these practical limitations and conducted a simulated study in a limited scenario arguing their negligible effect on the outcoming results. \n\n*** Other Comments ***\n- The authors categorize the problem in Eq. 1 as a multi-task regression problem, while originally the addressed problem diverges from the definition of multi-task learning in which the aim is to estimate $\\mathbf{L}$ and not $\\mathbf{X}$ (unlike brain source reconstruction).\n- I would call $\\Gamma$ as spatial and $\\mathbf{B}$ as temporal covariance matrix because even though the $\\Gamma$ is computed on time-series it at the end represents the spatial structure of the signal across different brain locations (and the same for $\\mathbf{B}$).\n- According to the ground-truth covariance matrix in Figure 2, the simulated signal can barely be considered structured. I suggest simulating data with a more salient structure in the signal and noise components to better evaluate the violation of assumptions.\n- Despite considering the spatial structure of data, I wonder why the reconstructed sources are so sparse in Figure 3. I would expect a more smooth localized effect. \n- It is not clear to me in what sense the proposed model is hierarchical. In the Bayesian context, hierarchical models have a specific definition (using the same prior in a hierarchy over a group of parameters). ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The proposed method in this study suffers from strong unrealistic assumptions on the distribution of noise and signal. The authors are aware of these limitations and conducted a simulation study to investigate their effect. However, the conducted simulated study in not well-designed with respect to actual realistic scenarios. I suggest:\n\n- to simulate data with a more salient spatio-temporal structure (the covariance matrix of the simulated data is almost diagonal in Figure 2) in both signal and noise components (of course different structures for signal and noise).", "review_text": "", "extra_scores": {"ethical_concerns": "NA", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "1tKUS8cTRI", "reviewer": "Reviewer_txKn", "summary": "This paper proposes an efficient method for source reconstruction from M/EEG data that estimates the spatiotemporal correlation structure in both the source space and the measurement error. The spatiotemporal pattern of the source and noise covariances is given by a Kronecker product structure with an independent (diagonal) spatial structure and either a full or Toeplitz temporal structure, where the temporal component is shared across the source and noise components. Since the full optimization problem is non-convex, the authors propose a majorization-minimization (MM) algorithm with provable convergence and derive simple, efficient iterative updates for each component. Two versions of the model are proposed: full Dugh, which uses a full positive definite matrix for the temporal covariance; and thin Dugh, which can be applied in the stationary data setting and uses a Toeplitz temporal covariance structure with an alternative set of efficient updates for the MM algorithm. On simulated data with known sources, thin Dugh achieves comparable or better performance as compared to existing source reconstruction methods (Champagne, eLORETA, and S-FLEX). When applied to simulated data that was generated while varying the true underlying temporal covariance structure, full Dugh is shown to outperform thin Dugh when the true matrix is full and vice versa when it has Toeplitz structure. Finally, experiments on real MEG data from subjects presented with visual or auditory stimuli show that full and thin Dugh can recover plausible and well-localized temporal and spatial patterns for the reconstructed sources, even with a very limited number of trials.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Overall, the paper is clearly written. The presentation of the method is sound, there is enough technical detail to understand the model and algorithm, and the experiments do a sufficient job of demonstrating the Dugh method's utility both when compared to other source estimation approaches on simulated data and for recovering plausible localized activity on real MEG data. I particularly appreciated the comparison between thin Dugh and full Dugh, both on the real MEG data as well as on synthetic data with known structure to demonstrate that each version of the method excels when the underlying structure aligns with exactly what it was designed for.\n\nOne issue I had with the paper was the limited discussion of and comparison to other source imaging methods that account for temporal structure in some way. This could include approaches such as mixed-norm inverse methods, which do encourage some temporal smoothness in the estimated source activations even if they don't model a full temporal covariance structure as the proposed method does. While such other methods are briefly referred to in the paper, comparing to other approaches that model temporal structure could further showcase the particular strength of the Kronecker-factored covariance structure and associated optimization algorithm that the authors present.\n\nAnother interesting experiment would have been a runtime comparison between the proposed method and the other methods used in the experiments. The authors claim their algorithm is efficient and derive a set of simple closed-form updates, but no direct comparison of efficiency is made against the other source estimation methods. In particular, this would validate the point the authors make about other methods that model temporal structure having \"computationally over-demanding inference schemes\" (line 6). A runtime vs. accuracy comparison to other spatio-temporal methods as well as the baselines already included in the paper would strengthen the presentation of the proposed Dugh method as an accurate source estimation method that is also less computationally intensive compared to alternatives that also model temporal structure.\n\nOther minor comments:\n- The caption for Figure 2 mentions \"Toeplitz structure in the first row, and full covariance structure in the second row\", when the figure shows the reverse", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors have discussed the technical limitations of their work, but not the potential negative societal impact. The most salient potential impact of this work is its use in neuroscience research, where any potential technical biases or shortcomings of the method could be reflected in the findings of any researcher who applies this method to their data. The authors could briefly discuss this potential impact in the final section of the paper. Beyond this, there are no obvious harmful societal consequences of this work to my knowledge.", "review_text": "", "extra_scores": {"ethical_concerns": "No ethical issues come to mind.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The authors present an efficient method to infer models with both spatial and temporal correlations by leveraging Kronoker, sparse, and toeplitz structures. The authors found the method interesting and highly practical as a way to bring additional structure into the regression framework without undue computational burden. While the reviewers did raise a few minor points related to the restrictive model assumptions and the potential for additional comparisons, the authors have provided thorough responses that the reviewers have found adequate. Thus I am happy to recommend this work for acceptance to NeurIPS.
Decision: Accept (Poster)