Efficient Hierarchical Bayesian Inference for Spatio-temporal Regression Models in Neuroimaging
Year: 2022

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.

heading: 1 Introduction
text: Probabilistic graphical models for regression problems, where both the model parameters and the observation noise can have complex spatio-temporal correlation structure, are prevalent in neuroscience, neuroimaging and beyond [1]. Example domains where these models are being applied include Gaussian process inference [2], sparse signal recovery and multi-task learning [3–6], array signal processing [7], fMRI analysis [8], climate science [9], computer hardware monitoring [10], and brain-computer interfaces (e.g. [11, 12]). 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- or electroetoencephalography (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–17] or have ignored temporal correlations overall [5], and therefore not fully addressed the inherent spatio-temporal problem structure.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
We propose a novel and efficient algorithmic framework using Type-II Bayesian regression that explicitly 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, we are then able to derive robust, fast and efficient majorization-minimization (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.

heading: 2 Spatio-temporal generative model
text: Let us consider a multi-task linear regression problem, mathematically formulated as
Yg = LXg + Eg for g = 1, . . . , G , (1)
in which a forward matrix, L ∈ RM×N , maps a set of coefficients or source components Xg ∈ RN×T to measurements Yg ∈ RM×T , with independent Gaussian noise Eg ∈ RM×T . G denotes the number of sample blocks in the multi-task problem, whileN , 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}Gg=1 given L and {Yg}Gg=1 can represent an inverse problem in physics, a multi-task regression problem in machine learning, or a multiple measurement vector (MMV) recovery problem in signal processing [18].
In the context of BSI, {Yg}Gg=1 refers to M/EEG sensor measurement data, and {Xg}Gg=1 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 as the finite element method for a given head geometry and known electrical conductivities using the quasi-static approximation of Maxwell’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 ill-posed, which can be dealt with by incorporating prior assumptions. Adopting a 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-Posteriori (MAP) estimation (Type-I Bayesian learning) [20–22] or, when the model has unknown hyperparameters, through Type-II MaximumLikelihood (ML) estimation (Type-II Bayesian learning) [23–25].
Here we focus on full spatio-temporal Type-II Bayesian learning, which assumes a family of prior distributions p(Xg|Θ) parameterized by a set of hyperparameters Θ. We further assume that the spatio-temporal correlation structure of the brain sources can be modeled by a Gaussian p(xg|Γ,B) ∼ N (0,Σ0), where xg = vec(X>g ) ∈ RNT×1 and where the covariance matrix Σ0 = Γ ⊗ B has Kronecker structure, meaning that the temporal and spatial correlations of the sources are modeled independently through matrices Γ and B, where ⊗ stands for the Kronecker product. Note that the decoupling of spatial and temporal covariance is a neurophysiologically 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 Γ = diag(γ), γ = [γ1, . . . , γN ]>, and B ∈ ST++, where ST++ denotes the set of positive definite matrices with size T × T . Note that using a diagonal matrix Γ amounts to assuming that sources in different locations are independent. This assumption may be relaxed in our future work. On the other hand, modeling a full temporal covariance matrix B amounts to assuming non-stationary 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((Xg)i.|γi,B) ∼ N (0, γiB), for i = 1, . . . , N , where (Xg)i. denotes the i-th row of source matrix Xg .
Analogous to the definition of the sources, we here also model the noise to have spatio-temporal structure. The matrix Eg ∈ RM×T represents T time instances of zero-mean Gaussian noise with full covariance, Σe = Λ⊗Υ, where eg = vec(E>g ) ∈ RMT×1 ∼ N (0,Σe), and where Λ and Υ 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 Λ = diag(λ), λ = [λ1, . . . , λM ]>. We consider both the general, heteroscedastic, case in which the noise level may be different for each sensor as well as the homoscedastic case in which the noise level is 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 := [Γ,0; 0,Λ], and Φ := [L, I]. These 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 ηg := [x>g , e > g ] to the measurements; i.e., yg = Φηg. Figure 1-(a) illustrates a probabilistic graphical model summarizing the spatio-temporal generative model of our multi-task linear regression problem.
The MMV model Eq. (1) can be formulated equivalently in single measurement vector (SMV) form by vectorizing the spatio-temporal data matrices and using the Kronecker product as follows: yg = Dxg+eg , where yg = vec ( Y>g ) ∈ RMT×1, and D = L⊗IT . In a Type-II Bayesian learning framework, the hyper-parameters of the spatio-temporal source model are optimized jointly with the model parameters {Xg}Gg=1. In our case, these hyper-parameters comprise the unknown source, noise and temporal covariance matrices, i.e., Θ = {Γ,Λ,B}.

