Joint Learning of Full-structure Noise in Hierarchical Bayesian Regression Models
Authors: Ali Hashemi
Chang Cai
Klaus Robert Muller
Srikantan Nagarajan
Stefan Haufe
Venue:  (2021)
Keywords: Full-structure Noise
Hierarchical Bayesian Regression Models
Sparse Bayesian Learning
Unsupervised Learning
Brain Source Imaging
Covariance Estimation.

Abstract: We consider hierarchical Bayesian (type-II maximum likelihood) models for observations with latent variables for source and noise, where both hyperparameters need to be estimated jointly from data. This problem has application in many domains in imaging including biomagnetic inverse problems. Crucial factors influencing accuracy of source estimation are not only the noise level but also its correlation structure, but existing approaches have not addressed estimation of noise covariance matrices with full structure. Here, we consider the reconstruction of brain activity from electroencephalography (EEG). This inverse problem can be formulated as a linear regression with independent Gaussian scale mixture priors for both the source and noise components. As a departure from  classical sparse Bayesan learning (SBL) models where across-sensor observations are assumed to be independent and identically distributed, we consider Gaussian noise with full covariance structure. Using  Riemannian geometry, we derive an efficient algorithm for updating both source and noise covariance along the manifold of positive definite matrices. Using the majorization-maximization framework, we demonstrate that our algorithm has guaranteed and fast convergence. We validate the algorithm both in simulations and with real data. Our results demonstrate that the novel framework significantly improves upon state-of-the-art techniques in the real-world scenario where the noise is indeed non-diagonal and fully-structured.

Under review as a conference paper at ICLR 2021
JOINT LEARNING OF FULL-STRUCTURE NOISE IN HIERARCHICAL BAYESIAN REGRESSION MODELS
Anonymous authors Paper under double-blind review
ABSTRACT
We consider hierarchical Bayesian (type-ll maximum likelihood) models for ob- servations with latent variables for source and noise, where both hyperparameters need to be estimated jointly from data. This problem has application in many domains in imaging including biomagnetic inverse problems. Crucial factors in- fluencing accuracy of source estimation are not only the noise level but also its correlation structure, but existing approaches have not addressed estimation noise covariance matrices with full structure. Here, consider the reconstruc- tion brain activity from electroenceepaalggrahyy (EEG). This inverse problem linear regression with independent Gaussian scale mixture
can be formulated priors for both the source and noise components. As departure from classical sparse Bayesan learning (SBL) models where across-sensor observations are as- sumed to be independent and identically distributed, we consider Gaussian noise with full covariance structure. Using Riemannian geometry. we derive an efficient algorithm for updating both source and noise covariance along the manifold of positive definite matrices. Using the demonstrate that our algorithm has guaranteed and fast convergence. We validate the algorithm both in simulations and with real data. Our results demonstrate that the novel framework significantly improves upon state-of-the-arr techniques in the real-world scenario where the noise is indeed non-diagonal and fully-structured.
crrrization-mxximizztte
framework, we
INTRODUCTION
Having precise knowledge of the noise distribution is a fundamental requirement for obtaining aceu- rate solutions in many regression problems (Bungert et al. 2020). In many applications however, it is impossible to separately estimate this noise distribution, as distinct "noise-only" (baseline) measure- ments are not feasible. An alternative, therefore, is to design estimators that jointly optimize over the regression coefficients as well as over parameters of the noise distribution. This has been pursued both in a (penalized) maximum wawnnnu--ikklhooooo likelihood settings (here referred to as Type-I approaches) (Petersen & Jung 2020 [Bertrand et al. 2019: Massias et al. (2018) as well as in hierarchical Bayesian settings (referred to Type-If) (Wipf & Rao [007]] Zhang & Rao 2011 Hashemi al., 2020 Cai al.,
2020a). Most contributions in the literature are, however, limited to the estimation of only diago- nal noise covariance (i.c., independent between different measurements) (Daye al. 2012 Van de
Geer et al.] 2013 Dalalyan et al., 2013 Lederer & Muller 2015). Considering diagonal noise covariance a limiting assumption practice as the noise interference in many realistic scenarios are highly correlated across measurements; and thus, have non-trivial off-diagonal elements.
This paper develops an efficient optimization algorithm for jointly estimating the posterior of regres- sion parameters as well as the noise distribution. More specifically, we consider linear regression with Gaussian scale mixture priors on the parameters and a full-structure multivariate Gaussian noise. We cast the problem problem, in which the variance hyperparameters and the noise covariance matrix are optimized by maximizing the Bayesian evidence of the model. Using Riemannian geometry, derive effi- cient algorithm for jointly estimating the source and noise covariances along the manifold of positive definite (P.D.) matrices.
hierarchical Bayesian (type-Il maximum-likelihoud likelihood) regression
To highlight the benefits of our proposed method in practical scenarios, we consider the problem of electromagnetic brain source imaging (BSI). The goal of BSI is to reconstruct brain activity

---

Under review as a conference paper at ICLR 2021
from magneto- or eleetroenceephallgraphy (M/EEG), which can be formulated sparse Bayesian learning (SBL) problem. Specifically, it can be cast as a linear Bayesian regression model with independent Gaussian scale mixture priors on the parameters and noise. As a departure from the classical SBI. approaches, here we specifically consider Gaussian noise with full covariance struc- ture. Prominent source of correlated noise in this context are, for example, eye blinks, heart beats, muscular artifacts and line noise. Other realistic examples for the need for such full-structure noise can be found in the areas of array processing (Li & Nehorai 2010) direction of arrival (DOA) estimation (Chen et al., 2008). Algorithms that can accurately estimate noise with full covariance structure are expected to achieve more accurate regression models and predictions in this setting.
2 TYPE-II BAYESIAN REGRESSION
We consider the linear model Y LX + E, in which a forward or design matrix, L E MxN, is mapped to the measurements, Y. by a set of coefficients or source components, X. Depending on the setting, the problem of estimating X given L and Y called an inverse problem in physics, a multi- task regression problem in machine learning. or a multiple measurement vector (MMV) recovery problem signal processing (Cotter et al. 2005). Adopting a signal processing terminology, the measurement matrix Y E RI- captures the activity of M sensors at T time instants, (() E
RMX1, = 1,..., T. while the source matrix, X R N xT" consists
the unknown activity of N
sources at the same time instants, x(t) E N 1...,. T. The matrix E = (e(1)...... ..... e(T)]
................................ IR RMX represents T time instances of zero-mean Gaussian noise with full covariance A, e(t) e
R A), 1,...: T, which is assumed to be independent of the source activations.
In this paper, we focus on M/EEG based brain source imaging (BSI) but the proposed algorithm can be used in general regression settings, in particular for sparse signal recovery (Candès et al. 2006 Donoho 2006) with a wide range of applications (Malioutov et al., 2005). The goal of BSI is to infer the underlying brain activity X from the BEG/MEG measurement Y given a known forward operator, called lead field matrix L. As the number of sensors is typically much smaller than the number of locations of potential brain sources, this inverse problem is highly ill-posed. This prob- lem is addressed by imposing prior distributions on the model parameters and adopting a Bayesian treatment. This can be performed either through Maximum-a -aattmumaaPooteriooii (MAP) estimation (Type-I Bayesian learning) (Pascual- Marqui et al. 1994, Gorodnitsky et al.) 1995 Haufe et al.) 2008 Gram- fort et al. 2012 Castanio-Canammi et al. 2015) or, when the model has unknown hyperparameters.
through Type-I1 Maximum-Likelihen estimation (Type-1l Bayesian learning) (Mika et al. 2000
Tipping 2001: Wipf & Nagarajan 2009
Seeger & Wipi 2010 Wu et al. (016))
In this paper, we focus on Type-II Bayesian learning. which assumes a family of prior distributions ((x(e) parameterized by set of hyperparameters @ These hyper-parameteen can be learned from the data along with the model parameters using hierarchical Bayesian approach (Tipping 2001
Wipf & Rao 2004) through the maximum-likelihod principle:
/ X e
e| arg max pl O = arg max
Here we assume zero-mean Gaussian prior with full covariance I for the underlying source dis- tribution, x(t) E R Nx1 ~ (), to I....,T. T. Just as most other approaches, Type-Il Bayesian learning makes the simplifying assumption of statistical independence between time samples. This leads to the following expression for the distribution of the sources and measurements:
T
T II NO.. r)
PX II
p (()
(2)
1=1
1=1
P Y X II I(x(())())) = II ((L(())
(3)
The parameters of the Type-Il model. + the unknown source and noise covariances, i.e.. θ A). The unknown parameters I and A are optimized based on the current estimates of the source and noise covariances in alternating iterative process. Given initial estimates of I and A,
2

---

Under review as a conference paper at ICLR 2021
the posterior distribution of the sources is Gaussian of the form (Sekihara & Nagarajan
2015
P r II N (1x (t),
where
(4)
[=1
Mx TL EEy) yy())
(5)
Σx = PL 2
(6)
My A+ LIL
(7)
The estimated posterior parameters (t) and Σx are then in turn used to update I and as the minimizers of the negative log of the marginal likelihood p Y YT, A), which is given by (Wipf al. 2010):
T Σ y(t) yy-y(t)
C (T,A)) = - log p( (Y|T, A) = log Eyl +
T
1=1
T
= log/ + LIL
Σ y(t) (A LF y())) t=1
(8)
T
where denotes the determinant of matrix. This process is repeated until convergence. Given the final solution of the hyperparameters @ (T"I. "), the posterior source distribution obtained by plugging these estimates into equations to
3 PROPOSED METHOD: FULL-STRUCTURE NOISE (FUN) LEARNING
Here we propose a novel and efficient algorithm, full-structure noise (FUN) learning, which is able to learn the full covariance structure of the noise jointly within the Bayesian Type-II regression framework. We first formulate the algorithm in its most general form, in which both the noise distribution and the prior have full covariance structure. Later, we make the simplifying assumption of independent source priors, leading to the pruning of the majority of sources. This effect, which has also been referred (SBL) is beneficial in our application of interest, namely the reconstruction of parsimonious sets of brain sources underlying experimental EEG measurements.
automatic relevance determination (ARD) or sparse Bayesian learning
Note that the Type-II cost function in equation is non-convex and thus non-trivial to optimize. number of iterative algorithms such as sngiorrzation-miinniaazn been proposed to address this challenge. Following the MM scheme, we first construct convex surrogate functions that majorizes [T A) in each iteration of the optimization algorithm. Then, we show the minimization equivalence between the constructed majoring functions and equation This result is presented in the following theorem:
(MM) (Sun et al. 2017) have
Theorem 1. Let A and Σ/ be fixed values obtained in the (k)--h iteration of the optimization
algorithm minimizing L" (T, A). Then, optimizing the non-convex type- ML cost function in equa- tion 8 (T, A), with respect to is equivalent to optimizing the following convex function, which majorizes equation
Leouy (T, (T,A²) (r((C}) -'r) (r(M
where C and M are defined as:
c = L 1 1)
M = T xt)xt))
(10)
Similarly, optimizing C (T, A) with respect to A is equivalent to optimizing the following comvex
majorizing function:
Coony @ ( tr((Ck) A) + ++
(11)
where C and My are defined as:
T M = T -x00) -xx)))))))x(0)))))))))))) L
: = (Σ*)
(12)
3

---

Under review as a conference paper at ICLR 2021
Proof. The proof is presented in Appendix A
We continue by considering the optimization of the cost functions Lourrc A+) and Loow (T* A) with respect to r and A, respectively. Note that in case of source covariances with full structure, the solution of Lanuree (T, A with respect to I lies in the (N2 N2 Riemannian manifold of positive definite (P.D.) matrices. This consideration enables us to invoke efficient methods from Riemannian geometry (see [Petersen al., 2006: Berger 2012 Jost & Jost, 2008), which ensures that the solution at each step of the optimization is contained within the lower-dimensional solu- tion space. Specifically, in order optimize for the source covariance, the algorithm calculates the geometric mean between the previously obtained statistical model source covariance, C and the source-space sample covariance matrix, M in each iteration. Analogously, to update the noise co- variance estimate, the algorithm calculates the geometric mean between the model noise covariance, C and the empirical censor-space residuals, M The update rules obtained from this algorithm are presented in the following theorem:
Theorem 2. The cost functions Loource (T, Ax) and poise (T&y R, A) are both strictly geodesically convex with respect to the P.D. manifold, and their optimal solution with respect to I and A. respec- tively, can be attained according the two following update rules:
p < (C$) ((C)) (C5) vV) (Cg)
(13)
A A+1 - (ck) )
(C* v+) (C*)
(14)
Proof. A detailed proof can be found in Appendix
Convergence of the resulting algorithm is shown in the following theorem.
Theorem 3. Optimizing the non-convex type-ll ML cost function in equation 8 C (T, A) with
alternating update rules for I and A in equation 13 and equation 14 leads with guaranteed convergence guarantees.
MM algorithm
Proof. A detailed proof can be found in Appendix C
While Theorems 3 reflect a general joint learning algorithm, the assumption of sources with full covariance structure often relaxed in practice. The next section will shed light on this important simplification by making a formal connection to SBL algorithms.
3.1 SPARSE BAYESIAN LEARNING WITH FULL NOISE MODELING
In brain source imaging, the assumption of full source covariance is often relaxed. Even if, tech- nically, most parts of the brain are active at all times, and the concurrent activations of different brain regions can never be assumed to be fully uncorrelated, there are many experimental settings in which it is reasonable to assume only small set of independent brain sources. Such sparse solutions are physiologically plausible in task-based analyses, where only a fraction the brain's 's macroscopic structures is expected to be consistently engaged. A common strategy in this case is to model independent sources through a diagonal covariance matrix. In the Type-ll Bayesian learning framework, this simplification interestingly leads to sparsity of the resulting source distributions, as, at the optimum, many of the estimated source variances are zero. This mechanism is known as sparse Bayesian learning and is closely related to the more general concept of automatic relevance determination. Here, we adopt the SBL assumption for the sources, leading to I -uppdates previously described in the BSI literature under the name Champagne (Wipf & Nagarajan 2009). As novelty and main focus of this paper, we here equip the SBL framework with the capability to jointly learn full noise covariances through the geometric mean based update rule in equation [4] In the SBL framework, the N modeled brain sources are assumed to follow independent univariate Gaussian distributions with zero mean and distinct unknown variances Tn (t) ~ (N 7n), n I....,N. In the SBL solution, the majority of variances is zero, thus effectively inducing spatial sparsity of the corresponding source activities. For FUN learning, we also impose a diagonal structure on the source covariance matrix, I diag(y), where 7 (71,-...70 YN] T. By constraining I' in equation 9
4

---

Under review as a conference paper at ICLR 2021
Algorithm 1: Full-structure noise (FUN) learning
Input: The lead field matrix L E RMXN and the measurement vectors y(t) E RMxx, - T. Result: The estimated prior source variances NN] noise covariance A, the posterior
mean (tx(t) and covariance Σx of the sources.
Set a random initial value for A well y- h....., 7w] and construct r diag(r).
Calculate the statistical covariance Σy A LIL
Repeat
Calculate the posterior mean as μ (1)
LL
(
Calculate C and M based on equation 10 and update Yn for n I,...,N NN based on
equation 15
Calculate C and M based on equation 12
and update A based on equation 14
Until stopping condition is satisfied;
Calculate the posterior covariance as Σx = r IL TL (Ly)
to the set of diagonal matrices, W, we can show that the update rule equation 13 for the source variances simplifies to the following form:
M
D (x), (t))
x+1 In
for I,...,N,
(15)
L E L
n.
where L, denotes the r-th column the lead field matrix. Interestingly, equation 15 is identical to the update rule of the Champagne algorithm. A detailed derivation of equation 15 can be found in Appendix D
Summarizing, the FUN learning approach, just like Champagne and other SBL algorithms, assumes independent Gaussian sources with individual variances (thus, diagonal source covariances), which updated through equation equation [5 Departing from the classical SBL setting, which as- sumes the noise distribution to be known, FUN models noise with full covariance structure, which summarizes the used update rules.
updated using equation 14 Algorithm [
Note that various recent Type-Il noise learning schemes for diagonal noise covariance matrices Hashemi et al., 2020 Cai et al. 2020a) that rooted in the concept of SBL be also derived as special cases of FÜN learning assuming diagonal source and noise covariances, i.e., T, A E W. Specifically imposing diagonal structure on the noise covariance matrix for the FUN algorithm, A, results in identical noise variance update rules as derived in Cai et al. (2020a) for heteroscedastic, and in Hashemii et al. (2020) for homoscedastic noise. We explicitly demonstrate this connection in Appendix E, Here, we note that heteroscedasticiry refers to the common phenomenon that mea- surements are contaminated with non-uniform noise levels across channels, while homoscedasticity only accounts for uniform noise levels.
4 NUMERICAL. SIMULATIONS AND REAL DATA ANALYSIS
Source, Noise and Forward Model: We simulated a sparse set of No = active brain sources that were placed random positions on the cortex. To simulate the electrical neural activity of these sources, T = 200 identically and independently distributed (i.i.d) points were sampled from Gaussian distribution, yielding sparse source activation vectors x(t). The resulting source distribution, represented as X = [x( x(T)). was projected to the EEG sensors through application of lead field matrix as the forward operator: ysignal LX. The lead field ma- trix, L E R58x2004 was generated using the New York Head model (Huang et al.) 2016) tak- ing into account the realistic anatomy and electrical tissue conductivities of an average human head. Further details regarding forward modeling is provided in Appendix F Gaussian additive noise was randomly sampled from a zero-mean normal distribution with full covariance matrix A: e(t) E IR Mx1 (0, A), = 1, .,.... T. This setting is further referred to full-structure noise. Note that we also generated noise with diagonal covariance matrix, referred to as heteroscedas- tic noise, in order to investigate the effect of model violation reconstruction performance. The
5

---

Under review as a conference paper at ICLR 2021
RMx was normalized by it Frobenius norm and added to the signal matrix Y as follows: Y Ysignal + (1-a)]|x"***II =/al| E, where a de- termines the signal-to-noise ratio (SNR) in sensor space. Precisely, SNR is obtained as follows: SNR = 10 (*/1-a). In the subsequently described experiments the following values of were used: (0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.65, 0.7, 0.8}. which correspond to the following SNRs:
noise matrix E
SNR=(-12, -7.4, -5.4, -3.5, -1.7, 0, 1.7, 3.5, 5.4, 7.4, 12) (dB). MATLAB codes for producing the results in the simulation study are uploaded here
We applied the full-structure noise learning approach on the synthetic datasets described above to recover the locations and time courses of the active brain sources. In addition to our proposed approach, two further Type-II Bayesian learning schemes, namely Champagne with homo- and heteroscedastic noise learning (Hashemi et al. 2020 Cai et al. 2020a), were also included as benchmarks with respect to source reconstruction performance and noise covariance estimation accuracy. Source reconstruction performance was evaluated accord- ing to the earth mover's distance (EMD) (Rubner et al. ) the error in the reconstruction of the source time courses, the average Euclidean distance (EUCL) (in mm) between each simulated source and the best (in terms of absolute correlations) matching reconstructed source, and finally Fl-measure score (Chinchor & Sundheim 1993). A detailed definition of evaluation metrics is provided in [] To evaluate the accuracy of the noise covariance matrix estimation, the fol- lowing two metrics were calculated: the Pearson correlation between the original and reconstructed noise covariance matrices, A and A, denoted by A sim, and the normalized mean squared error (NMSE) between A and Å, defined as NMSE || All///|A|1 Note that NMSE measures the reconstruction of the true scale of the noise covariance matrix, while A sim is scale-invariant and hence only quantifies the overall structural similarity between simulated and estimated noise covariance matrices. Each simulation was 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 rep- etitions was calculated. Convergence of the optimization programs for each run was defined if the relative change of the Frobenius-none norm of the reconstructed sources between subsequent iterations was less than 10-8 maximum of 1000 iterations was carried out if no convergence was reached beforehand.
Evaluation Metrics and Simulation Set-up:
Figure [ shows two simulated datasets with five active sources in presence of full-structure noise (upper panel) as well as heteroscedastic noise (lower panel) at 0 (dB) SNR. Topographic maps depict the locations of the ground-truth active brain sources (first column) along with the source reconstruction result of three noise learning schemes assuming noise with homoscedastic (second column), heteroscedastic (third column), and full (fourth column) structure. For each algorithm, the estimated noise covariance matrix is also plotted above the topographic map. Source reconstruction performance was measured in terms of EMD and time course correlation (Corr), and is summarized in the table next each panel. Besides, the accuracy of the noise covariance matrix reconstruction was measured on terms of A sim and NMSE. Results are included in the same table. Figure | (upper panel) allows for a direct comparison of the estimated noise covariance matrices obtained from the three different noise learning schemes. It can be seen that FUN learning can better capture the overall structure of ground truth structurr noise as evidenced by lower NMSE and similarity errors compared to the heteroscedastic and homoscedastic algorithm variants that are only able to recover a diagonal matrix while enforcing the off-diagonal elements to zero. This behaviour results in higher spatial and temporal accuracy (lower EMD and time course error) for FUN learning compared to competing algorithms assuming diagonal noise covariance. This advantage is also visible in the topographic maps. The lower-panel of Figure presents analogous results for the setting where the noise covariance is generated according to a heteroscedasticc model. Note that the superior spatial and temporal reconstruction performance of the heteroscedastic noise learning algorithm compared to the full-structure scheme is expected here because the simulated ground truth noise is indeed heteroscedastic. The full-structure noise learning approach, however, provides fairly reasonable performance in terms of EMD, time course correlation (corT), and Asim, although it is designed to estimate full-structure noise covariance matrix. The convergence behaviour of all three noise learning variants is also illustrated in Figure [] Note that the full-structure noise learning approach eventually reaches lower negative log- likeliholood values in both scenarios, namely full full-structure and heteroscedastic noise.
Figure 2 shows the EMD, the time course reconstruction error, the EUCL and the Fl measure score incurred by three different noise learning approaches assuming homoscedastic (red), heteroscedastic

---

Under review as a conference paper at ICLR 2021
100
Ground Truth
Homoscedastic
Heteroscedastic
Full Structure
Noise Cov
Loss
Full-structural Noise
ocs
12
ca
-100
-200
10
10
Number of iterations
Brain Maps
Ground Truth
Homoscedastic
Heteroscedastic
Full Structure
100
Heteroscedastic Noisa Brain Maps
Noise Cov
6a n.s
20
en
one
12N2S 21
10:0000-300
10
10a
Number of Iterations
Figure I: Two examples of the simulated data with five active sources in presence of full-structure noise (upper panel) as well as heteroscedastic noise (lower panel) at 0 (dB) SNR. Topographic maps depict the locations of the ground- active brain sources along with the source reconstruction results of three noise learning schemes. For each algorithm, the estimated noise covariance matrix is also plotted above the topographic maps. The source reconstruction performance of these examples terms of EMD and time course correlation (Corr) is summarized in the associated table next to each panel. We also report the accuracy with which the ground-truth noise covariance was estimated in terms of the sim and NMSE. The convergence behaviour all three noise estimation approaches also shown.
+ Hetproaaasstsseeeeeessssccssaai O Full Structure
Full-structure Noise
0.4
Time Course Error
0.3
0.3
F1-Measure
EUCL (mm)
20
EMD
0.2
0.5
0.1
0.1
10
-5
5
10
5
10
0
5
10
SNR(dB)
SNR(dB)
SNR(dB)
SNR(dB)
Heteroscedastic Noise
0.4
0.4
Time Course Error
0.3
0.3
F1-Measure
EUCL (mm)
20
EMD
0.2
0.2
0.5
0.1
0.1
0
-5
0
10
-5
5
10
-5
0
5
10
SNR(dB)
SNR(dB)
SNR(dB)
SNR(dB)
Figure 2: Source reconstruction performance (mean + SEM) of the three different noise learning schemes for data generated by realistic lead field matrix. Generated sensor signals were super- imposed by either full-structure or heteroscedastic noise covering a wide range of SNRs. Perfor- mance was measured in terms of the earth mover's distance (EMD), time-course correlation error, F1-measure and Euclidean distance (EUCL) (mm) between each simulated source and the recon- structed source with highest maximum absolute correlation.
7

---

Under review as a conference paper at ICLR 2021
120 trials
63 trials
12 trials
2 trials
1 trials
L
R
Figure 3: Auditory evoked field (AEF) localization results versus number of trials from rep- resentative subject using FUN learning algorithm. All reconstructions show focal sources at the expected locations in the left (L: top panel) and right (R: bottom panel) auditory cortex. As a result, the limited number of trials does not influence the reconstruction results of FUN learning algorithm.
(green) and full-structure (blue) noise covariances for a range of 10 SNR values. The upper panel represents the evaluation metrics for the setting where the noise covariance is full-structure model, while the lower-panel depicts the same metric for simulated noise with heteroscedastic diagonal co- variance. Concerning the first setting, FUN learning consistently outperforms its homoscedastic and heteroscedastic counterparts according to all evaluation metrics in particular in low-SNR settings. Consequently, as the SNR decreases, the gap between FUN learning and the two other variants increases. Conversely, heteroscedastic noise learning shows improvement over FUN learning according to all evaluation metrics when the simulated noise is indeed heteroscedastic. However, note that the magnitude of this improvement is not as large as observed for the setting where the noise covariance is generated according to a full-structure model and then is estimated using the FUN approach.
Analysis of Auditory Evoked Fields (AEF): Figure shows the reconstructed sources of the Au- ditory Evoked Fields (AEF) versus number of trials from single representative subject using FUN learning algorithm. Further details on this dataset can be found in Appendix [ We tested the re- construction performance of FUN learning with the number of trials limited to 2, 12, 63 and 120. Each reconstruction was performed 30 times with the specific trials themselves chosen as a random subset of all available trials. As the subplots for different trials demonstrate, FUN learning algorithm is able to correctly localize bilateral auditory activity to Heschel's gyrus, which is the characteristic location of the primary auditory cortex, under a few trials even a single trial.
5 DISCUSSION
This paper focused on sparse regression within the hierarchical Bayesian regression framework and its application in EEG/MEG brain source imaging. To this end we developed an algorithm, which is, however, suitable for much wider range of applications. What is more, the same concepts used here for full-structure noise learning could be employed in other contexts where hyperparameters like kernel widths in Gaussian process regression [Wu al. 2019) or dictionary elements in the dic- tionary learning problem (Dikmen & Févotte 2012) are to be inferred. Besides, using FUN learning algorithm may also prove useful for practical scenarios in which model residuals are expected to be correlated, e.g., probabilistic canonical correlation analysis (CCA) (Bach & Jordan 2005), spectral independent component analysis (ICA) (Ablin et al.) 2020). wireless communication (Prasad et al. 2015 Gerstoft al. 2016 Haghighatshoar & Caire 2017 Khalilsarai al. 2020), robust portfo- lio optimization in finance (Feng et al. (016), graph learning (Kumar et al., 2020), thermal field reconstruction (Flinth & Hashemi] 2018), and brain functional imaging (Wei et al. 2020).
8

---

Under review as a conference paper at ICLR 2021
Noise learning has also attracted attention in functional magnetic resonance imaging (fMRI) Cai
eet al., 2016; Shvartsman et al.) 2018 Cai et al. 2019b; 2020b; Wei et al. 2020). where various models like matrix-norma (MN), factor analysis (FA), and Gaussian-prooe processs (GP) regression have been proposed. The majority of the noise learning algorithms in the fMRI literature rely on the EM framework, which is quite slow in practice and has convergence guarantees only under certain strong conditions. In contrast to these existing approaches, our proposed framework not only applies to the models considered in these papers, but also benefits from theoretically proven convergence guarantees. To be more specific, we showed in this paper that FUN learning is an instance of the wider class of majorization majorization-minimizzaoo (MM) framework, for which provable fast convergence is guaranteed. It is worth emphasizing our contribution within the MM optimization context as well. In many MM implementations, surrogate functions are minimized using an iterative approach. Our proposed algorithm, however, obtains a closed-form solution for the surrogate function in each step, which further advances its efficiency.
In the context of BSI, Engemann & Gramfort (2055) proposed method for selecting single regu- larization parameter based on cross-validation and maximum-likeliban estimation, Huizenga et al. ) De Munck et al. (2002); Bijmaa et al. (2003); De Munck et al. (2004; Ahn & Jun (2011); ; Jun et al. (2006) and Plis et al. (2006) assume more complex spatiotemporal noise covari- ance structures. A common limitation of these works is, however, that the noise level is not estimated as part of the source reconstruction problem on task-related data but from separate noise recordings. Our proposed algorithm substantially differs in this respect, it learns the noise covariance jointly with the brain source distribution. Note that The idea of joint estimation of brain source activity and noise covariance has been previously proposed for Type-I learning methods in (Massias al. 2018 Bertrand et al. 2019). In contrast to these Type-I methods, FUN is Type-Il method, which learns the prior source distribution as part of the model fitting. Type-II methods have been reported to yield consistently superior results than Type-I methods (Owen et al., 2012: Cai et al. 2019a 2020a Hashemi et al. 2020). Our numerical results show that the same hold also for FUN learning, which performs on par better than existing variants from the Type-Il II family (including conventional Champagne) in this study. We plan to provide a formal comparison of the performance of noise learning within Type-I and Type-Il estimation in our future work.
While being broadly applicable, our approach is also limited by a number of factors. Although Gaussian noise distributions are commonly justified, it would be desirable to also include more robust (e.g.. heavy-tailed) non-Gaussian noise distributions in our framework. Another limitation is that the superior performance of the full-structure noise learning technique comes the ex- pense of higher computational complexity compared to the variants assuming homoscedastic or heteroscedastic strucutre. Besides, signals in real-world scenarios often lie in a lower-dimensional space compared to the original high-dimensional ambient space due to the particular correlations that inherently exist in the structure of the data. Therefore, imposing physiologically plausible con- straints on the noise model, e.g., low-rank or Toeplitz structure, not only provides side information that can be leveraged for the reconstruction but also reduces the computational cost in two ways: a) by reducing the number of parameters and b) by taking advantage of efficient implementations using circular embeddings and the fast Fourier transform (Babu) 2016). Exploring efficient ways to incorporate these structural assumptions within Riemannian framework is another direction of future work.
6 CONCLUSION
This paper proposes an efficient optimization algorithm for jointly estimating Gaussian regression parameter distributions as well as Gaussian noise distributions with full covariance structure within a hierarchical Bayesian framework. Using the Riemannian geometry of positive definite matrices, we derived an efficient algorithm for jointly estimating source and noise covariances. The benefits of our proposed framework were evaluated within extensive set of experiments in the context of electromagnetic brain source imaging inverse problem and showed significant improvement upon state-of-te oftte----hhe-rr techniques in the realistic scenario where the noise has full covariance structure. The performance of our method is assessed through a real data analysis for the auditory evoked field (AEF) dataset.