heading: 3 Proposed method — full Dugh
text: The unknown parameters Γ, Λ, and B can be optimized in an alternating iterative process. Given initial estimates, the posterior distribution of the sources is a Gaussian of the form p(xg|yg,Γ,Λ,B) ∼ N (x̄g,Σx), whose mean and covariance are defined as follows:
x̄g = Σ0D >(Λ⊗B + DΣ0D>)−1yg = Σ0D>Σ̃−1y yg , (2)
Σx = Σ0 −Σ0D>Σ̃−1y DΣ0 , (3)
where Σy = LΓL> + Λ, and where Σ̃y = Σy ⊗B denotes the spatio-temporal statistical model covariance matrix. The estimated posterior parameters x̄g and Σx are then in turn used to update Γ, Λ, and B as the minimizers of the negative log of the marginal likelihood p(Y|Γ,Λ,B), which is given by:
Lkron(Γ,Λ,B) = T log |Σy|+M log |B|+ 1
G G∑ g=1 tr(Σ−1y YgB −1Y>g ) , (4)
where | · | denotes the determinant of a matrix. A detailed derivation is provided in Appendix A. Given the solution of the hyperparameters Γ, Λ, and B, the posterior source distribution is obtained by plugging these estimates into Eqs. (2)–(3). This process is repeated until convergence.
The challenge in this high-dimensional 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 majorizationminimization (MM) [27] is a powerful non-linear optimization framework that can be leveraged to solve similar Bayesian Type-II inference problems. Here we extend this work to our spatiotemporal setting. Building on the idea of majorization-minimization, we construct convex surrogate functions that majorize Lkron(Γ,Λ,B) in each iteration of the proposed algorithm. Then, we show the minimization equivalence between the constructed majorizing functions and Eq. (4). These results are presented in theorems 1 and 3. Theorems 2 and 4 propose an efficient alternating optimization algorithm for solving Lkron(Γ,Λ,B), which leads to update rules for the spatial and temporal covariance matrices Γ and B as well as the noise covariance matrix Λ.
Starting with the estimation of the temporal covariance based on the current source estimate, we can state the following theorem:
Theorem 1. Optimizing the non-convex Type-II cost function in Eq. (4), Lkron(Γ,Λ,B), with respect to the temporal covariance matrix B is equivalent to optimizing the following convex surrogate function, which majorizes Eq. (4):
Ltimeconv(Γk,Λk,B) = tr ( (Bk)−1B ) + tr(MktimeB −1), (5)
where Mktime is defined as
Mktime := 1
MG G∑ g=1 Y>g ( Σky )−1 Yg , (6)
and where Γk, Λk, and Σky 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 Ltimeconv(Γk,Λk,B) is strictly geodesically convex with respect to the P.D. manifold and its minimum with respect to B can be attained by iterating the following update rule until convergence:
Bk+1 ← (Bk)1/2 (
(Bk) −1/2Mktime(B
k) −1/2 )1/2 (Bk) 1/2 . (7)
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, Mktime, and the model temporal covariance matrix from the previous iteration, Bk. A geometric representation of the geodesic path between the pair of matrices {Bk,Mktime} on the P.D. manifold and the geometric mean between them, representing the update for Bk+1, is provided in Figure 1-(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 a convex surrogate function that majorizes Lkron(Γ,Λ,B) in each iteration of the optimization algorithm by considering augmented variables comprising both the source and noise covariances, i.e., H := [Γ,0; 0,Λ], and Φ := [L, I]. The following theorem holds. Theorem 3. Minimizing the non-convex Type-II cost function in Eq. (4), Lkron(Γ,Λ,B), with respect to the spatial covariance matrix H is equivalent to minimizing the following convex surrogate function, which majorizes Eq. (4):
Lspaceconv (Γ,Λ,Bk) = Lspaceconv (H,Bk) = tr ( Φ>(Σky) −1ΦH ) + tr ( MkSNH −1) , (8)
where MkSN is defined as
MkSN := H kΦ>(Σky) −1Mkspace(Σ k y) −1ΦHk ,with
Mkspace := 1
TG G∑ g=1 Yg(B k)−1Y>g , (9)
and where Bk 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), h = [γ1, . . . , γN , σ21 , . . . , σ 2 M ] >, and Σy = ΦHΦ>. The update rule for H then takes a simple form, as stated in the following theorem:
Theorem 4. The cost function Lspaceconv (H,Bk) is convex in h, and its minimum with respect to h can be obtained according to the following closed-form update rule, which concurrently estimates the scalar source variances and heteroscedastic noise variances:
Hk+1 = diag(hk+1), hk+1i ← √ gki zki for i = 1, . . . , N +M ,where (10)
g := diag(MkSN) (11) z := diag(Φ>(Σky) −1Φ) (12)
Proof. A detailed proof can be found in Appendix E.
The final estimate of the posterior distribution (2) can be obtained by starting from any initial guess for the model parameters, H0 = {Γ0,Λ0} and B0, and iterating between update rules (7), (10), and (2) 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. (4) with alternating update rules for B and H = {Γ,Λ} in Eq. (7) and Eq. (10) defines an MM algorithm, which is guaranteed to converge to a stationary point.
Proof. A detailed proof can be found in Appendix G.