---

Under review as a conference paper at ICLR 2021
REFERENCES
Pierre Ablin, Jean-François Cardoso, and Alexandre Gramfort. Spectral independent component analysis with noise modeling for M/EEG source separation. arXiv preprint arXiv:2008. 09693, 2020.
Minkyu Ahn and Sung Jun. MEG/EEG Spatiotemporal Noise Covariance estimation the Least Squares Sense. In The Asia-Paciff Signal and Information Processing Association (APSIPA), 2011.
Prabhu Babu. MELT- www.-mmmii-i---iiiiiiooteooooooo estimation of low-rank Toeplitz covariance matrix. IEEE Signal Processing Letters, 23(11):1577-151, 2016.
Francis R Bach and Michael Jordan. A probabilistic interpretation of canonical correlation analysis. 2005.
Sylvain Baillet, John Mosher, and Richard M Leahy. Electromagnetic brain mapping. IEEE Signal Processing Magazine, 18(6):14-30, 2001.
Yousra Bekhti, Felix Lucka, Joseph Salmon, and Alexandre Gramfort. A hierarchical Bayesian per- spective majorization-ninimizza for non-convex sparse regression: application to M/EEG source imaging. Inverse Problems, 34(8):085010, 2018.
Aharon Ben-Tal. On generalized means and generalized convex functions. Journal of Optimization Theory and Applications, 21(1):1-13, 1977.
Marcel Berger. A panoramic view of Riemannian geometry. Springer Science & Business Media, 2012.
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, pp. 3959-3970, 2019.
Rajendra Bhatia. Positive definite matrices, volume 24. Princeton University Press, 2009.
Fetsje Bijma, Jan De Munck, Hilde M Huizenga, and Rob M Heethaar. mathematical approach to the temporal stationarity of background noise in MEG/EEG measurements. Neurolmage, 20 (1):233-243, 2003.
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(3): 1055-1070, 2009.
Lcon Bungert, Martin Burger, Yury Korolev, and isation for inverse problems with imperfect forward operators and general noise models. arXiv preprint arXiv:2005. /4131, 2020.
Schoenlicb. Variational regular-
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, 2019a.
Chang Cai, Ali Hashemi, Mithun Diwakar, Stefan Haufe, Kensuke Sekihara, and Srikantan S Na- garajan. Robust estimation of noise for electromagnetic brain imaging with the Champagne algo- rithm. Neurolmage, 225: 17411, 2020a.
Ming Bo Cai, Nicolas 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, 1005)e9 2019b.
Ming Bo Cai, Michael Shvartsman, Angi Wu. Hejia Zhang, and Xia Zhu. Incorporating structured assumptions with probabilistic graphical models in fMRI data analysis. Neuropsychologia. pp. 107500, 2020b.
10

---

Under review as a conference paper at ICLR 2021
Mingbo 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, pp. 4951-4959, 2016.
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.
Sebastián Castaño-Candamill Johannes Höhne, Juan-David Martínez- Vargas, Xing-Wei An, Ger- man Castellanos-Dominguz,, and Stefan Haufe. Solving the EEG inverse problem based on space- time- frequency structured sparsity constraints. Neurolmage, 118:598 612, 2015.
Chiao En Chen, Flavio Lorenzelli, Ralph Hudson, and Kung Yao. Stochastic maximum- likellhood
DOA estimation in the presence of unknown nonuniform noise. IEEE Transactions on Signal Processing, 56(7):3038330 0048,,304,220 2008
Naney Chinchor and Beth M Sundheim. Muc-5 evaluation metrics. Fifth Message Understanding Conference (MUC-5), 1993.
Shane F Cotter, Bhaskar D Rao, Kjersti Engan, and Kenneth Kreutz-Delgado. Sparse solutions to linear inverse problems with multiple measurement vectors. IEEE Transactions on Signal Processing, (7777)247722888 2005.
Sarang S Dalal, JM Zumer, V Agrawal, KE Hild, Kensuke Sekihara, and SS Nagarajan. NUTMEG: neuromagnetic source reconstruction toolbox. Neurology & Clinical Neurophysiology: NCN, 2004:52, 2004.
Sarang Dalal, Johanna M Zumer, Adrian Guggisberg, Michael Trumpis, Daniel DE Wong, Kensuke Sekihara, and Srikantan S Nagarajan. MEG/EEG source reconstruction, statistical eval- uation, and visualization with NUTMEG. Computational Intelligence and Neuroscience, 2011, 2011.
Arnak Dalalyan, Mohamed Hebiri, Katia Meziani, and Joseph Salmon. Learning heteroscedastic models by convex programming under group sparsity. In International Conference on Machine Learning, pp. 379-387, 2013.
Jason V Davis, Brian Kulis, Prateek Jain, Suvrit Sra, and Inderjit S Dhillon. Information-theorett metric learning. In Proceedings of the 24th International Conference on Machine Learning, pp. 209-216. 2007.
Z John Daye, Jinbo Chen, and Hongzhe Li. High-dimensional heteroscedastic regression with an application to eQTL data analysis. Biometrics, 68(1):3163366, 2012.
Jan Casper De Munck, Hilde M Huizenga, Lourens Waldorp, and RA Heethaar. Estimating sta- tionary dipoles from MEG/EEG data contaminated with spatially and temporally correlated back- ground noise. IEEE Transactions on Signal Processing, 50(7):1565-1522 2002.
Jan Casper De Munck, Fetsje Bijma, Pawel Gaura, Cezary Andrzej Sieluzycki, Maria Inês Branco, and Rob M Heethaar. A maximum- likelihood MEG/EEG data sets. IEEE Transactions on Biomedical Engineering, (1(1)):213321222 2004.
maximum-likelihood estimator for trial-to-trial variations in noisy
Onur Dikmen and Cédric Févotte. Maximum marginal likelihood estimation for nonnegative dic- tionary learning in the Gamma-Poisson model. JEEE Transactions on Signal Processing, 60(10): 5163-5175, 2012.
David L Donoho. Compressed sensing. IEEE Transactions Information Theory, 52(4): 1306, 2006.
Denis A Engemann and Alexandre Gramfort. Automated model selection in covariance estimation and spatial whitening of MEG and EEG signals. Neurolmage, (08:328-342. 2015.
Dylan Fagot, Herwig Wendt, Cédric Févotte, and Paris Smaragdis. Majorization minimization algo- rithms for convolutive NMF with the beta-divergence. In ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 8202-8206. IEFF, 2019.
Majorization-minnnizzttoon
11

---

Under review as a conference paper at ICLR 2021
Yiyong Feng, Daniel P Palomar, et al. A signal processing perspective financial engineering. Foundations and Trends@ in Signal Processing, (11-))11 2016.
Axel Flinth and Ali Hashemi. Approximate recovery of initial point-like and instantaneous sources from coarsely sampled thermal fields via infiniteddimensional compressed sensing. In 2018 26th European Signal Processing Conference (EUSIPCO), pp. 1720-1724. IEEE, 2018.
Peter Gerstoft, Christoph F Mecklenbräuker, Angeliki Xenaki, and Santosh Nannuru. Multisnapshot sparse Bayesian learning for DOA. IEEE Signal Processing Letters, 23(10): 146100)::466111 1469-1473, 2016.
Irina Gorodnitsky, John S George, and Bhaskar Rao. Neuromagnetic source imaging with FOCUSS: recursive weighted minimum norm algorithm. Electroencephalogguppy and Clinical Neurephysiology. 95(4):231-251) 1995.
Alexandre Gramfort, Matthieu Kowalski, and Matri Hämäläinen. Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods. Physics in Medicine and Biology, 1777):197,,011
Saeid Haghighatshoar and Giuseppe Caire. Massive MIMO channel subspace estimation from low- dimensional projections. IEEE Transactions on Signal Processing, 65(2):303-315 2017.
Matti Hämäläinen, Riita Hari, Risto J Ilmoniemi, Jukka Knuutila, and Olli Lounasmaa. Magne- foencephalography--tee theory, instrumentation, and applications to noninvasive studies of the work- ing human brain. Reviews of modern Physics, 65(2):413, 1993.
Ali Hashemi and Stefan Haufe. Improving EEG source localization through spatio-temporal sparse Bayesian learning. In 2018 26th European Signal Processing Conference (EUSIPCO), pp. 1935- 1939. IEEE, 2018.
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. bioRxiv, 2020.
Stefan Haufe and Arne Ewald. A simulation framework for benchmarking BEG-based brain con- nectivity estimation methodologies. Brain topography, pp. 1-18, 2016.
Stefan Haufe, Vadim V Nikulin, Andreas Ziche, Robur- Müller, and Guido Nolte. Combining sparsity and rotational invariance in EEG/MEG source reconstruction. Neurolmage, 42(2):726~ 738, 2008
Yu Huang, Lucas Parra, and Stefan Haufe. The New York head - precise standardized volume conductor model for EEG source localization and tES targeting. Neurolmage, 140:150-1633 2016.
Hilde M Huizenga, Jan C De Munck, Lourens J Waldorp, and Raoul PPP Grasman. Spatiotemporal EEG/MEG source analysis based on a parametric noise covariance model. IEEE Transactions on Biomedical Engineering, 49(6):533,533 2002.
David R Hunter and Kenneth Lange. A tutorial on MM algorithms. The American Statistician, 58 (1):30-37, 2004.
Matthew W Jacobson and Jeffrey A Fessler. An expanded theoretical treatment of iteration- dependent majirize-minimize algorithms. 2422, 2007.
a IEEE Transactions on Image Processing, 16(10):2411-
Jürgen Jost and Jeurgen Jost. Riemannian geometry and geometric analysis, volume 42005. Springer, 2008
Sung C Jun, Sergey M Plis. Doug M Ranken, and David M Schmidt. Spatiotemporal noise covari- ance estimation from limited empirical magnetoenceehaalggaphic data. Physics Medicine & Biology, 51(21):5549, 2006.
Mahdi Barzegar Khalilsarai, Tianyu Yang, Saeid Haghighatshoar, and Giuseppe Caire. Structured channel covariance estimation from limited samples Massive MIMO. In ICC 2020-2020 IEEE International Conference Communications (ICC), pp. 1-7. IEEF, 2020.
12