heading: 4 Efficient inference for stationary temporal dynamics — thin Dugh
text: With full Dugh, we could present a 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 on 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 (eventrelated 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) to 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 xn(t+ 1) = βxn(t) + √ 1− β2ξn(t), n = 1, . . . , N ; t = 1, . . . , T with AR coefficient β ∈ (−1, 1) and innovation noise ξn(t). It can be shown that the temporal correlation matrix corresponding to this AR
model has Toeplitz structure, Bi,j = β|i−j| [30]. Consequently, we now constrain the cost function Eq. (5) to the set of Toeplitz matrices:
Bk+1 = arg min B∈B, H=Hk tr((Bk)−1B) + tr(MktimeB −1) , (13)
where set B denotes the set of real-valued positive-definite Toeplitz matrices of size T × T . In order to be able to derive efficient update rules for this setting, we bring B into a diagonal form. The following proposition holds.
Proposition 1. Let P ∈ RL×L with L > T , be the circulant embedding of matrix B ∈ 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 B can be diagonalized as follows:
B = QPQH with Q = [IM ,0]FL, and (14)
[FL]m,l = 1√ L e i2π(l−1) L (m−1) , (15)
where P = diag(p) = diag(p0, p1, . . . , pL−1) is a diagonal matrix. The main diagonal coefficients of P are given as the normalized discrete Fourier transform (DFT), represented by the linear operator FL, of the first row of B: p = FLB1. (Note that a Toeplitz matrix can be represented by its first row or column, B1.).
Proof. This is a direct implication of the fact that Toeplitz matrices can be embedded into circulant matrices of larger size [31] and the result that circulant matrices can be approximately diagonalized using the Fourier transform [32]. Further details can be found in Appendix H.
The solution of Eq. (13) can be obtained by direct application of Proposition 1. 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 p can be obtained by iterating the following closed-form update rule until convergence:
pk+1l ← √ ĝkl ẑkl for l = 1, . . . , L ,where (16)
ĝ := diag(PkQH(Bk)−1Mktime(B k)−1QPk) (17) ẑ := diag(QH(Bk)−1Q) . (18)
Proof. A detailed proof is provided in Appendix I.

heading: 4.1 Efficient computation of the posterior
text: Given estimates of the spatial and temporal covariance matrices, we can efficiently compute the posterior mean by exploiting their intrinsic diagonal structure.
Theorem 7. Given the diagonalization B = QPQH of the temporal correlation matrix and the eigenvalue decomposition LΓL> = UxDxU>x of LΓL
>, where Dx = diag(d1, . . . , dM ), the posterior mean is efficiently computed as
x̄g = (Γ⊗B) D>Σ̃−1y yg = tr ( QP ( Π QHY>g Ux ) ( U>x LΓ >)) , (19) where denotes the Hadamard product between the corresponding elements of two matrices of equal size. In addition, Π is defined as follows: [Π]l,m = 1σ2m+pldm for l = 1, . . . , L and m = 1, . . . ,M .
Proof. A detailed proof is provided in Appendix J.
The resulting algorithm, obtained by iterating between update rules (16), (10) and (19), is called thin Dugh (as opposed to full Dugh introduced above).