---

Under review as a conference paper at ICLR 2021
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(22):1-60, 2020.
Johannes Lederer and Christian L Muller. Don't t fall for tuning parameters: tuning-fiie free variable se- lection in high dimensions with the TREX. In Proceedings of the Twenty-Ninth AAAI Conference on Artificial Intelligence, pp. 2729-2735, 2015.
Tao Li and Arye Nehorai. Maximum likelihood direction finding in spatially colored noise fields using sparse sensor arrays. IEEE Transactions on Signal Processing, 593):1048-1062, 2010.
Leo Liberti. On class of nonconvex problems where all local minima are global. Publications de l'Insiiuut Mathémathique, 76(90):101-109, 2004.
Thomas Lipp and Stephen Boyd. Variations and extension of the convex-concand procedure. Opti- mization and Engineering, 17(2):63-2277, 2016.
Dmitry Malioutov, Müjdat Cetin, and Alan S Willsky. sparse signal reconstruction perspective for source localization with sensor arrays. IEFE Transactions on Signal Processing, 53(8):3010- 3022, 2005.
Mudassir Masood, Ali Ghrayeh, Prabhu Babu, Issa Khalil, and Mazen Hasna. A minorization- maximization algorithm for an-based MIMOE secrecy rate maximization. In 2016 IEEE Global Conference Signal and Information Processing (GlobalSIP), pp. 975-980. IEEE, 2016.
Mathurin Massias, Olivier Fercoq, Alexandre Gramfort, and Joseph Salmon. Generalized concomi- tant multi-task lasso for sparse multimodal regression. In International Conference Artificial Intelligence and Statistics, pp. 1189999-190 1007, 2018.
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.
Eric Mjolsness and Charles Garrett. Algebraic transformations of objective functions. Neural Net- works, 3(6):651-669, 1990.
Maher Moakher. A differential geometric approach to the geometric mean of symmetric positive- definite matrices. SIAM Journal Matrix Analysis and Applications, 26(3):735-777, 2005.
Julia P Owen, David P Wipf, Hagai T Attias, Kensuke Sekihara, and Srikantan Nagarajan. Per- formance evaluation of the Champagne source reconstruction algorithm on simulated and real M/EEG data. Neuroimage, 60(1):305 -233005333,,0 2012.
Diethard Ernst Pallaschke and Stefan Rolewiez. Foundations of mathematical optimization: convex analysis without linearity, volume 388. Springer Science & Business Media, 2013.
Athanase Papadopoulos. Metric spaces, convexity and nonpositive curvature, volume European Mathematical Society, 2005.
Roberto D Pascual-Marqui, Christoph Michel, and Dietrich Lehmann. Low resolution electro- magnetic tomography: new method for localizing electrical activity in the brain. International Journal of psychophysiology, 18(1):49-65, 1994.
Hendrik Bernd Petersen and Peter Jung. Robust instance-optimal recovery of sparse signals at unknown noise levels. arXiv preprint arXiv:2008.08885, 2020.
Peter Petersen, S Axler, and KA Ribet. Riemannian geometry, volume 171. Springer, 2006.
Sergey M. Plis, David M. Schmidt, Sung C. Jun, and Doug M. Ranken. A generalized spatiotemporal covariance model for stationary background in analysis of MEG data. In 2006 International Conference of the IEEE Engineering in Medicine and Biology Society, pp. 3680-3683, IEEE, 2006.
13

---

Under review as a conference paper at ICLR 2021
Ranjitha Prasad, Chandra R Murthy, and Bhaskar D Rao. Joint channel estimation and data detection in MIMO-OFDM systems: A sparse Bayesian learning approach. IEEE Transactions on Signal Processing, 63(20):5369-533) 5382, 2015.
Tamas Rapesak. Geodesic convexity in nonlinear optimization. Journal of Optimization Theory and Applications, 69(1):199-133, 1991.
Meisam Razaviyayn, Mingyi Hong, and Zhi-Quan Luo. A unified convergence analysis of block successive minimization methods for nonsmooth optimization. SIAM Journal on Optimization, 1262):1126-1153, 2013.
Yossi Rubner, Carlo Tomasi, and Leonidas Guibas. The earth mover's distance as metric for image retrieval. International Journal of Computer Vision, 40(2):99-121, 2000.
Matthias W Seeger and David Wipf. Variational Bayesian inference techniques. IEEE Signal Processing Magazine, 27(6):81-91, 2010.
Kensuke Sekihara and Srikantan S Nagarajan. Electromagnetic brain imaging: a Bayesian perspec- rive. Springer, 2015.
Michael Shvartsman, Narayanan Sundaram, Mikio Aoi, Adam Charles, Theodore Willke, and Jonathan Cohen. Matrix-normal models for fMRI analysis. In International Conference on Arti- ficial Intelligence and Statistics, pp. 1914-1923. PMLR, 2018.
Suvrit Sra and Reshad Hosseini. Conic geometric optimization on the manifold of positive definite matrices. SIAM Journal on Optimization, 25(1):713-739, 2015.
Suvrit Sra and Reshad Hosseini. Geometric optimization in machine learning. In Algorithmic Advances in Riemannian Geometry and Applications, pp. 73-91. Springer, 2016.
Ying Sun, Prabhu Babu, and Daniel Palomar.
algorithms in signal
processing, communications, and machine learning. IEEE Transactions on Signal Processing, 65 19::79-8116, 2017.
Michael E Tipping. Sparse Bayesian learning and the relevance vector m Journal of Machine Learming Research, 1(Jun):211 244, 2001.
Sara Van de Geer, Johannes Lederer, al. The Lasso, correlated design, and improved oracle in- equalities. In From Probability to Statisties and Back: Hfigh-Dimensional Models and Processes- Festschrift in Honor of Jon A. Wellner, pp. 303-316. Institute of Mathematical Statisties, 2013.
CJ van Rijsbergen. Information retrieval. 1979.
Nisheeth K Vishnoi. Geodesic convex optimization: Differentiation on manifolds, geodesics, and convexity. arXiv preprint arXiv:1806.06377. 2018.
Huilin Wei, Amirhossein Jafarian, Peter Zeidman, Vladimir Litvak, Adeel Razi, Dewen Hu, and Karl J Friston. Bayesian fusion and multimodal DCM for EEG and fMRI. Neurolmage, 211: 116595, 2020.
Ami Wiesel, Teng Zhang, et al. Structured robust covariance estimation. Foundations and Trends@ in Signal Processing, 8(3):127-216, 2015.
David Wipf and Srikantan Nagarajan. A unified Bayesian framework for MEG/BEG source imaging. Neurolmage, 44(3):947-966, 2009.
David P Wipf and Bhaskar D Rao. Sparse Bayesian learning for basis selection. IEEE Transactions on Signal Processing, 52(3):1133-214, 2004.
David P Wipf and Bhaskar D Rao. An empirical Bayesian strategy for solving the simultancous sparse approximation problem. IEEE Transactions Signal Processing, 55(7):3704-3711 2007.
David P Wipf, Julia P Owen, Hagai T Attias, Kensuke Sekihara, and Srikantan S Nagarajan. Robust Bayesian estimation the location, orientation, and time course of multiple correlated neural sources using MEG. Neurolmage, 49(1 (111))641-655. 2010.
14