heading: 5 Simulations
text: 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 [33] and two other competitive methods – eLORETA [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.653s, and S-FLEX: 20.963s. The codes are publicly available at https://github.com/AliHashemi-ai/Dugh-NeurIPS-2021.

heading: 5.1 Pseudo-EEG signal generation and benchmark comparison
text: We simulate a sparse set of N0 = 3 active sources placed at random locations on the cortex. To simulate the electrical neural activity, we sample time series of length T ∈ {10, 20, 50, 100} from a univariate linear autoregressive AR(P) process. We use stable AR systems of order P ∈ {1, 2, 5, 7}. The resulting source distribution is then projected to the EEG sensors, denoted by Ysignal, using a realistic lead field matrix, L ∈ R58×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 a 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 Ysignal as follows: Y = Ysignal + (1−α)‖Ysignal‖ F
α‖E‖F E, where α determines the signal-to-noise ratio (SNR) in
sensor space. Precisely, SNR is defined as SNR = 20log10 (α/1−α). In this experiment the following values of α are used: α ∈ {0.55, 0.65, 0.7, 0.8}, which correspond to the following SNRs: SNR ∈ {1.7, 5.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 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 to investigate the impact of model violations on the temporal covariance estimation, we generate a random Gaussian source matrix, X ∈ R2004×30×G representing the brain activity of 2004 brain sources at 30 time instances for different numbers of trials G ∈ 10, 20, 30, 40, 50. In all trials, sources are randomly sampled from a zero-mean normal distribution with spatio-temporal covariance matrix Γ ⊗ B, where B ∈ R30×30 is either a full-structural random PSD matrix or a Toeplitz matrix with β = 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 0 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̂ − B||2F /||B||2F . 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 scale-invariant and hence only quantifies the overall structural similarity between simulated and estimated noise covariance matrices.
In Figure 2, 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 is 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. Figure 3 demonstrates the estimated temporal covariance matrix obtained 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 a
setting where the ground-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 to thin Dugh that is only able to recover a 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 a full-structure temporal covariance matrix. Finally, the performance of both learning schemes are significantly improved by increasing the number of trials.

heading: 6 Analysis of real MEG and EEG data recordings
text: We now demonstrate the performance of the novel algorithms on two MEG datasets (see Figure 4 and

heading: 5) and one EEG dataset (see Appendix L). 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
text: 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. No 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 a 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 to 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 to 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 at 100 ms and a later peak around 200 ms (similar performance was observed for full Dugh).

heading: 7 Discussion
text: Inverse modeling is a challenging problem that can be cast as a multi-task regression problem. Machine learning-based methods have contributed by systematically incorporating prior information into such models (e.g. [13, 25, 39–41]). While, so far, explicit models of spatio-temporal dynamics 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, we have demonstrated that the inclusion of both spatial and temporal model parameters in the model indeed leads to significantly improved reconstruction performance. Finally, the utility of our algorithms is showcased in challenging real-world applications by reconstructing expected sources from real M/EEG data based on very few experimental trials.
The M/EEG BSI problem has been previously framed as a 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-task 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 a Type-II 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, as in [43]. A number of alternative approaches have been proposed to estimate the spatio-temporal correlation structure of the noise [44–48]. These works, however, do not estimate the noise characteristics 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-end learning philosophy as opposed to a step-wise independent estimation process that can give rise to error 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 to achieve tractable inference due to the fact that the spatio-temporal statistical model covariance matrix can be formulated in a separable form, i.e., Σ̃y = Σy ⊗ 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 a 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 regression problems with spatio-temporal covariances. Incorporating prior knowledge, namely, constraining the solutions to Riemannian manifolds, and using ideas of geodesic convexity, circulant embeddings, and majorization-minimization, 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.
Acknowledgments and Disclosure of Funding
This result is part of a 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/Intelligent Data Analysis research group at Technische Universität Berlin. He further wishes to thank the Charité – Universitätsmedizin 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 Education and Research as BIFOLD – Berlin Institute for the Foundations of Learning and Data (ref. 01IS18025A and ref. 01IS18037A), and the German Research Foundation (DFG) as Math+: Berlin Mathematics Research Center (EXC 2046/1, project-ID: 390685689), Institute of Information & Communications Technology Planning & Evaluation (IITP) grants funded by the Korea Government (No. 2019-0-00079, Artificial Intelligence Graduate School Program, Korea University). SSN was funded in part by National Institutes of Health grants (R01DC004855, R01EB022717, R01DC176960, R01DC010145, R01NS100440, R01AG062196, and R01DC013979), University of California MRPI MRP-17–454755, the US Department of Defense grant (W81XWH-13-1-0494).

Review summaries:
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  Kronecker 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.

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.

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.  

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.