---

Under review as a conference paper at ICLR 2021
Anqi Wu, Oluwasanmi Koyejo, and Jonathan W Pillow. Dependent relevance determination for smooth and structured sparse regression. J. Mach. Learn. Res., 20(89): 43, 2019.
Tong Tong Wu, Kenneth Lange, et al. The MM alternative to EM. Statistical Science, 25(4))492- 505, 2010.
Wei Wu, Srikantan Nagarajan, and Zhe Chen. Bayesian machine learning: EEG\,EG signal pro- cessing measurements. IEEE Signal Processing Magazine, 33(1):14-36, 2016.
Alan Yuille and Anand Rangarajan. The concave-convex procedure. Neural computation, 15(4): 915-936, 2003.
Pourya Zadch, Reshad Hosseini, and Suvrit Sra. Geometric mean metric learning. In International Conference on Machine Learning, pp. 2464-2471, 2016.
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-926, 2011.
A PROOF OF THEOREM
Proof. We start the proof by recalling equation
T C² (T, = - log Y = logg Σy + Σ ( E Zyy'(())())
(16)
The upper bound on the log Eyl term can be directly inferred from the concavity of the log- determinant function and its first-order Taylor expansion around the value from the previous iter- ation, Dÿ. which provides the following inequality (Sun et al. 2017 Example 2):
log Eyy S log
Σ
tr
(>)) Ey
-x/]]
=
log Σ + tr (x+))
-1 Σy tr (
>
(17)
Note that the first and last terms in equation 17 do not depend on r; hence, they can be ignored in the optimization procedure. Now, we decompose My into two terms, each of which only contains either the noise or source covariances:
tr [ Σy tr -
-1 ( (A LA+
(M+)
A
tr
LLL
(18)
In next step, we decompose the second term in equation 8 each of which function of either only the noise only the source covariances. To this end, exploit the following relationship between sensor and source space covariances:
ET y(t) Ey yt),. into two terms,
T y(t) my y(t)
1 Σ
Tt T t=1
x
1 - (v(t)) -xxt)-Lxx((8 (t)) A (y((t) Lx" (t))
Lx)(
(19)
By combining equation 18 and equation 19 depend on I, we have:
rearranging the terms, and ignoring all terms that do not
(T(T)) < tr
T X (t) + const
T
= =Cr((C$) r' tr ( + const = Leouve (T, A A + const
(20)
Note that constant values in equation 20 do not depend r: hence, they be ignored in the optimization procedure. This
where C2 (x+) ( and M = Σt x
15

---

Under review as a conference paper at ICLR 2021
proves the equivalence of equation 8 and equation 9 when the optimization is performed with respect
to r.
The equivalence of equation 8 and equation 11 can be shown analogously, with the difference that
we only focus on noise-related terms in equation 18 and equation T9
(A) - tr (2) A
(1()) Lx
Lxx(x)) (t)) A - (y(() Lx (t) const
-Ir(( A) +
const Cnowe (T* A) + const
(21)
where C Σy. and M = (y - Lx (x)) - Lx* (t)) Constant values in
equation do not depend on A; hence, they can again be ignored in the optimization procedure. Summarizing, we have shown that optimizing equation 8 is equivalent to optimizing Luoie (T A)
and S (T,*) which concludes the proof.
B PROOF OF THEOREM 2
Before presenting the proof, the subsequent definitions and propositions are required:
Definition 4 (Geodesic path). Let M be Riemannian manifold, i.e., differentiable manifold whose tangent space endowed with an inner product that defines local Euclidean structures. Then, geodesic between points M, denoted by Po, P1 .M, is defined the shortest connecting path between those two points along the manifold, (i P1) E M for LE [0,1], 1], where I = and 1 -1 defines the starting and end points of the path, respectively.
In the current context, (Po, P1) defines a geodesic curve on the positive definite (P.D.) mani- fold joining two P.D. matrices, Po, P1 > The specific pairs of matrices will deal with are {C M and (
Definition 5 (Geodesic on the P.D. manifold). Geodesics on the manifold of P.D. matrices can be shown to form a cone within the embedding space. We denote this manifold by S44. Assume two P.D. matrices Po, P1 E S44. Then,for [0, 1]. the geodesic curve joining Po to P. is defined as
Bhutia 2009
Chupter. 6):
(i 0. P1) = (Po) ( (PP)) (Po) 1E [0,
(22)
Note that Po and P1 are obtained as the starting and end points of the geodesic path by choosing I 0 and / = 1, respectively. The midpoint of the geodesic, obtained by setting = }.2 is called the geometric mean. Note that, according to Definition ()) the following equality holds
(T - (To) (IT)) (
(TO))*)
(To) 1/2
(To) 1/2) -
(23)
Definition 6 (Geodesic convexity). Let po and P1 be arbitrary points on subset A of
Riemannian manifold M. Then real-valued function f with domain c. M with f A + R is called geodesic convex (g-convex) if the following relation holds:
f (C((Po-P)) S Iff(p)) + +11
(24)
where E j,, 1] and C(po. pi) denotes the geodesic path connecting two points Po and P1 as defined in 4 Thus, analogy classical convexity, the function f is g-convex if every geodesic ((po. P1) of M between Pu, E A, lies in the g- g-convex set A. Note that the set A C MM is called g-convex, any geodesics joining an arbitrary pair of points lies completely in A.
Remark 7. Note that g-convexity is generalization of classical (linear) convexity non non-Euclidean (non-linear) geometry and metric spaces. Therefore, is straightforward to show that all convex functions in Euclidean geometry also g-convex, where the geodesics between pairs of matrices are simply line segments:
C((po, PI) Ipo 1)p
(25)
16

---

Under review as a conference paper at ICLR 2021
For the sake of brevity, omit detailed theoretical introduction of g-convexity, and only borrow
a result from Zadeh et al. (2016); Sra & Hosseini (2015). Interested readers are referred to Wiesel
et al.. (2015 Chapter 1) for gentle introduction to this topic, and Papadopoulos (2005 Chapter. 2)
(1 Ben- (977); Liberti 2004); |Pallaschke & Rolewicz (203)): Bonnabel & Sepul-)
chre (2009);
Moakher (2005): Sra & Hosseini (2016): Vishnoi
(2sh for more in-depth technical
details.
Now we are ready state the proof, which parallels the one provided in Zadeh el. 2016 Theorem. 3).
Proof. We only show the proof for Lource (T, Ax). The proof for n (T A) can be presented analogously; and therefore, is omitted here for brevity. We proceed in two steps. First, we limit our attention to P.D. manifolds and express equation 24 in terms of geodesic paths and functions that lie this particular space. We then show that L (T, AAx)) is strictly g-convex this
specific domain. the second step, we then derive the updates rules proposed in equation and equation 14
B. PART I: PROVING G-CONVEXITY THE MAJORIZING COST FUNCTIONS
We consider geodesics along the P.D. manifold by setting ((Po. P1) to (T F1) as presented in Definition 5 and define f(.) be f(T) tr(CkT) tr )) representing the cost function
R.
source
We now show that ((T) strictly g-convex on this specific domain. For continuous functions as considered in this paper, fulfilling equation 24 for fTT) and (To, [1) with 1 = 1/2 is sufficient to prove strict g-convexity:
tr Sy (Tu [1)) tr (M 11)
< tr (C&To) + z (MSTO ") + tr (CTT,) + 2th (M
2
(26)
Given C E S++, i.c, C2 > 0 and the operator inequality Bhatia 2009 Chapter. 4)
Cya(To, (1) 2 Ti,
(27)
we have:
tr (C (To, 2th (CHTo) + (C&TT)
(28)
2
which is derived by multiplying both sides of equation 27 with C followed by taking the trace on both sides.
Similarly, we can write the operator inequality for [To² I using equation 23 as:
E/2 (To E1/2 (T 1 r-1 r-
(29)
Multiplying both sides of equation 29 by M E S44. and applying the trace operator on both sides leads to:
ir (M -1) + 1 - (M
tr MMEEyaaToo (T)
(30)
Summing up equation 28 and equation 30 proves equation 26 and concludes the first part of the proof.
B.2 PART II: DETAILED DERIVATION OF THE UPDATE RULES IN EQUATIONS 13 AND 14
We now present the second part of the proof by deriving the update rules in equations I and 14 Since the cost function Leone (T, Ag) strictly g-convex, its optimal solution in the k t-th iteration
17

---

Under review as a conference paper at ICLR 2021
is unique. More concretely, the optimum can be analytically derived by taking the derivative of equation 9 and setting the result to zero as follows:
Source (T,As) (Ck) - T- 0,
(31)
which results in
r (C M
(32)
This solution is known as the Riccati equation, and is the geometric mean between Cg and M Davis et al. 2007 Bonnabel & Sepulchre 2009):
r (C+) ( (CM))
(C))
The update rule for the full noise covariance matrix can be derived analogously:
A = (C*)
(
Remark 8. Note that the obtuined update rules are closed-form solutions for the surrogate cost functions, equations 9 and 17, which stands contrast conventional majorization minimization algorithms (see section in the appendix), which require iterative procedures in each step of the optimization.
Deriving the update rules in equation 13 and equation 14 concludes the second part of the proof of Theorem
C PROOF OF THEOREM 3
In the following, we provide proof for Theorem 3 by showing that alternating update rules for I and A in equation 13 and equation 14 are guaranteed to converge to local minimum of the Bayesian Type-II likelihood equation the general class of majorization mitiiiiitoommmiimiattioon (MM) algorithms, for which this property follows by construction. To this end, we first briefly review theoretical concepts behind the majorization- minimization (MM) algorithmic framework (Hunter & Lange 2004 Razaviyayn et al. 2013 Ja- cobson & Fessler 2007 Wu et al., 2010)).
In particular, we will prove that FUN learning is an instance of
C.I REQUIRED CONDITIONS FOR MAJORIZATION-MINIIIZATTII ALGORITHMS
MM encompasses family of iterative algorithms for optimizing general non-linear cost functions. The main idea behind iM to replace the original cost function in each iteration by upper bound, also known as majorizing function, whose minimum easy to find. The MM class covers broad range of common optimization algorithms such convex-concave procedures (CCCP) and proximal methods (Sun et al. 2017 Section IV), (Mjolsness & Garrett 1990 Yuille & Rangarajan 2003 Lipp & Boyd (016). Such algorithms have been applied in various domains such brain source imgging Haahemi & Haufe 2018; Bekhhii et al.] 2018 Cai et al. 2020a Hashemi et al. 2020), wireless communication systems with massive MIMO technology (Masood al. 2016]] Haghighatshoar & Caire 2017: Khalilsarai et al. (2020). and non-negative matrix factorization (Fagot et al. 2019). Interested readers are referred to Sun et al. (2017) for an extensive list of applications on MM.
The problem of minimizing a continuous function f(u) within a closed convex set C R":
min f(u) subject to u e uu
(33)
within the MM framwork can be summarized as follows. First, construct a continuous surrogate function g(ula²) that majorizes, upper-bounds, the original function fu) and coincides with fu) at given point s
[A1]
uu f(u²)
Yu Vu Eu.
[A2]
g(ulu*) 2 fu)
18

---

Under review as a conference paper at ICLR 2021
Second, starting from initial value u². generate
sequence of feasible points ...., u u ++l as solutions of a series of successive simple optimization problems, where
[A3]
u ++1
arg min g(ulu²). uel
If a surrogate function fulfills conditions [A1]-[3], then the value of the cost function f decreases each iteration: f( S f(u²). For the smooth functions considered this paper, we further require that the derivatives of the original and surrogate functions coincide at u²:
[A4]
Vg(u"|u>> vf("²)
uu Eu.
We can then formulate the following theorem:
Theorem 9. Assume that an MM algorithm fulfills conditions [A]]-[4]. Then, every limit point of the sequence of minimizers generated in [A3], is a stationary point of the original optimization problem in equation 33
Proof. A detailed proof is provided in Razaviyayn et al. 2013 Theorem 1).
C.2 DETAIL DERIVATION OF THE PROOF OF THEOREM
3
We now show that FUN learning is an instance of majorization-mininizztooo defined above, which fulfills Theorem 9
Proof. We need to prove that conditions [A1]-[A4] are fulfilled for FUN learning. To this end, we recall the upper bound on log Eyy in equation 17 which fulfills condition [A2] since it majorizes log |yyl as a result of the concavity of the log-determinant function and its first-order Taylor ex- pansion around Myy Besides, it automatically satisfies conditions [AI] and [A4] by construction, because the majorizing function in equation 17 is obtained through Taylor expansion around Σ/. Concretely, [A1] satisfied because the equality in equation 17 holds for Ey Σ Similarly, [A4] satisfied because the gradient of log EEy at point Σý, (Σş) defines the linear Taylor approxi- mation log b + tr (x*) -1 (Ey - } Thus, both gradients coincide in M by construction. Now, we prove that [A3] can satisfied by showing that C (T,A*) reaches its global min- imum in each MM iteration. This guaranteed if Leouree (T, A*) can be shown to be convex or g-convex with respect to r. To this end, we first require the subsequent proposition:
Proposition 10. Any local minimum of g fonvettticction over a g-convex set a global minimum.
Proof. A detailed proof is presented in Rapesak 1991. Theorem 2.1).
Given the proof presented in appendix B.I B.I we can conclude that equation 20 g-convex; hence, any local minimum of Ceouv (T) A is global minimum according to Proposition 10 This proves that condition [A3] fulfilled and completes the proof that the optimization of equation X with respect to I using the convex surrogate cost function equation 9 leads to an MM algorithm. For the sake of brevity, we omit the proof for the optimization with respect to A based on the convex surrogate function in equation Ceonv (T' AA), it can be presented, analogously.
D DERIVATION OF CHAMPAGNE AS A SPECIAL CASE OF FUN LEARNING
We start the derivation of update rule equation 15 by constraining I' to the set of diagonal matrices W: r diag(y), where 7 l?...... 7Nv] TT We continue by rewriting the constrained optimization with respect to the source covariance matrix,
p arg min (r(CTT) +r((MT rew, A=
(34)
19

---

Under review as a conference paper at ICLR 2021
as follows:
*+| arg min 7. A=
diag (CS) -1 | 7 diag M y-
(35)
Cdias
(7/l1/*)
where -- =
by ..... IN IT is defined as the element-wise inversion of y. The optimization with respect the scalar source variances is then carried out by taking the derivative of equation 33 with respect to 7n. for nt = 1, ....,, N, and setting it to zero,
a ([(C))
I In + [M mo
8m
[(C
-1 n.n
M
(%))
for n -I...,,N
where L, denotes the n-th column of the lead field matrix. This yields the following update rule
M
E (xh L E3) Ln
for 1, I....,N.
(c))
which is identical to the update rule of Champagne (Wipf & Nagarajan 2009).
E DERIVATION OF CHAMPAGNE WITH HETEROSCEDASTIC NOISE LEARNING AS SPECIAL CASE OF FUN LEARNING
Similar to Appendix D we start by constraining A to the set of diagonal matrices W: A = diag( (A), where À [Ai,.. As]T. IT We continue by reformulating the constrained optimization with respect to the noise covariance matrix,
A A+1
arg min tr(Ck + (M
(36)
AEW, r=
follows:
A k++ arg min
diag [C*) -11 "I + diag [M]] X
(37)
A. P=P
Aeliog (A[**)
where A-1 [X]] ..... IT is defined the element-wise inversion of A. The optimization with respect to the scalar noise variances then proceeds by taking the derivative of equation 37 with respect Amy for 1, ....., M, and setting it zero,
([(ck) "] Am + [M]]
Xm)
axm.
[(kk)
-1
1
[M m.m
(Am) )2
m,m
0 for m I,..,,M
M
This yields the following update rule:
1 1x(()( - (xx( [x)
M
T
k++
m
(Σ))
m.m
for m 1....,M M
(38)
which is identical to the update rule of the Champagne with heteroscedastic noise learning as pre- sented in Cai et al. (2020a).
20

---

Under review as a conference paper at ICLR 2021
Figure 4: Accuracy of the noise covariance matrix reconstruction in- curred by three different noise learn- ing approaches assuming homoscedas- tic (red), heteroscedastic (green) and full-structure (blue) noise covariances. The ground-trutt noise covariance ma- trix is either full-structure (upper row) or heteroscedastic diagonal (lower row). Performance is assessed in terms of the Pearson correlation between the en- tries of the original and reconstructed noise covariance matrices, A and Å, de- noted by Asim (left column). Shown the similarity error 1 A sim, Fur- ther, the normalized mean squared error (NMSE) between A and A, defined as NMSE d| - Al////|||/ is reported (right column).
- Homoscedastic - Heteroscadastic o Full Structure
Full-structure Noise
sim
NMS
0.5
0.5
0
0
-6
-4
-2
-6
-4
-2
SNR(dB)
SNR(dB)
Heteroscedastic Noise sm
NMSE
00.5
0.5
0
0
-6
-4
-2
-6
-4
-2
SNR(dB)
SNR(dB)
F PSEUDO-EEG SIGNAL. GENERATION
Our simulation setting is an adoption of the EEG inverse problem, where brain activity is to be reconstructed from simulated pseudo-EEG data (Haufe & Ewald 2016))
Forward Modeling: Populations of pyramidal neurons in the cortical gray matter are known be the main drivers of the EEG signal (Hämäläinen et al.) 1993: Baillet et al. 200T). Here, we use realistic volume conductor model of the human head model the linear relationship between pri- mary electrical source currents generated within these populations and the resulting scalp surface potentials captured by EEG electrodes. The lead field matrix, E R58x2004 was generated using the New York Head model (Huang et al. 2016) taking into account the realistic anatomy and elec- trical tissue conductivities of an average human head. In this model, 2004 dipolar current sourees were placed evenly on the cortical surface and 58 sensors were considered. The lead field matrix, E IR R58x2004 was computed using the finite element method. Note that the orientation of source currents was fixed to be perpendicular to the cortical surface, so that only scalar source amplitudes needed to be estimated.
Evaluation Metries: Source reconstruction performance was evaluated according to the following metrics. First, the earth mover's distance (EMD) (Rubner et al.] 2000; Haufe et al. 2008) was used quantify the spatial localization accuracy. The EMD measures the cost needed to transform two probability distributions defined on the same metric domain (in this case, distributions of the true and estimated sources defined in 3D Euclidean brain space) into each other. EMD scores were normalized to J0, 1]. Second, the error in the reconstruction of the source time courses was measured. To this end, Pearson correlation between all pairs of simulated and reconstructed (i.e.. those with non-zero activations) source time courses was assessed as the mean of the absolute correlations obtained for each source, after optimally matching simulated and reconstructed sources based on maximal absolute correlation. We also report another metric for evaluating the localization error as the average Euclidean distance (EUCL) (in mm) between each simulated source and the best (in terms of absolute correlations) matching reconstructed source. For assessing the recovery of the true support, we also compute F1 Fi-meauure scores (Chinchor & Sundheim 1993 van Rijsbergen 1979): where P denotes the number of true active sources, while TP and HP are the
F /xTT///TTPPP. numbers of true and false positive predictions. Note that perfect support recovery, i.c., F, = 1, is only achieved when there is perfect correspondence between ground-truth and estimated support.
21

---

Under review as a conference paper at ICLR 2021
To evaluate the accuracy of the noise covariance matrix estimation, the following two metrics were calculated: the Pearson correlation between the original and reconstructed noise covariance ma- trices, A and A, denoted by A sim, and the normalized mean squared error (NMSE) between A and Å, defined as: NMSE = І| A II Similarity error was then defined as one minus the Pearson correlation:
A Asim, Note that NMSE measures the reconstruction of the true scale of the noise covariance matrix, while A sim is scale-invariant and hence only quantifies the overall structural similarity between simulated and estimated noise covariance matrices.
Evaluating the accuracy of the noise covariance matrix estimation: deepees the accuracy with which the covariance matrix is reconstructed by three different noise learning approaches as- suming noise with homoscedastic (red), heteroscedastic (green) and full (blue) structure. The ground truth noise covariance matrix either had full (upper row) or heteroscedastic (lower row) structure. Performance was measured terms of similarity error and NMSE. Similar to the trend observed in Figure 2, full-structure noise learning leads to better noise covariance estimation accuracy (lower NMSE and similariy error) for the full-structure noise model, while superior reconstruction perfor- mance is achieved for heteroscedastic noise learning when true noise covariance is heteroscedastic.
G FURTHER DETAILS ON AUDITORY EVOKED FIELDS (AEF) DATASET
The MEG data used in this article were acquired in the Biomagnetic Imaging Laboratory at the University of California San Francisco (UCSF) with a CTF Omega 2000 whole-head MEG system from VSM MedTech (Coquitlam, BC, Canada) with 1200 Hz sampling rate. The lead field for each subject was calculated with NUTMEG (Dalal et al. 2004) using single-sphere head model (two spherical orientation lead fields) and 8 mm voxel grdd. Each column normalized to have norm of unity. The neural responses one subject to an Auditory Evoked Fields (AFF) stimulus localized. The AEF response was elicited with single 600 ms duration tones kHz) presented binaurally. 120 trials were collected for AEF dataset. The data were first digitally filtered from 1 to 70 Hz remove artifacts and DC offset, time-aligned to the stimulus, and then averaged across the following number of trials:(1,5 ms and the stimuuu time window was selected to be 60 ms to 180 ms, where 0 the onset of the tone Wipf et al. 2010 Dalal et al. Owen
63, 120). The pre-stimulus window was selected to be 100 ms
Owen et al. 2012 Cal et al. 2019a).
22


Reviews:
[{"review_id": "Vh2i8ncYCg4", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 4, "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The authors propose a methodology for type-II maximum likelihood on a hierarchical Bayesian model for EEG signals. The particular feature of the model, which separates it from other EEG models, is the consideration of a full covariance matrix which makes the noise correlated and heteroskedastic. \n\nThe model, as claimed by the authors, is fully Gaussian and therefore tractable. As a consequence, the inference poses no challenges other than the computational complexity. To address this, the authors propose a mechanism for, what they claim is, efficient optimisation. This contribution alone is not sufficient (over the standard literature) for publication as a theoretical improvement. \n\nGiven the lack of a theoretical advancement, I was hoping that the contribution of the article came in the experimental treatment, however, it was not the case.  A single set of experiments using synthetic data was considered, where the proposed method was compared against other benchmark. It is far form surprising when the authors deal with exact inference on a model where the observations where produced under the same statistical assumptions.\n\nI also would like to emphasise that the discussion of the paper states that  \"This paper proposes an efficient optimization algorithm for jointly estimating....\" and \"The benefits of our proposed framework were evaluated within an extensive set of experiments \". None of these claims are true or at least they not validated by any supporting evidence in the paper. \n\nPerhaps with the stated future work and stronger experimental results (real data), this paper can be improved. \n\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "Gwg2uPpqiwT", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 5, "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Joint Learning of Full-structure Noise in Hierarchical Bayesian Regression Models\n\nSummary:\n\nThe paper argues that modeling the full covariance structure in a sparse bayessian learning setting leads to significantly better results in eeg inverse problems. The paper details a majorization-minimization type algorithm leading to a set of fairly simple update rules. The proposed method is evaluated on simulated data.\n\nPositive:\n\n1. The proposed method is well motivated and the problem is highly relevant.\n2. The mathematical details regarding the algorithm are presented in sufficient detail.\n3. The paper is well written and easy to follow for the most part.\n4. Experiments are reasonable and presented clearly.\n\nNegative:\n\n1. The abstract could be improved to more clearly describe the problem and contributions of the paper in a self-contained manner.\n2. Has this particular problem (sparse bayesian regression with full covariance noise) not been considered by others? The main contribution, in my view, is algorithmic; which other algorithms have been used previously for this type of problem? (I think both ML-II and MCMC and possibly other methods have previously been used.) I would have liked a review and experimental comparison.\n3. While the experimental evaluation is reasonable, I think the paper would benefit from a demonstration and benchmarking with competing approaches on a real data task.\n4. Experiments on simulated data highlighting more clearly the *algorithmic* advantages of the proposed method would be appreciated.\n5. I did not notice a link to software implementing the proposed method? Sharing software implementations will significantly strengthen the contribution and allow the community to reproduce the results.\n\nRecommendation:\n\nWeak reject.\n\n\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "9AN9CUTcTfF", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 4, "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Summary:\nThe paper generalised Type-II ML regression models for scenarios where different noise dimensions cannot be assumed independent, but instead one needs to model the full covariance structure. This is clearly an important problem and it is well motivated in the work.\n\nReasons for score:\nI recommend rejecting the paper even though it represents high-quality work in statistics, because I think it is somewhat tangentially related to ICLR and the contribution would be better appreciated in a different venue.\n\nStrong points:  (1) Addresses an important problem. (2) Seems to work well in practice\n\nWeaknesses: (1) Limited conceptual novelty. (2)Technical contribution hidden in Appendix\n\nDetailed review:\nThe work addresses a relevant statistical question of accounting for correlated noise in hierarchical linear regression, but feels somewhat of a poor fit for ICLR. It formally fits within the scope, but still feels out of place in the sense that neither readers interested in the theoretical contributions nor people looking to apply these methods would consider ICLR as a natural venue to look for the information. The development is restricted to a specific, relatively simple, model family that is frequently used in several fields but that is not at the core of the ICLR community. This is highlighted also by the fact that the technical contribution is largely in statistical properties of the covariance estimator, and for this audience gets hidden in the Appendix. Consequently, I believe that paper would much more naturally fit into a publication forum in statistics.\n\nThe proposed approach itself is sound and well developed. Accounting for correlated noise is a very obvious thing to do, but the technical details are non-trivial. The authors rely on Riemannian optimisation for covariance matrices and are able to use the recent Champagne algorithm for SBL. The detailed derivation of Theorem 2 shows non-trivial technical contribution, but remains somewhat isolated as it is hidden in the Appendix. For example, there is no discussion on whether the result derived here would have uses also in other model families. I can see several potential uses for better tools for learning full covariance noise e.g. in matrix factorisation models (e.g. probabilistic CCA relies on covariance estimates) or non-linear regression models, but the authors do not discuss this at all. A proper discussion on this would be important to link the work more closely to the broader activities in the field, to extend the contribution beyond the current viewpoint of a very specific model.\n\nThe empirical experiments are well carried out and demonstrate the value of learning the full covariance matrix compared to methods that only operate with diagonal noise. This is sufficient, since no clear comparison methods accounting for full covariance are available.", "extra_scores": {}, "extra_text": {}}, {"review_id": "-aG7bvrGso", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 4, "confidence": 5, "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper proposes an efficient optimization method for estimating the full noise covariance in a hierarchical Bayesian framework. It's shown in the experiment that the optimization method could recover the true noise covariance in a simulated example and estimating the full covariance has better performance than homo- and heteroscedastic covariance.\n\nI think the proposed method is an effective tool to estimate the full noise covariance especially for the problem setting in this paper. But the overall novelty and contribution are not strong enough for the ICLR community.\n\nPapers in fMRI literature [Michael Shvartsman et al 2017, Anqi Wu et al 2019] have proposed to work with full noise covariance in more complicated models such as factor analysis, Gaussian process regression. The basic model in this paper is a bit too simple compared with other models preventing from making significant methodological contributions. It might fit a signal processing or brain source imaging specialized publication better.\n\nAlso in many applications (especially with brain data), it's shown that a full rank noise covariance is not always preferable given that there are usually some correlations among measurements that lead to lower dimensional subspace at the noise level. So I'm not quite sure whether a full covariance without any structural or subspace assumption would really outperform low-rank full covariance when applying to the real data.\n\nAnother issue in this paper is there is no real data application. I'm not very convinced that simulated data generated from a realistic lead field matrix is considered as the real-world data. \n", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper presents hierarchical Bayesian methods for modelling the
full covariance structure in cases where noise dimensions cannot be
assumed independent.

This is an important problem with potential practical importance. The
work is solid.

Conceptual novelty in the work is somewhat limited.

The method is applied in the paper on hierarchical linear
regression. It is claimed to be applicable to other methods as well,
and the claim is plausible, but to be fully convincing, results and
comparisons would need to be shown. The new extended discussion does
help somewhat.

There was also discussion about whether ICLR is the best match for
this work. This is not a strereotypical ICLR paper though is relevant.

Authors are encouraged to continue this line of work.

Decision: Reject