# Combining spatial wavelets and sparse Bayesian learning for extended brain sources reconstruction

Samy Mokhtari, Jean-Michel Badier, Christian G. Bénar, Bruno Torrésani

<figure><img src="image_1.png" /></figure>

## To cite this version:

Samy Mokhtari, Jean-Michel Badier, Christian G. Bénar, Bruno Torrésani. Combining spatial wavelets and sparse Bayesian learning for extended brain sources reconstruction. 2024. hal-04558288

HAL Id: hal-04558288

https://amu.hal.science/hal-04558288

Preprint submitted on 24 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Combining spatial wavelets and sparse Bayesian learning for extended brain sources reconstruction

Samy Mokhtari, Jean-Michel Badier, Christian G. Bénar, Bruno Torrésani

► To cite this version:

Samy Mokhtari, Jean-Michel Badier, Christian G. Bénar, Bruno Torrésani. Combining spatial wavelets and sparse Bayesian learning for extended brain sources reconstruction. IEEE Transactions on Biomedical Engineering, In press, 10.1109/TBME.2025.3629010. hal-04558288v3

HAL Id: hal-04558288

https://amu.hal.science/hal-04558288v3

Submitted on 7 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Combining spatial wavelets and sparse Bayesian learning for extended brain sources reconstruction

Samy Mokhtari, Jean-Michel Badier, Christian G. Bénar, Bruno Torrésani

**Abstract—**Objective: the accurate reconstruction of extended cortical activity from M/EEG data is a difficult, ill-conditioned problem. This work proposes to model unknown sources as expansions on a wavelet system defined on the cortical surface, and addresses resulting numerical optimization problems. The objective is to obtain accurate source localization, together with quantitatively relevant amplitude and time course.

**Approach:** Unknown sources are expanded on a system of spectral graph wavelets (SGW) defined on the cortical surface. Unknown wavelet coefficients are estimated using either variational or Bayesian formulations, involving priors that favor extended source through sparsity in the wavelet domain: sparsity-inducing regularization, or sparse Bayesian learning (SBL). These approaches are tested and compared with concurrent approaches on numerical simulations. The quality of reconstructions is assessed using a set of complementary metrics.

**Results:** SGW-based approaches are able to identify accurately extended sources. The combination with SBL is particularly attractive, as it doesn't involve hyperparameter tuning. It yields accurate and robust results with respect to all considered metrics, and performs remarkably well in terms of depth bias.

**Conclusion:** This paper demonstrates the usefulness of cortical wavelets for reconstructing cortical activity from M/EEG data, and the impact of sparse Bayesian learning in this context.

**Significance:** Being able to identify localization, depth and time course of brain activity from M/EEG data is important in clinical applications such as epilepsy, as it can improve the detection of potential sources of seizures.

**Index Terms**—M/EEG inverse problem, extended sources, spectral graph wavelets, Sparse Bayesian Learning.

## I. INTRODUCTION

Magneto- and Electro-encephalography (M/EEG) are the main (non-invasive) imagery modalities which allow, thanks to their high temporal resolution, physicians and clinicians to detect and localize cortical activities. From a given set of sensors and measurements (which contain cortical activity, brain background activity and sensor noise), the aim is to reconstruct the sources activity. This is done by solving the M/EEG inverse problem [1], [2], [3]. To address this challenge, multiple methods have been designed, generally classified in parametric versus imaging/distributed methods. While the former considers a few equivalent current dipoles of unknown locations, the latter considers a fine spatial discretization of the

brain cortical surface, with ~ $10^4$ unknown cortical sources with known and fixed locations. While inherently leading to very under-determined inverse problems, distributed methods can be expected to be better suited to large cortical activities. However, results tend to depend strongly on prior information, which is injected into the chosen approach through regularization (in variational methods) or prior (in Bayesian approaches).

The goal of this paper is to introduce and study variational and Bayesian approaches which rely on a spatial transformation, and aim at promoting extended solutions compressed in a relevant transformed domain. We particularly focus on a spatial wavelet transform which, combined with sparsity-promoting regularizations/priors is indeed able to produce extended solutions. While variational approaches suffer, as often, from the difficulty of choosing appropriate values for hyperparameters, we show that the integration of wavelet-based regularization into the sparse Bayesian learning (SBL) paradigm yields more robust solutions. Wavelet-based solutions had already been studied in [4], where spherical wavelet bases were used in the framework of variational approaches. We rather rely here on spectral graph wavelets [5], which provide expansions on an overcomplete spatial wavelet family, and investigate corresponding variational and SBL solutions. These solutions are compared with reference solutions using extensive numerical simulations designed to mimic real evoked potentials. Performance assessment is done through a set of metrics, which include standard ones (implemented in specialized computer packages such as MNE-Python [6], which we therefore use) as well as new ones, which we found convenient in the context of sparse solutions.

Besides methodological contributions (in particular the use of cortical wavelet frames and additional comparison metrics), the main outcome of this paper is the fact that the combination of cortical wavelet expansions for the brain activity with the sparse Bayesian learning framework, provides relevant and robust source estimates in the case of extended activity, which do not rely on difficult parameter choices, and outperform significantly concurrent approaches of the same family (variational and Bayesian, in distributed source models).

The paper is organized as follows. After the current introduction, Section II presents reference solution that will be used in this paper. The construction of wavelets and the corresponding solvers are described in Section III, and evaluation metrics are discussed in Section IV. Section V presents the simulation protocol, and discusses and compares results obtained using various solvers. Finally, Section VI is devoted to discussion and conclusions. Preliminary results of this work were presented in the conference paper [7].

S. Mokhtari and B. Torrésani are with Aix Marseille Univ, CNRS, I2M, Marseille, France. email: samy.mokhtari@univ-amu.fr, bruno.torresani@univ-amu.fr

J.M. Badier and C. Bénar are with Aix Marseille Univ, Inserm, INS, Institut de Neurosciences des Systèmes, Marseille, France. email: jean-michel.badier@univ-amu.fr, christian.benar@univ-amu.fr

This work has benefited from support from the French government through the BMWs project (ANR-20-CE45-0018), and the *Institut Convergence ILCB* (ANR-16-CONV-0002).
# Combining spatial wavelets and sparse Bayesian learning for extended brain sources reconstruction

Samy Mokhtari, Jean-Michel Badier, Christian G. Bénar, Bruno Torrésani

**Abstract—Objective:** The accurate reconstruction of extended cortical activity from M/EEG data is a difficult, ill-conditioned problem. This work proposes to model distributed sources through spectral graph wavelets on the cortical surface, and addresses resulting numerical optimization problems. The objective is accurate localization, especially for extended sources, together with quantitatively relevant amplitude and time course.

**Approach:** Unknown sources are expanded on a system of spectral graph wavelets (SGW) defined on the cortical surface. Unknown wavelet coefficients are estimated using either variational or Bayesian formulations, involving priors that favor extended sources through sparsity in the wavelet domain: sparsity-inducing regularization, or sparse Bayesian learning (SBL). These approaches are tested and compared with concurrent approaches on real (open-access) data and numerical simulations. The quality of reconstructions is assessed using complementary metrics.

**Results:** SGW-based approaches are able to identify accurately extended sources. The combination with SBL is particularly attractive, as it doesn't involve hyperparameter tuning and automatically adapts to the signal to noise ratio (SNR). It yields accurate and robust results with respect to all considered metrics, and performs remarkably well in terms of depth bias.

**Conclusion:** This paper demonstrates the usefulness of spectral graph cortical wavelets for reconstructing cortical activity from M/EEG data, especially when coupling spatial wavelets with SBL.

**Significance:** Being able to identify localization, depth, amplitude and time course of brain activity from M/EEG data is important in clinical applications such as epilepsy, as it can improve the detection of potential sources of seizures.

**Index Terms**—M/EEG inverse problem, extended sources, spectral graph wavelets, Sparse Bayesian Learning.

## I. INTRODUCTION

Magneto- and Electro-encephalography (M/EEG) are the main (non-invasive) imagery modalities which allow, thanks to their high temporal resolution, physicians and clinicians to detect and localize cortical activities. From a given set of sensors and measurements, the aim is to reconstruct the sources activity. This is done by solving the M/EEG inverse problem [1], [2], [3], [4]. To address this challenge, multiple methods have been proposed, generally categorized as parametric versus distributed methods. While the former models the source activities with a few equivalent current dipoles of unknown locations and is therefore well adapted for

identifying focal, localized activities, the latter usually models the source space with a fine spatial discretization of the cortical surface, with ~ $10^4$ unknown cortical sources with known and fixed locations. While leading to very under-determined inverse problems, distributed methods are expected to be better suited to extended cortical activities. They lead to an optimization problem on a suitable objective function, which is solved using variational methods (e.g. Tikhonov-like regularizations). This optimization problem may also be interpreted, with a Bayesian point of view, as a Maximum-A-Posteriori estimate on the cortical sources distributions, given a specified prior. In both cases, the precise estimation of the spatial extent of the activity (a major feature for distributed methods) remains difficult. This paper introduces and studies variational and Bayesian approaches whose regularizations/priors rely on a spatial transformation, and aim at promoting extended solutions compressed in a relevant transformed domain. Indeed, as highlighted in [5]-[6], enforcing sparsity directly in the source domain may prove inadequate for extended sources. In [6], the spatial transformation relied on the difference of current density between neighboring vertices (namely VB-SCCD), while in [5], the transformation was implemented through a spatial wavelet transform based on a multi-resolution model of the brain cortical surface (with $l^1$ penalization). Following this idea of sparsity in a relevant transformed domain, we focus here on an overcomplete spatial wavelet transform (i.e. a wavelet frame, which differs from the multi-resolution approach), and combine it with dedicated sparsity-promoting regularizations/priors. The motivation behind a spatial wavelet transform originates from image processing, where sparsity in the wavelet domain is widely recognized as a relevant model for natural images, which often feature smooth regions separated by edges [7]. Overcompleteness means that the number of wavelets is larger than the dimension of the space, which often results in increased sparsity [8]. Variational approaches often suffer from the difficulty of choosing appropriate values for regularization parameters, and we show that the integration of these overcomplete wavelets into the Sparse Bayesian learning (SBL) paradigm yields more robust solutions. As introduced above, wavelet-based solutions (by which we mean solutions obtained with a penalization/prior defined in a wavelet domain) had already been studied in the framework of variational approaches, for instance in [9] where spherical (multi-resolution) wavelet bases were used, or in [10], [5], [11], with face-based multi-resolution wavelets. Here, the spatial wavelet transform relies on the Spectral Graph Wavelets (SGW) framework [12], which only requires a graph adjacency matrix, and yields expansions

S. Mokhtari and B. Torrésani are with Aix Marseille Univ, CNRS, I2M, Marseille, France. email: firstname.name@univ-amu.fr

J.M. Badier and C. Bénar are with Aix Marseille Univ, Inserm, INS, Institut de Neurosciences des Systèmes, Marseille, France. email: firstname.name@univ-amu.fr

This work has benefited from support from the French government through the BMWs project (ANR-20-CE45-0018), and grants (ANR-22-CE19-0012-01, ANR-11-INBS-0006) to France Life Imaging network, and (ANR-16-CONV-0002) to Institute of Language, Communication, and the Brain.
## II. STATE OF THE ART AND REFERENCE SOLUTIONS

Suppose M/EEG data are measured at $J_0$ sensors during $L$ time steps, and the cortical surface is discretized on $N$ source points. The M/EEG inverse problem starts with the observation equation connecting cortical sources to measurements :

$$Z_0 = G_0S + B_0, \quad (1)$$

where $Z_0 \in \mathbb{R}^{J_0 \times L}$ denotes space-time measurements, $S \in \mathbb{R}^{N \times L}$ the unknown sources activity, $B_0 \in \mathbb{R}^{J_0 \times L}$ the space-time matrix of both observation noise and background brain activity (zero-mean, with known covariance matrix $\Sigma_{B_0}$), and $G_0 \in \mathbb{R}^{J_0 \times N}$ is the leadfield matrix which encodes electromagnetic propagation from cortex to sensors.

Given $G_0$ and $\Sigma_{B_0}$ (and some additional assumptions), M/EEG inverse problem aims to estimate sources $S$ from observations $Z_0$. Before explicitly addressing observation equation (1), data are often pre-processed to correct for noise correlations between sensors; this can be done with a whitening and dimension reduction of the noise covariance matrix. The resulting matrix $B = \Upsilon B_0 \in \mathbb{R}^{J \times L}$ (where $\Upsilon \in \mathbb{R}^{J \times J_0}$ denotes the whitener) is thus white with zero mean and unit variance. For the sake of simplicity, we will hereafter consider whitened noise $B$, leadfield $G$ and data $Z$ matrices.

Several approaches are routinely used in literature to tackle observation equation (1), which include among others dipole localization [8], [9] or various forms of beamformers [10], [11]. This article will be mainly concerned with variational and Bayesian approaches for distributed source models (see [1]). Most variational methods boil down to the following generic optimization problem (2), which involves a quadratic data fidelity term and a penalty term $f(S)$:

$$S_* = \arg \min_{S \in \mathbb{R}^{N \times L}} \left[ \frac{1}{2} \|Z - GS\|_F^2 + f(S) \right], \quad (2)$$

where $\| \cdot \|_F$ stands for the Frobenius norm, i.e. the square root of the sum of squared matrix coefficients.

A wide variety of penalizations/regularizations can be found in literature. While the current paper does not aim to provide a comprehensive review, we nevertheless recall the penalizations that will be of most interest for the present study.

### A. Closed-form solutions with quadratic regularizations

The most commonly used penalization is quadratic:

$$f(S) = \|DS\|_F^2, \quad (3)$$

for some symmetric, generally diagonal, invertible matrix $D$. With the above penalization, the solution of the optimization problem (2) admits a closed-form expression:

$$S_* = D^{-2}G^T (I_J + GD^{-2}G^T)^{-1} Z. \quad (4)$$

Special cases of (3) include the MNE (minimal norm estimate) solution [12] ($D = \sqrt{\lambda}|I_N$), and the wMNE (weighted minimal norm estimate) solution [13], [14], where $D$ is a weight matrix intended to correct MNE depth bias. In e-LORETA [15], an iterative algorithm yields a data driven estimate for $D$ and therefore a corresponding source estimate.

While closed form solutions are very convenient, MNE and wMNE are however known to overestimate the spatial extent of sources $S$, and underestimate corresponding amplitudes. In addition, the parameter $\lambda$, linked to the source spatial covariance matrix $\Gamma = D^{-2}$, is an extra parameter whose fine tuning is not obvious. In many situations, the interest is more focused on the spatial extent, namely source localization. A number of post-processings are used to yield more accurate support estimates. dSPM [16] and s-LORETA [17] are among the most popular.

### B. Sparse solutions through non-differentiable penalizations

To reduce the blurring effect, other penalizations $f$ have been proposed in literature, often promoting some form of sparsity. There is now a vast literature on the subject, we refer to [18] for a review and a comparison of some recent contributions in the context of M/EEG inverse problem. In these approaches, no closed form solution exists as the penalization is no longer differentiable, and the optimization problem (2) has to be solved numerically. One may first think of the MCE (minimal current estimate) solution [19], [20], where $f(S) = \lambda \|S\|_1$ tends to yield very "spiky" source estimates, with quite irregular time courses. One may also cite the VB-SCCD approach [21], where $f(S) = \lambda \|S\|_{TV} = \lambda \|\nabla S\|_1$. The latter is based on total variation (TV), where $\nabla$ is a discrete gradient operator which has to be properly defined on the surface mesh. TV-based solutions tend to enforce sparsity in a transformed domain, here the gradient. A slightly modified and sparser version, called sVB-SCCD [21], [22] is based upon a combined $\ell^1$-TV penalization $f(S) = \lambda (\alpha \|S\|_1 + \|S\|_{TV})$. Hereafter, we shall rather consider the following form:

$$f(S) = \lambda (\gamma \|S\|_1 + (1 - \gamma) \|S\|_{TV}), \quad \gamma \in [0, 1] \quad (5)$$

In [22], $\alpha \in [0.01, 1]$, and $\alpha = 0.67$ is advocated as a reasonable choice for the tested scenarios. For both MCE and sVB-SCCD, numerical algorithms have to be carefully chosen, [22] advocates for ADMM. Regularization (5) allows for an easier interpretation of the hyperparameters, where $\lambda$ plays the same role as for the MCE solver, while $\gamma$ drives the convex combination of $\ell^1$ and TV-penalizations. Generalized TV penalization have been studied in [23].

Lastly, mixed norms (MxNE) combining sparsity in space with persistence in time were proposed and studied in [24], [25]. The regularization term is of the form $f(S) = \lambda \|S\|_{21} = \lambda \sum_k \|S[:k]\|_2$. Mixed norms were also used in combination with time-frequency dictionaries [26], and combined with TV spatial regularization [22].

### C. Bayesian point of view

A major weakness of the above mentioned approaches is the difficulty of tuning the hyperparameters involved in the formulations (e.g. regularization parameter $\lambda$). In simple situations, Bayesian interpretations yield a relationship between these parameters and quantities such as signal to noise ratio (SNR), which can be guessed or measured from sensor signals. For instance, wMNE may be interpreted in terms of a Maximum A Posteriori (MAP) estimation, modeling columns of $B$ (baseline) as iid samples from a normal distribution $N(0, I)$,
on an overcomplete spatial wavelet family. We investigate corresponding variational and SBL solutions. These are compared with real open-access data (auditory evoked potentials, median nerve stimulation), and reference simulated solutions using extensive numerical simulations designed to mimic real evoked potentials. Performance assessment is done using a set of metrics, which include standard ones (implemented in specialized computer packages such as MNE-Python [13], which we therefore use) as well as new ones which we found convenient in the context of sparse solutions.

Besides methodological contributions (e.g. the use of SGW frames and additional comparison metrics), the main outcome of this paper is to show that the combination of such overcomplete wavelets with the SBL framework provides relevant and robust source estimates for extended activities. The approach does not require difficult parameter choices, adapts automatically to the SNR, and outperforms significantly concurrent distributed sources reconstruction methods.

The paper is organized as follows. Section II presents the reference M/EEG inverse problem solvers used in this paper. The construction of wavelets and corresponding solvers are described in Section III, and evaluation metrics are discussed in Section IV. Section V presents the tested real open-access data and the simulation protocol, and discusses results obtained using various solvers. Finally, Section VI is devoted to conclusions. Preliminary results of this work were presented in the conference paper [14].

## II. STATE OF THE ART, DISTRIBUTED SOLUTIONS

Consider M/EEG data measured at $J_0$ sensors and $L$ time samples, with brain cortical surface discretized on $N$ points. The observation equation relating cortical sources $\mathbf{S} \in \mathbb{R}^{N \times L}$ (here with constrained orientations) and data $\mathbf{Z}_0 \in \mathbb{R}^{J_0 \times L}$ is:

$$ \mathbf{Z}_0 = \mathbf{G}_0\mathbf{S} + \mathbf{B}_0, \quad (1) $$

given the leadfield matrix $\mathbf{G}_0 \in \mathbb{R}^{J_0 \times N}$, and unknown baseline $\mathbf{B}_0 \in \mathbb{R}^{J_0 \times L}$ (zero-mean, with known covariance matrix $\Sigma_{\mathbf{B}_0}$).

Data are often pre-processed to correct for noise correlations between sensors; this can be done by whitening and dimension reduction of the noise covariance matrix. The resulting matrix $\mathbf{B} = \mathbf{Y}\mathbf{B}_0 \in \mathbb{R}^{J \times L}$ ($\mathbf{Y} \in \mathbb{R}^{J \times J_0}$ denotes the whitener) is thus spatially white. For simplicity, we will hereafter consider whitened noise $\mathbf{B}$, leadfield $\mathbf{G}$ and data $\mathbf{Z}$ matrices.

Approaches routinely used in literature include among others dipole localization [15], [16] or various forms of beamformers [17], [18], [19]. This article is mainly concerned with variational and Bayesian approaches for distributed source models (see [1]). By variational methods we mean solving a suitable optimization problem, involving a quadratic data fidelity term and a regularization term $f(\mathbf{S})$:

$$ \mathbf{S}_* = \arg \min_{\mathbf{S} \in \mathbb{R}^{N \times L}} \left[ \frac{1}{2} \|\mathbf{Z} - \mathbf{G}\mathbf{S}\|_F^2 + f(\mathbf{S}) \right], \quad (2) $$

where $\| \cdot \|_F$ stands for the Frobenius norm.

A wide variety of regularizations $f(\mathbf{S})$ can be found in literature. They embed one or several parameters, which drive the compromise between data fidelity and regularizations.

While the current paper does not aim to provide a self-contained review, we nevertheless recall regularizations that will be of most interest for the present study.

### A. Regularizations in the source domain

The most commonly used regularizations in the source domain are quadratic: $f(\mathbf{S}) = \|\mathbf{D}\mathbf{S}\|_F^2$ for some symmetric, generally diagonal, invertible matrix $\mathbf{D}$. These allow a closed-form solution (see for example [1]). Special cases include MNE [20] ($\mathbf{D} = \sqrt{\lambda}|I_N$), and wMNE [21], [22], where $\mathbf{D}$ aims to correct MNE depth bias. In e-LORETA [23], an iterative algorithm yields a data driven estimate for $\mathbf{D}$ and therefore a corresponding source estimate. MNE and wMNE are known to overestimate the spatial extent of sources $\mathbf{S}$, and underestimate corresponding amplitudes. Also, fine tuning the parameter $\lambda$, or the matrix $\mathbf{D}$ is often difficult.

L1-norm based regularizations, i.e. $f(\mathbf{S}) = \lambda \|\mathbf{S}\|_1$, enforce sparsity and allow reducing the blurring effect. There, no closed form solution exists, and the problem (2) has to be solved numerically. MCE [24], [25] tends to yield very "spiky" source estimates, with quite irregular time courses, and is therefore more adapted to focal sources. Once again, parameter tuning remains complex. We refer to [26] for a review of recent contributions, including variants such as MxNE that combine spatial sparsity with time persistence.

### B. Regularizations in a transformed domain

As highlighted in [27], straightforward sparsity in the source domain may be counterproductive for recovering extended sources, whose relevance has been shown in [28], [29]. To bypass this drawback, one may look for sparsity in a transformed domain: corresponding solutions are sparse in the transformed domain, but not in the source domain any more.

For example, VB-SCCD [6] minimizes the total variation (TV) $\lambda \|\mathbf{S}\|_{\text{TV}} = \lambda \|\nabla \mathbf{S}\|_1$, using a discrete gradient operator $\nabla$ defined on the surface mesh. TV-based solutions tend to enforce sparsity in the gradient domain, which favours solutions that are constant over patches in the source domain. We use here a slightly modified and sparser version, called sVB-SCCD [6], [30], based upon a combined $l^1$-TV regularization

$$ f(\mathbf{S}) = \lambda (\gamma \|\mathbf{S}\|_1 + (1 - \gamma) \|\mathbf{S}\|_{\text{TV}}), \quad \gamma \in [0, 1] \quad (3) $$

For both MCE and sVB-SCCD, numerical algorithms must be carefully chosen, and parameter tuning remains complex. In [6], VB-SCCD was only compared to wMNE and cLORETA, on simulated data, with one evaluation metrics: ROC-AUC, corrected to limit bias inherent to unbalanced sets. The issue of amplitude thresholding inherent to ROC curves was outlined as a limitation, as well as the underestimation (respectively overestimation) of sources amplitudes (respectively extents).

We especially highlight compressible spatial transforms based on the Laplacian matrix or spatial (spherical) wavelets (with $l^1$ penalization both in the original and transformed domains, see [9]). The guiding idea behind such a compressible representation was to be able to address both focal and extended sources. This is done by tuning a regularization parameter in the cost function, balancing between penalizations in the transformed and original domains (thus adding a parameter
and columns of $\mathbf{S}$ as iid samples from a normal distribution $N(0, \mathbf{\Gamma})$, with $\mathbf{\Gamma} = \lambda^{-1}|_N$ in the case of MNE. The posterior distribution of sources $p(s(l)|z(l))$ is multivariate normal, with covariance matrix

$$ \Sigma_S = \mathbf{\Gamma} - \mathbf{\Gamma} \mathbf{G}^T (\Sigma_B + \mathbf{G} \mathbf{\Gamma} \mathbf{G}^T)^{-1} \mathbf{G} \mathbf{\Gamma}. \quad (6) $$

In the MCE case, modeling columns of $\mathbf{S}$ as iid samples from a product of 1D Laplacian distributions with parameter $\lambda^{-1}$, the sources covariance matrix is given by $\mathbf{\Gamma} = 2\lambda^{-2}|_N$. For both MNE and MCE solutions, heuristics can be derived to connect the regularization parameter $\lambda$ to the square signal to noise ratio $\rho^2 = \text{Tr}(\Sigma_{Z_0}) / \text{Tr}(\Sigma_{B_0})$ (signal and noise levels in the sensor domain are defined by the traces of the covariance matrices $\Sigma_{Z_0}$ and $\Sigma_{B_0}$):

$$ \lambda_{\text{MNE}} = \frac{\|\mathbf{G}_0\|_F^2}{(\rho^2 - 1)\text{Tr}(\Sigma_{B_0})}; \quad \lambda_{\text{MCE}} = \sqrt{2\lambda_{\text{MNE}}} \quad (7) $$

Such heuristics are more difficult to obtain for penalization terms $f(\mathbf{S})$ that involve TV or mixed norms, and even worse for combinations of norms that involve several parameters.

**Remark 1:** The above heuristics differ slightly from expressions routinely used in reference softwares, which perform similar calculations from whitened observations $\mathbf{Z}$ instead of original ones $\mathbf{Z}_0$ as we do. We believe the expressions we use are more relevant when interpreted in a Bayesian context.

### D. Sparse Bayesian Learning

Beside MAP approaches, Sparse Bayesian Learning (SBL, see for example [27], [28]) allows estimating simultaneously the sources amplitudes and their covariance matrices, assumed to be diagonal $\mathbf{\Gamma} = \text{diag}(\gamma_1, \dots, \gamma_N)$. SBL-based inverse problems have recently been put into a unifying context [29], for which most algorithms fall into the class of Majorization-Minimization algorithms (MM, see [30] for a review). They essentially boil down to solving the following non-convex optimization problem, with respect to the source variances $\underline{\gamma} = (\gamma_1, \dots, \gamma_N) \in (\mathbb{R}_+)^N$

$$ \underline{\gamma}_* = \arg\min_{\underline{\gamma} \in (\mathbb{R}_+)^N} \left[ \frac{1}{L} \text{Tr}(\mathbf{Z}^T \mathbf{\Sigma}_{\mathbf{Z}}^{-1} \mathbf{Z}) + \ln \det(\mathbf{\Sigma}_{\mathbf{Z}}) \right], \quad (8) $$

where $\mathbf{\Sigma}_{\mathbf{Z}} = (|\mathbf{J}| + \mathbf{G} \mathbf{\Gamma} \mathbf{G}^T)$ is the posterior covariance matrix. Once $\mathbf{\Gamma}$ has been estimated, its value can be plugged in the (wMNE closed-form) expression of $\mathbf{S}$.

For the record, let us mention the alternative Bayesian approach, based on the MEM (maximum entropy on the mean) principle [31], [32], [33]. Remarkably, the MEM approach leads to a smooth optimization problem in the domain of sensors (instead of sources), hence of much lower complexity, and for which simple and faster differentiable optimization tools can be used. The wavelet techniques to be explained below would easily be included in this approach.

### III. USE OF A FAMILY OF CORTICAL WAVELETS

Approaches that exploit sparsity in the spatial domain, such as dipole models or MCE/MxNE are well adapted to identify focal activity. However, when the activity of interest has a significant spatial extent, they produce estimates that are often

too sparse. In such situations, it is wise to search for solutions that are sparse in a transformed domain. As an example, VB-SCCD [18] optimizes the sparsity of a spatial gradient of the solution. Similar approaches based on Laplacian and spherical wavelet transform have been considered in [4], where authors report improved results for localizing extended activity.

#### A. Problem formulation

Bases or redundant systems of waveforms have often been used to describe time dependence and perform dimension reduction in the time domain. Examples include wavelets [32], [33], time-frequency systems, or waveforms estimated from data. Such approaches do not seem to have been considered for the spatial dependence of sources (except [4], and approaches based upon cortex parcellation [32]), due to the lack of such simple waveform families. In this section, we introduce and use spatial wavelets for the M/EEG inverse problem.

We define wavelets as functions $w_k \in \mathbb{R}^N$ on the cortex, labelled by a vertex and a scale index. Concatenating these vectors into a matrix $\mathbf{W} = (w_1 \dots w_{N_w})^T \in \mathbb{R}^{N_w \times N}$, the family of wavelets forms a frame for functions on the surface if $\mathbf{W}$ is left invertible (in other words the family of wavelets is a complete set in $\mathbb{R}^N$). The construction of such a frame will be detailed in the next section. $\mathbf{W}$ is the analysis operator of the wavelet frame. A space-time matrix $\mathbf{S} \in \mathbb{R}^{N \times L}$ can be written as $\mathbf{S} = \mathbf{W}^T \mathbf{A}$, in terms of a wavelet-time coefficient matrix $\mathbf{A} \in \mathbb{R}^{N_w \times L}$. Putting the focus on the matrix $\mathbf{A}$ is called a synthesis-based approach, as opposed to the analysis-based approach which rests on the matrix $\mathbf{W}\mathbf{S} \in \mathbb{R}^{N_w \times L}$. In this synthesis-based setting, the generic inverse problem (2) may be written as

$$ \mathbf{A}_* = \arg\min_{\mathbf{A} \in \mathbb{R}^{N_w \times L}} \left[ \frac{1}{2} \|\mathbf{Z} - \mathbf{G}_W \mathbf{A}\|_F^2 + g(\mathbf{A}) \right] \quad (9) $$

with $\mathbf{G}_W = \mathbf{G} \mathbf{W}^T$ the transformed leadfield matrix. This formulation is very close to the usual one, the difference being that the optimization is done in the synthesis coefficients domain. Notice that when solving (9), the waveform system is used only twice, for the computation of the transformed leadfield matrix $\mathbf{G}_W$ and for the synthesis of the solution $\mathbf{S}_*$. All reference solvers for the M/EEG inverse problem can be transposed mutatis mutandis to this new setting. We stress that the dimension of the transformed space $(N_w \times L)$ may be significantly larger than the one of the original space $(N \times L)$, as the number $N_w$ of waveform atoms used to decompose the solution may be, in the case of a redundant wavelet family, several times larger (depending on the number of wavelet scales considered) than the number of sources $N$.

#### B. Signals, Fourier transform and wavelets on graphs

We now detail the construction of the wavelet frame $\mathbf{W}$ for the graph describing cortical sources. In the present work, we decided to follow the spectral graph wavelets (SGW) framework developed in [5].

Consider the cortical mesh graph $(\mathcal{V}, \mathcal{E}, \omega)$ (with vertex set $\mathcal{V}$, edge set $\mathcal{E}$ and weight function $\omega$, the latter allowing to define the graph adjacency matrix $\mathcal{A}$). The graph Laplacian is given by $\mathcal{L} = \mathcal{D} - \mathcal{A}$, where $\mathcal{D}$ is the diagonal degree
tuning issue, as for VB-SCCD). Suggested methods showed improved results compared to MNE and the sparse source imaging method of [31], both for simulated and somatosensory data. Evaluation metrics such as Spatial Dispersion, Dipole Localization Error and ROC curves were considered, once again with the issue of amplitude thresholding.

In [10], [5] face-based spatial wavelets were introduced (with $l^1$ penalization on wavelets), then extended in [11] ($l^1$ penalization on both variation, see [6], and wavelets). [5] reports better performances of face-based wavelets for sources reconstruction, compared to vertex-based wavelets, MNE and cLORETA, with respect to classical metrics such as ROC curves and AUC, Spatial Dispersion, and Distance of Localization Error (with, again, the issue of amplitude thresholding). This literature also suggests that face-based wavelets achieve better compression performance. Furthermore, while still challenging, adaptive estimation of sources extents could be tackled by wavelets compression level. In [11], penalization on the variation term is designed to mitigate the impact of the wavelets compression level.

Whatever the wavelet design, such wavelet transforms face a common weakness and issue : in addition to the compression level, remains the challenge of tuning the regularization parameter (which for instance may drive the tradeoff between penalization in the original and transformed domains), which is acknowledged as a main limitation in [11].

### C. Bayesian point of view

In simple situations, Bayesian interpretations yield a relationship between these parameters and quantities such as signal to noise ratio (SNR), which can be guessed or measured from sensor signals. For instance, wMNE may be interpreted in terms of a Maximum A Posteriori (MAP, also called Type I) estimation, modeling columns of B (resp. S) as iid samples from a normal distribution $N(0, I)$ (resp. a normal distribution $N(0, \Gamma)$). Similarly, MCE is obtained by modeling columns of S as iid samples from a product of 1D Laplacian distributions with parameter $\lambda^{-1}$ (i.e. sources covariance matrix $\Gamma = 2\lambda^{-2}I_N$). For both MNE and MCE, heuristics can be derived to connect the parameter $\lambda$ to the squared signal to noise ratio $SNR^2 = \text{Tr}(\Sigma_{Z_0})/\text{Tr}(\Sigma_{B_0})$ (signal and noise levels in the sensor domain are defined by the traces of the covariance matrices $\Sigma_{Z_0}$ and $\Sigma_{B_0}$):

$$ \lambda_{\text{MNE}} = \frac{\|\mathbf{G}_0\|_F^2}{(\text{SNR}^2 - 1)\text{Tr}(\Sigma_{B_0})}; \quad \lambda_{\text{MCE}} = \sqrt{2\lambda_{\text{MNE}}} \quad (4) $$

Such heuristics are more difficult to obtain for penalization terms $f(S)$ that involve TV or mixed norms, and even worse for combinations of norms that involve several parameters.

**Remark 1:** The above heuristics differs slightly from expressions routinely used in reference softwares, which perform similar calculations on whitened observations Z instead of original measurements, $Z_0$, as done here. We believe the expressions we use are more relevant when interpreted in a Bayesian context.

### D. Sparse Bayesian Learning

In MAP (Type I) approaches, the prior probability distribution of sources is fixed. Type II approaches extend the

problem to parametric families of priors, that depend on hyper-parameters to be estimated in a hierarchical Bayes framework [32]. In the context of interest here, the hyper-parameters $\gamma \in \mathbb{R}_+^N$ are the diagonal elements of the prior source covariance matrix $\Gamma = \text{diag}(\gamma)$, and are estimated by maximizing the conditional density $p(Z|\gamma) = \int p(Z|S, \gamma)p(S|\gamma) dS$. Once $\Gamma$ has been estimated, its value can be plugged into the MAP (Type I) estimate of S.

Taking $p(S|\gamma) = \prod_{l=0}^{L-1} N(0, \Gamma)$, with L the number of time samples, the MAP estimate is simply the wMNE solution based upon $\Gamma$. The corresponding type II approach goes under the name of Sparse Bayesian Learning [32] (SBL), and is conveniently reformulated in [33] as the non-convex optimization problem

$$ \underline{\gamma}^* = \arg \min_{\underline{\gamma} \in \mathbb{R}_+^N} \left[ \text{Tr} (C_z \Sigma_z (\underline{\gamma})^{-1}) + \ln \det (\Sigma_z(\underline{\gamma})) \right] \quad (5) $$

where $C_z$ is an estimate of the data covariance matrix, and $\Sigma_z(\underline{\gamma})$ is the posterior covariance matrix of observations

$$ \Sigma_z(\underline{\gamma}) = I_J + G \Gamma G^T . \quad (6) $$

The problem (5) is solved using numerical algorithms, which are described in details in Section 2 of [34].

SBL-based inverse problems have recently been put into a unifying context in [34] (see also [35]), which gives a review of several algorithms that can be used to solve (5). Most algorithms fall into the class of Majorization-Minimization family (MM, see [36] for a review). In the next Section, SBL is formulated in a transformed domain, in the same spirit as [37].

## III. USE OF A FAMILY OF CORTICAL WAVELETS

In the current work, we suggest to couple another construction of spatial wavelets (different from the one in [9], [10], [5]), namely Spectral Graph Wavelets (SGW), with the previously described Sparse Bayesian Learning (SBL) framework. The first foundation, SGW, which only requires a graph adjacency matrix, allows for a spatial wavelets framework compatible with any spatial discretization. Furthermore, spatial wavelets built within SGW framework result in overcomplete wavelet dictionaries, property which may be of interest for recovering extended sources and their redundant spatial features.

As for the second foundation, SBL, it allows tackling the regularization hyperparameter conundrum, by "learning", through a non-convex optimization problem and all the time samples at hand, N different parameters (most of which will be zero, thus enforcing sparsity) for the N unknown sources.

### A. Problem formulation

The spectral graph wavelets (SGW) framework of [12] defines wavelets as functions on the cortex (i.e. N-dimensional vectors), labelled by a vertex and a scale index. Concatenating these vectors into a $N_w \times N$ matrix W, the family of $N_w$ wavelets forms a frame of $\mathbb{R}^N$ if W is left invertible (in other words the family of wavelets is a complete set in $\mathbb{R}^N$). The construction of such a frame will be detailed in the next section. W is the analysis operator of the wavelet frame. A space-time matrix $S \in \mathbb{R}^{N \times L}$ can be written as $S = W^T A$, in terms of a wavelet-time coefficient matrix $A \in \mathbb{R}^{N_w \times L}$.
matrix, whose diagonal elements $d_i$ equal the number of edges connected to edge $i$. $\mathcal{L}$ is a real symmetric matrix (for undirected graphs), and therefore possesses a complete set of orthonormal eigenvectors $\{\chi_1, \dots, \chi_N\}$, associated to real, non-negative and increasing eigenvalues $\{\lambda_1, \dots, \lambda_N\}$. The minimal eigenvalue $\lambda_1$ equals zero, with a multiplicity equal to the number of connected components of the graph. The linear transform $\mathcal{F}: X \in \mathbb{R}^N \mapsto \{\mathcal{F}[X](n) = \langle \chi_n, X \rangle\}_{1 \le n \le N}$ may be interpreted as the graph analogue of the Fourier transform. The orthonormality of the basis $\{\chi_1, \dots, \chi_N\}$ gives the Parseval formula and the inverse transform: for any $X \in \mathbb{R}^N$, $\|\mathcal{F}[X]\| = \|X\|$ and $X = \sum_{n=1}^N \mathcal{F}[X](n)\chi_n$.

Following [5], a graph analogue of the wavelet transform can be derived from the above graph Fourier transform. Given a dilation parameter $s > 0$ and a vertex $n$, the graph wavelet transform of a signal $X \in \mathbb{R}^N$ may be defined as follows:

$$ \mathcal{W}[X](s, n) = \langle \Psi_{s,n}, X \rangle, \quad (10) $$

with $\Psi_{s,n}$ the wavelet atom associated to the vertex $n$ and dilation parameter $s$. In [5], these atoms $\Psi_{s,n}$ are generated in the spectral domain, as the graph Fourier transforms of Dirac functions centered on each vertex, multiplied by a band-pass positive kernel $g$: for any $l \in \{1, \dots, N\}$:

$$ \mathcal{F}[\Psi_{s,n}](l) = g(s\lambda_l) \langle \chi_l, \delta_n \rangle \quad (11) $$

As outlined in [5], such a construction may lead to an ill-conditioned matrix for the synthesis operator (frame vectors associated to small scales $s$ may become numerically parallel). To bypass this difficulty, [5] suggests to encapsulate the low-frequency component of signals through an ad hoc "scaling function" $\Phi_n$ (labelled by a vertex index), defined independently, in the spectral domain, via a low-pass kernel $h$:

$$ \mathcal{F}[\Phi_n](l) = h(\lambda_l) \langle \chi_l, \delta_n \rangle \quad (12) $$

For the specific design of both low-pass $h$ and band-pass $g$ kernels, we refer the reader to [5].

### C. Wavelet-sparse solutions

We now consider the inverse problem formulated in the wavelet domain, using the synthesis point of view. This requires computing the transform-domain leadfield matrix $G_W = GW^T$. This is the most costly part of the approach.

1) **Variational formulations:** Denote by $A = [A_\Psi, A_\Phi]$ the vector of multiscale (wavelet+scaling) synthesis coefficients:

$$ S = W^T A = \Psi^T A_\Psi + \Phi^T A_\Phi . $$

In addition to classical regularizations in the sources or gradient domains, the wavelet domain optimization problem (9) will hereafter be solved using two different penalizations:

$$ g_{\text{MNE}}(\mathbf{A}) = \|\mathbf{D}\mathbf{A}\|_{F}^{2}, \quad g_{\text{MCE}}(\mathbf{A}) = \lambda_{\text{MCE}} \|\mathbf{A}\|_{1} \quad (13) $$

The first one will be coupled with a wavelet-SBL approach to estimate the weight matrix $\mathbf{D}$, while the second aims at promoting sparsity in the wavelet coefficient domain.

2) **Sparse Bayesian Learning approaches:** The wavelet version of SBL is also a simple transposition of classical SBL approach, and yields the non convex problem

$$ \frac{\gamma_A}{\gamma_A \in (\mathbb{R}_+)^N} = \arg \min_{\gamma_A \in (\mathbb{R}_+)^N} \left[ \text{Tr}(\mathbf{C}_Z \mathbf{\Sigma}_Z (\gamma_A)^{-1}) + \ln \det(\mathbf{\Sigma}_Z (\gamma_A)) \right] \quad (14) $$

where $\Gamma_A = \text{diag}(\gamma_A)$ is the prior covariance matrix of wavelet synthesis coefficients, $\mathbf{\Sigma}_Z(\gamma_A) = I_J + \mathbf{G}_W \Gamma_A \mathbf{G}_W^T$ is the posterior covariance matrix of observations, and $\mathbf{C}_Z = \frac{1}{L} \mathbf{Z}\mathbf{Z}^T$ is the sample covariance matrix. Once $\gamma_A^*$ has been obtained, estimated sources are obtained via

$$ \mathbf{A}_* = \Gamma_A^* \mathbf{G}_W^T [\mathbf{\Sigma}_Z(\gamma_A^*)]^{-1} \mathbf{Z}, \quad \mathbf{S}_* = \mathbf{W}^T \mathbf{A}_* \quad (15) $$

This approach will be called sgw-SBL or w-SBL.

## IV. ANALYSIS OF EVALUATION METRICS

To assess, qualitatively and quantitatively, the performances of any inverse problem model for the distributed sources framework, evaluation metrics are mandatory. Given a simulation case where the ground truth is known, one may try to encapsulate, with well-chosen metrics:

* purely-distance based information (spatial extent of true/estimated sources, spatial deviation or discrepancy between these two sets...);

* information related to the number of true/estimated active sources, the number of true positives/negatives and false positives/negatives;

* amplitude-based information;

* or a mix of these features.

Beside these various choices and associated metrics, one shall first state which physical value will be studied under these scopes. Indeed, one could study straightforwardly the sources time-dependent amplitudes, their time-dependent absolute or squared value, or rather a time-average of these amplitudes on a given time segment of interest. One may also consider either normalized amplitudes or not. In the following, to facilitate interpretation and comparisons among all considered inverse problem solvers, all chosen metrics/scopes will be applied on time-averaged values, as hereafter detailed: given a space-time matrix $\mathbf{S}_* \in \mathbb{R}^{N \times L}$ gathering the amplitude maps for all time snapshots, a time-average $\langle \mathbf{S}_*^2 \rangle$ is computed on a prescribed *time segment of interest* denoted by soi, thus leading to an amplitude map $\mathcal{A} \in \mathbb{R}^N$ describing the amount of energy contained in the *segment of interest*:

$$ \mathcal{A} = \sqrt{\langle \mathbf{S}_*^2 \rangle_{\text{soi}}} \quad (16) $$

In this section, we lay down and describe the several evaluation metrics that will be used to assess the models performances. To this end, let us first detail some notations.

**Notations:** from now on,

* $I$ denotes the index set of true active sources; its cardinality is denoted by $|I|$.

* $\hat{I}_\alpha$ denotes the index set of estimated active sources, defined with respect to a threshold parameter $\alpha > 0$;

* $|\tilde{s}_i| = \frac{1}{\sum_{k=1}^N |s_k|} |s_i|$ denotes the normalized (positive, sum to 1) sources amplitudes;

* $\| \cdot \|$ denotes the 3D Euclidean norm;

* $r_i$ denotes the 3D vector position of source $i$.
As opposed to the *analysis-based* approach which rests on the matrix $\mathbf{WS} \in \mathbb{R}^{N_w \times L}$, the *synthesis-based* setting focuses on $\mathbf{A}$ and the generic inverse problem (2) writes:

$$ \mathbf{A}_{*} = \arg \min_{\mathbf{A} \in \mathbb{R}^{N_w \times L}} \left[ \frac{1}{2} \|\mathbf{Z} - \mathbf{G}_W \mathbf{A}\|_F^2 + g(\mathbf{A}) \right] \quad (7) $$

with $\mathbf{G}_W = \mathbf{G}\mathbf{W}^T$ the transformed leadfield matrix. This formulation is very close to the usual one, the difference being that regularization, and then optimization are done on synthesis coefficients. Notice that when solving (7), the wavelet matrix is used only twice, for the computation of the transformed leadfield matrix $\mathbf{G}_W$ and for the synthesis of the solution $\mathbf{S}_*$.

All reference solvers for the M/EEG inverse problem can be transposed *mutatis mutandis* to this new setting. We stress that the dimension $N_w \times L$ of the transformed space (with $L$ the number of time samples) may be significantly larger than the one of the original space ($N \times L$), as the number $N_w$ of wavelets may be, in the case of a redundant wavelet family, several times larger than the number of sources $N$ ($N_w = (N_s + 1)N$, with $N_w$ the number of wavelets).

## B. Signals, Fourier transform and wavelets on graphs

We now detail the construction of the wavelet frame $\mathbf{W}$ for the graph associated with a triangulated cortical surface: graph vertices and edges are the vertices and edges of triangles. Generalizations include weighted graphs, with weights depending on edge lengths, or originating from connectivity data as in [38]. We stick here to the simplest situation, where all edges have unit weight.

Consider the cortical mesh graph $(\mathcal{V}, \mathcal{E})$, with vertex set $\mathcal{V}$ and edge set $\mathcal{E}$ and denote by $\mathbf{Ad}$ the graph adjacency matrix. The graph Laplacian is given by $\mathbf{L} = \mathbf{D} - \mathbf{Ad}$, where $\mathbf{D} = \text{diag}(d_1, \dots, d_N)$ is the degree matrix, whose diagonal element $d_i$ equals the number of edges connected to edge $i$. $\mathbf{L}$ is a real symmetric matrix, and therefore admits an eigenvalue decomposition $\mathbf{L} = \mathbf{F}\mathbf{A}\mathbf{F}^T$. $\mathbf{F}$ is the orthogonal matrix of Laplacian eigenvectors, and $\mathbf{A} = \text{diag}(\lambda_1, \dots, \lambda_N)$ encodes the real, non-negative and increasing eigenvalues $\{\lambda_1, \dots, \lambda_N\}$. The minimal eigenvalue $\lambda_1$ equals zero, with a multiplicity equal to the number of connected components of the graph. By analogy with the usual Laplacian on the cyclic graph $\mathbb{Z}_N$ (or other classical settings), $\mathbf{F}$ is called a Fourier matrix, its columns playing a similar role to that of sine waves in conventional Fourier analysis, but defined on the graph. The linear transform $X \in \mathbb{R}^N \to \mathbf{F}^T X$ is the graph analogue of the classical Fourier transform, and enjoys the same properties (norm conservation, and inversion). Eigenvalues $\lambda_k \in \mathbb{R}_+$ correspond to spatial frequencies.

Following [12], a graph analogue of the wavelet transform can be derived, using suitably rescaled kernel functions $g, h$ defined on $\mathbb{R}_+$, interpreted as frequency responses of band-pass and low-pass filters respectively. Given a dilation parameter $s > 0$, introduce the $N \times N$ matrices

$$ \mathbf{\Psi}_s = \mathbf{F} \text{diag}(g(s\lambda_1), \dots, g(s\lambda_N)) \mathbf{F}^T, \quad (8) $$

$$ \mathbf{\Phi}_s = \mathbf{F} \text{diag}(h(s\lambda_1), \dots, h(s\lambda_N)) \mathbf{F}^T. \quad (9) $$

$\mathbf{\Psi}_s$ is the wavelet transform at scale $s$, and $\mathbf{\Phi}_s$ is a scaling function transform at scale $s$, similar to what is done in

multiresolution wavelets. Wavelets and scaling functions at scale $s$ are the columns of $\mathbf{\Psi}_s$ and $\mathbf{\Phi}_s$ respectively.

Given scale values $s_1 > \dots > s_{N_s}$, let $\mathbf{\Psi}$ be the horizontal concatenation of wavelet matrices $\mathbf{\Psi}_{s_k}$ and $\mathbf{\Phi} = \mathbf{\Phi}_{s_1}$. The complete transform is given by the matrix

$$ \mathbf{W} = [\mathbf{\Phi} \ \mathbf{\Psi}]^T = [\mathbf{\Phi}_{s_1} \ \mathbf{\Psi}_{s_1} \ \dots \ \mathbf{\Psi}_{s_{N_s}}]^T. \quad (10) $$

For suitable choices of kernels $g, h$ and scales $s_k$, the transform is invertible (one says that the rows of $\mathbf{W}$ form a frame in $\mathbb{R}^N$): for every $X \in \mathbb{R}^N$, there exists a coefficient vector $A \in \mathbb{R}^{N_w}$, with $N_w = N(N_s + 1)$, such that $X = \mathbf{W}^T A$.

We refer to Equation (69) of Section 8 in [12] for the specific design of low-pass and band-pass kernels $h, g$, based upon piecewise polynomials. Examples are plotted in Fig. 1. Fig. 2 displays (false) color images of some wavelets and scaling functions on the cortical surface.

**Remark 2:** Spectral graph wavelets differ significantly from wavelet bases already used in M/EEG source localization. They do not form a basis of the source space, but an over-complete dictionary, and do not involve subsampling at each scale: to each vertex of the mesh graph is associated a family of $N_s$ wavelets and a scaling function. Redundancy induces an extra computational cost, but can be expected to yield sparser solutions: the more redundant the wavelet family, the more likely it contains waveforms that better fit the solution.

## C. Wavelet-sparse solutions

We now consider the inverse problem formulated in the wavelet domain, using the synthesis point of view. This requires computing the transform-domain leadfield matrix $\mathbf{G}_W = \mathbf{G}\mathbf{W}^T$. This is the most costly part of the approach.

1) **Variational formulations:** Denote by $\mathbf{A} = [\mathbf{A}_{\Phi}, \mathbf{A}_{\Psi}]$ the vector of multiscale (wavelet+scaling) synthesis coefficients:

$$ \mathbf{S} = \mathbf{W}^T \mathbf{A} = \mathbf{\Phi}^T \mathbf{A}_{\Phi} + \mathbf{\Psi}^T \mathbf{A}_{\Psi}. $$

In addition to classical regularizations in the sources or gradient domains, the wavelet domain optimization problem (7) will hereafter be solved using two different penalizations:

$$ g_{\text{MNE}}(\mathbf{A}) = \|\mathbf{D}\mathbf{A}\|_F^2, \quad g_{\text{MCE}}(\mathbf{A}) = \lambda_{\text{MCE}} \|\mathbf{A}\|_1 \quad (11) $$

The first one is used in the wavelet Sparse Bayesian Learning (w-SBL) approach that estimates the weight matrix $\mathbf{D}$, the second aims at promoting sparsity in the wavelet domain.

2) **Sparse Bayesian Learning approaches:** The wavelet version of SBL is also a simple transposition of classical SBL approach, and yields the non convex problem

$$ \gamma_A^* = \arg \min_{\gamma_A \in (\mathbb{R}_+)^{N_w}} \left[ \text{Tr}(\mathbf{C}_Z \mathbf{\Sigma}_Z (\gamma_A)^{-1}) + \ln \det(\mathbf{\Sigma}_Z(\gamma_A)) \right] \quad (12) $$

where $\Gamma_A = \text{diag}(\gamma_A)$ is the prior covariance matrix of wavelet synthesis coefficients, $\mathbf{\Sigma}_Z(\gamma_A) = I_J + \mathbf{G}_W \Gamma_A \mathbf{G}_W^T$ is the posterior covariance matrix of observations, and $\mathbf{C}_Z = \frac{1}{L} \mathbf{Z}\mathbf{Z}^T$ is the sample covariance matrix. Once $\gamma_A^*$ has been obtained, estimated sources are computed via

$$ \mathbf{A}_* = \Gamma_A^* \mathbf{G}_W^T [\mathbf{\Sigma}_Z(\gamma_A^*)]^{-1} \mathbf{Z}, \quad \mathbf{S}_* = \mathbf{W}^T \mathbf{A}_*. \quad (13) $$

This approach will be called sgw-SBL or w-SBL. The numerical computation of $\gamma_A^*$ in (12) is done using the CHAM-PAGNE algorithm, with McKay updates ([34], section 2.3.3).
A. Region Localization Error : purely distance-based metrics

We start with the classical Region Localization Error (RLE). Given a thresholding parameter $\alpha > 0$ (by default a percentage of the maximum value), and assuming $I \cap \hat{I}_{\alpha} \neq \emptyset$ and $\hat{I}_{\alpha} \neq I$, the RLE metrics is computed as follows :

$$ RLE(\alpha) = \frac{1}{2\#I} \sum_{k \in I \setminus (I \cap \hat{I}_{\alpha})} \min_{l \in \hat{I}_{\alpha}} \|r_k - r_l\| + \frac{1}{2\#\hat{I}_{\alpha}} \sum_{l \in \hat{I}_{\alpha} \setminus (I \cap \hat{I}_{\alpha})} \min_{k \in I} \|r_k - r_l\| \quad (17) $$

RLE combines equi-balanced information related to, on the one hand, the average smallest distance between the set of unexplained true sources/false negatives ($I \setminus (I \cap \hat{I}_{\alpha})$) (and its relative population) and the estimated sources $\hat{I}_{\alpha}$, and on the other hand, the average smallest distance between the set of "false" estimated sources/false positives ($\hat{I}_{\alpha} \setminus (I \cap \hat{I}_{\alpha})$) (and its relative population) and the true sources.

Notice that RLE does not measure any "distance" between true and estimated sources amplitudes, and is of course sensitive to the amplitude thresholding parameter $\alpha$.

B. Spatial Deviation (SD) error

Spatial Deviation (SD) is another distance-based metrics, which interestingly involves an amplitude-weighting. Given a true active source in $r_i$, $SD_i$ measures the spatial spread of the estimated sources around $r_i$, following

$$ SD_i = \sum_{p=1}^{N} |\tilde{s}_p| \|r_p - r_i\|^2, \quad (18) $$

Despite amplitude-weighting, SD is, like RLE, blind to a possible discrepancy between reference and estimated amplitudes.

In the case of a single true source, this metrics reaches its optimal value (zero, without any thresholding) when the estimated sources are concentrated in the true location $r_i$.

In the case of spatially-extended sources, the optimal value of the SD metrics is no longer zero. Indeed, for all true point sources $i_k \in I$ ($\#I \ge 2$), one has, for an optimal reconstruction (with respect to localization) :

$$ SD_{i_k}^{opt} = \sum_{p \in I} |\tilde{s}_p| \|r_p - r_{i_k}\|^2 \quad (19) $$

One may derive lower and upper bounds for this optimal value $SD_{i_k}^{opt}$ from the term $\frac{1}{\#I} \sum_{p \in I} \|r_p - r_{i_k}\|^2 = SD_{i_k}^{ref}$, which corresponds to the SD error for spatially-extended sources with identical amplitude. For a sub-optimal reconstruction with respect to localization, one may show that :

$$ SD_{i_k} = (1 - \beta)SD_{i_k}^{opt} + \beta SD_{i_k} |I^c \quad (20) $$

where $\beta = \frac{(N-\#I)S_{ave}|I^c|}{NS_{ave}}$, and $S_{ave}|I^c$ refers to the average absolute amplitude on the complementary index set $I^c$. Thus, the SD metrics for each point source $i_k \in I$ of a spatially-extended activity is the convex combination of :

* an optimal SD error $SD_{i_k}^{opt}$, where only true sources locations are reconstructed by the model (i.e. no false positives);

* an SD error $SD_{i_k|I^c}$ computed by only taking into account false positives sources $i \notin I$.

The weight $\beta$ in (20) takes into account the relative population of "false positives" ($N-\#I$), and their average amplitude (relatively to the global average). We stress that, unlike the RLE metrics, SD is blind to the position of false negatives, as their zero-amplitude do not affect the SD error value.

To derive a global metric from the local SD value, we shall hereafter consider the maximum difference :

$$ \delta_{SD} = \max_{i} |SD_{i}^{est} - SD_{i}^{ref}| \quad (21) $$

C. Depth metrics and distance between centers of mass

To investigate the possible depth bias phenomenon, well-known for the MNE solution, one may compare the local orthogonal distance to the inner-skull for both estimated and reference sources centers of mass, assuming that the brain curvature has a minimal impact on these measurements.

We stress that these centers of mass may either be restricted or not to the graph describing cortical sources, thus affecting the interpretation of the depth bias. In addition to depth, it is also interesting to keep track of the 3D Euclidean distance between estimated and reference centers of mass, as a large distance may affect the interpretation of the depth metrics. Notice that, even if these metrics take into account the respective estimated and reference amplitudes to compute centers of mass, they do not measure any "distance" between the two amplitudes.

D. Entropy, number of active sources, and spatial extent

Shannon entropy has been used in signal processing literature as a measure of dispersion. Unlike SD, entropy does not measure dispersion away from a reference location, and rather focuses on spreading/concentration. It is well suited to deal with a discrete set of sources amplitudes and possibly not-connected active cortical regions. Indeed, let us consider the discrete set of probabilities $\{\mathbb{P}(X^\Delta \equiv \text{source}_i) = |\tilde{s}_i|\}_{1 \le i \le N}$. Its entropy is given by:

$$ H_{X^\Delta} = - \sum_{i=1}^{N} |\tilde{s}_i| \ln |\tilde{s}_i| \quad (22) $$

and its link with the number of active sources $N$ can be made explicit (using the normalized amplitudes $\tilde{s}_i$) as:

$$ H_{X^\Delta} = \ln(N) + \frac{1}{NS_{ave}} \sum_{i=1}^{N} [S_{ave} \ln(S_{ave}) - |s_i| \ln|s_i|] \quad (23) $$

Thus, Shannon entropy mixes information related to the number $N$ of considered active sources, their average amplitude value $S_{ave}$, and the deviation of each amplitude term $|s_i| \ln|s_i|$ from its average value. Shannon entropy may then be compared to its "reference" counterpart, where $N_{ref}$ active sources are for instance known to share the same amplitude. Indeed, in such a case, one has (with $H_{X^\Delta}^{ref} = \ln(N_{ref}))$:
IV. ANALYSIS OF EVALUATION METRICS

To assess, qualitatively and quantitatively, performances of inversion methods in the distributed sources framework, evaluation metrics are mandatory. A review of some of these can be found in [39]. We will consider here :

* absolute metrics, which bring information on reconstructed sources independently from any ground truth or reference.

* dissimilarity measures, which measure a distance between reconstructed sources and ground truth, or sources reconstructed by different methods.

In both cases, it is important to keep in mind that each metric is sensitive to specific features/information : purely spatial distance-based information (spatial extent of true/estimated sources, spatial dispersion or discrepancy between these two sets...), number of true/estimated active sources, amplitude-based information, or a mix of these features.

Source estimates generally take the form of multivariate time series. In the following, we avoid considering solutions at a fixed time, and focus on time-averaged values, as hereafter detailed : given a space-time matrix $S_* \in \mathbb{R}^{N \times L}$ gathering the amplitude maps for all time snapshots, a time-average $\langle S_*^2 \rangle$ is computed on a prescribed *time segment of interest* denoted by soi, thus leading to an amplitude map $A \in \mathbb{R}^N$ describing the amount of energy contained in the *segment of interest*. We also introduce a normalized version $\tilde{A}$, that sums to one.

$$A = \sqrt{\langle \mathbf{S}_*^2 \rangle_{\text{soi}}}, \quad \tilde{A} = A / \sum_{p=1}^{N} A_p . \quad (14)$$

It is customary to apply a thresholding operation to estimated sources or amplitude maps, to limit the impact of small values. We will refrain as much as possible from doing so, except when using metrics that explicitly involve thresholding. Indeed, one of the key advantages of the w-SBL solver is that it does not require any *a posteriori* thresholding or correction.

In this section, we describe several evaluation metrics that will be used to assess the performance of models. We first introduce some notations. We will denote by $r_i$ the 3D vector position of source vertex $i$ on the discretized cortical surface, and by $|r|$ the 3D Euclidean norm of a vector $r$.

A. Absolute metrics

1) Spatial Dispersion (SD): SD is a distance-based metric, which involves an amplitude-weighting (see [40] for a discussion). Given a cortex location $r_i$, SD$_i$ measures a spatial spread of the estimated sources around $r_i$. We limit ourselves to the center of mass $\langle r \rangle$ of the sources, and define

$$SD = \sum_{p=1}^{N} \tilde{A}_p |r_p - \langle r \rangle|^2 . \quad (15)$$

Remark 3 (Center of mass): One may naturally want to define geometrically the center of mass as $\langle r \rangle_g = \sum_{p=1}^N \tilde{A}_p r_p$ (slight generalization of the definition used in [39]). However, when the amplitude map has large support (as is the case in a large number of methods), this geometrical center of mass can be significantly deeper than the surface, resulting in an overestimation of the spatial dispersion. For that reason, we rather use the center of mass constrained to the discretized

cortical surface (i.e. a center vertex on the mesh grid), as implemented in the MNE-Python software. We denote by $\langle r \rangle$ the center of mass thus defined.

SD vanishes in the case of a point source. For a spatially-extended source, if it happens to exhibit a connected spatial support (e.g. a patch-like shape), SD may be interpreted, with certain care, in terms of the "radius" of an equivalent disk

$$\Delta_{SD} = \sqrt{2SD} \quad (16)$$

2) Entropy: The interpretation of the SD metric as a spatial extent is much less straightforward when it comes to extremely sparse reconstructed sources (e.g. a set of few distant focal sources), or extended sources with disconnected spatial supports. This calls for complementary metrics to measure spreading/concentration, for which we use the entropy.

$$H(\mathcal{A}) = - \sum_{p=1}^{N} \tilde{A}_p \ln \tilde{A}_p . \quad (17)$$

While SD measure dispersion from a reference location, entropy focuses on the actual spreading/concentration without reference to any center of mass or whatsoever. For example, $H(\mathcal{A})$ is maximum when all $\tilde{A}_i$ are equal (maximum spread) and minimum when all $\tilde{A}_i$ vanish except a single one (maximum sparsity). Entropy is complementary to SD: a sparse source with few distant peaks may have a large SD (like a smooth, extended source) but will have a small entropy (unlike the smooth extended source).

3) Center of mass depth: Depth bias, which tends to favor superficial sources against deeper ones, is a common artifact in source localization methods. It is therefore important to be able to estimate the depth of brain sources, and to compare it with ground truth in simulations, and expert knowledge on real data. Like Spatial Dispersion, depth metrics must be treated carefully if reconstructed sources have large support, or exhibit disconnected spatial supports. In this work, we will use the depth of the center of mass of the amplitude map $\mathcal{A}$ (estimated or true) as a metric (see remark 3). We will limit ourselves to centers of mass restricted to the discretized surface.

B. Dissimilarity metrics

A first assessment can be done by comparing the metrics discussed in section IV-A above. For example, compare SD, entropy, center of mass depth or overall amplitude values $\sum_p A_p$, for the tested and reference solutions. We also use standard measures such as binary classification-based measures, Region Localization Error, and investigate a less common complementary measure, based on optimal transport.

Given a reference amplitude map $\mathcal{A}$, denote by $I$ the index set of corresponding active sources ($\mathcal{A}_i \neq 0$), and by $\sharp I$ its cardinality. Given an estimated amplitude map, $\hat{I}_\alpha$ is the index set of active sources after thresholding with threshold $\alpha > 0$ (such that $\tilde{A}_i \ge \alpha$).

1) Region Localization Error: RLE [41] is a commonly used, distance-based dissimilarity metric, based on a threshold $\alpha > 0$, by default a percentage of the maximum value:

$$\text{RLE}(\alpha) = \frac{1}{2} \left[ \frac{1}{\sharp I} \sum_{k \in I} \min_{l \in \hat{I}_\alpha} |r_k - r_l| + \frac{1}{\sharp \hat{I}_\alpha} \sum_{k \in I} \min_{l \in \hat{I}_\alpha} |r_k - r_l| \right] \quad (18)$$
$$ H_{X^\Delta} = H_{X^\Delta}^{\text{ref}} + \ln \left( \frac{N}{N_{\text{ref}}} \right) + \frac{\ln(S_{ave})}{N} \sum_{i=1}^{N} \left( 1 - \frac{|s_i| \ln(|s_i|)}{S_{ave} \ln(S_{ave})} \right) \quad (24) $$

For connected active cortical regions, to link the Shannon entropy to an estimate of the sources spatial extent, one may consider the previous discrete set of probabilities as sorted with respect to the increasing Euclidean distance to the sources center of mass. Then, considering that the spatial discretization steps $\Delta x_i$ between two consecutive sources are uniform, denote by $d_N = N \Delta x$ the distance of the furthest estimated active source with respect to the center of mass. Thus, one may re-write equation (24) as follows:

$$ H_{X^\Delta} = H_{X^\Delta}^{\text{ref}} + \ln \left( \frac{d_N}{d_{\text{ref}}} \right) + \frac{\Delta x \ln(S_{ave})}{d_N} \sum_{i=1}^{N} \left( 1 - \frac{|s_i| \ln(|s_i|)}{S_{ave} \ln(S_{ave})} \right) \quad (25) $$

Shannon entropy thus encapsulates both population-based (or distance-based) and amplitude-based information. Notice, however, that it does not measure any "distance" between $S_{ave}$ and its reference counterpart $S_{ave}^{\text{ref}}$, which does not take part in $H_{X^\Delta}^{\text{ref}}$ computation for identical sources.

### E. Binary classifier metrics: Precision, Recall, F1 score

The models performances may also be assessed using binarization of amplitude maps, through classical scores such as: the precision score, which measures the ability of a classifier not to label as false positive (FP) a sample that is truly negative (TN); the recall, which measures the ability of the classifier to find all the positive samples:

$$ \text{Precision} = \frac{\text{TP}}{\text{TP} + \text{FP}}, \quad \text{Recall} = \frac{\text{TP}}{\text{TP} + \text{FN}} \quad (26) $$

Precision and Recall range between 0 and 1. The F1 score is the harmonic mean of Precision and Recall, and thus deteriorates as soon as one of them decreases:

$$ F1 = 2 \times \frac{\text{Precision} \times \text{Recall}}{\text{Precision} + \text{Recall}} \in [0, 1] \quad (27) $$

These scores are known to be more robust to unbalanced sets than Receiver Operator Curves (and their associated Area Under the Curve).

### F. Wasserstein distance: an Optimal Transport metrics measuring a distance between estimated and reference amplitudes

To take into account a "distance" between estimated and reference sources amplitudes, one may turn towards an Optimal Transport metrics, such as Wasserstein distance [34]. Given:

* a reference distribution $\mu_a = \sum_{i=1}^{N} a_i \delta_{x_i}$, with weights $a_i$ (positive, sum to 1) associated with true sources amplitudes;

* an estimated distribution $\mu_b = \sum_{i=1}^{N} b_i \delta_{y_i}$ with weights $b_i$ (positive, sum to 1) associated with estimated amplitudes; one may compute the minimum cost (with respect to a distance $d$) to move the estimated distribution $\mu_b$ such that it coincides

with the reference one $\mu_a$. This minimum cost is given by the following optimization problem. Given $p \in [1, +\infty[$, set

$$ W_p(d; a, b) = \min_{\gamma \in \mathbb{R}^{N \times N}} \left( \sum_{i,j} \gamma_{i,j} d(x_i, y_j)^p \right)^{\frac{1}{p}} \quad (28) $$

under the constraints that $\gamma.1_{\mathbb{R}^N} = a$ and $\gamma^T.1_{\mathbb{R}^N} = b$; $\gamma \ge 0$. The solution of the above optimization problem (28) results in the Wasserstein distance $W_p(d; a, b)$. In the above definition, $\gamma$ refers to any transportation map from $\mathbb{R}^N$ to $\mathbb{R}^N$ which preserves mass between the distributions $\mu_a$ and $\mu_b$. The "distance" between estimated and reference amplitudes is thus encapsulated within the constraints on $\gamma$, while the geometrical distance between estimated and reference sources is tracked by $d$. No thresholding is required in the definition of the Wasserstein distance, only a normalization of reference and estimated amplitudes. In the following numerical experiments, the 1-Wasserstein distance $W_1$ is used, with respect to the Euclidean distance, to assess the model performances in terms of both amplitude and localization.

## V. NUMERICAL EXPERIMENTS

### A. Wavelet Frame design

For the current paper and the following performance evaluation section, we hereafter consider up to 3 wavelet scales, in addition to the low-frequency component recovered by the scaling function. Following [5], wavelets resulting from the band-pass kernel $g$ are defined in a piecewise manner: cubic spline within the bandwidth, and 2nd degree polynomial increase and decrease outside. This results in a common quality factor ($Q \approx 1.38$), and rather close frame bounds 0.71 and 1.41. As can be seen on Fig. 1, the finest wavelet is designed to that its bandwidth catches the largest eigenvalue $\lambda_{max}$. All frequencies below the green bandwidth are then either encapsulated with the low-frequency component, or taken into account with additional wavelet scales, whose dilation parameters $s$ are spread along a logarithmic scale.

<figure><img src="image_3.png" /><figcaption>Fig. 1. Distribution of spectral graph wavelets along the graph Laplacian spectrum - vertical lines indicate the bandwidth</figcaption></figure>

### B. Experimental setup

In the current paper, we focus on simulated data, and for the sake of reproducibility and future comparisons, we decided to stick as much as possible to datasets and resources included in the MNE-Python software package [6].
RLE only exploits source supports, not amplitudes themselves, and is of course sensitive to the threshold parameter $\alpha$. It may also be shown that RLE only accounts for the set differences $I \setminus \hat{I}_{\alpha}$ and $\hat{I}_{\alpha} \setminus I$ of these supports.

2) **Wasserstein distance:** To better exploit true and estimated source amplitudes (not only their support), and avoid the dependency on a user-defined thresholding parameter, methods based upon optimal transport theory (see [42] for a review) represent, in our opinion, a valuable alternative to RLE.

To make a long story short, an optimal transport distance evaluates the cost of moving quantities from one place to another and seeks the optimal way of doing so. In the context of interest here, we will try to evaluate the cost of moving the estimated amplitude map to the ground truth. This requires

* A prior normalization of amplitudes (from $\mathcal{A}$ to $\mathcal{A}'$, see above), for preventing from global amplitude biases.

* A distance on the domain of interest, here the discretized cortical surface. The Euclidean distance leads to the so-called Wasserstein distances $W_p$.

* The exponent $p$ controls the relative weight of small and large distance values. Here we choose $p = 1$.

Wasserstein distances are evaluated via dedicated numerical algorithms (we refer to [43] for a detailed account).

3) **Binary classifier metrics : Precision, Recall:** Precision and Recall are classical performance scores, defined on binary amplitude maps, which measure either (Precision) the ability of a classifier not to label as false positive (FP) a sample that is truly negative (TN); either (Recall) the ability of the classifier to find all the positive samples. Precision-Recall is a useful measure of success of prediction when the classes are very imbalanced [44]. Precision and Recall range between 0 and 1. The F1 score is the harmonic mean of Precision and Recall, and thus deteriorates as soon as one of them decreases.

In practice, these quantities are computed for comparing estimated maps to ground truth, and therefore require a threshold $\alpha$. They do not fully exploit the amplitude values $\mathcal{A}$.

# V. NUMERICAL EXPERIMENTS

## A. Wavelet Frame design

We hereafter consider wavelet frames with up to 3 wavelet scales, in addition to the low-frequency component recovered by the scaling function. Following [12], wavelets resulting from the band-pass kernel $g$ are defined piecewise : cubic spline within the bandwidth, and 2nd degree polynomial increase and decay outside. This yields a constant quality factor ($Q \approx 1.38$), and rather close frame bounds 0.71 and 1.41 for the considered cortex (MNE-Python open-access sample audivis database). As shown in Fig. 1, the finest wavelet is designed so that its bandwidth catches the largest eigenvalue $\lambda_{max}$ of the cortex graph Laplacian. All frequencies below the green bandwidth are then either encapsulated with the low-frequency component, or taken into account with additional wavelet scales, whose dilation parameters $s$ are spread along a logarithmic scale.

Fig. 2 also displays the normalized amplitudes of the scaling function and wavelets (centered on a vertex within the right primary auditory cortex), for a wavelet frame with 3 levels of

<figure><img src="image_2.png" /><figcaption>Fig. 1. Distribution of spectral graph wavelets along the graph Laplacian spectrum - vertical lines indicate the bandwidth. wavelets are scaled copies of the $g$ kernel (see blue, red and green curves), scaling function is the $h$ kernel (black curve).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Normalized scaling function (top left) and wavelets (with decreasing scales, from top right to bottom right) associated to a frame with 3 levels of resolution (green line: 5% peak-amplitude level curve).</figcaption></figure>

resolution. Scaling functions and wavelets are here displayed with a 5% peak-amplitude thresholding : beyond the green level curve, their amplitude are uniformly below 5% of their peak value, and thus set to zero for visualization purposes. This green contour line allows estimating (through MNE-Python embedded metrics) a label area for each scaling function and wavelets, and thus an equivalent radius. These scaling function and wavelets spatial supports are summarized in Table I.

<table><caption>TABLE I EQUIVALENT RADIUS (MM) FOR SPATIAL SUPPORT OF SPECTRAL GRAPH WAVELETS CENTERED ON A VERTEX WITHIN THE RIGHT PRIMARY AUDITORY CORTEX (5% PEAK-AMPLITUDE THRESHOLDING)</caption><thead><tr><th>Equiv. radius (mm)</th><th>frame N<sub>s</sub> = 1</th><th>frame N<sub>s</sub> = 2</th><th>frame N<sub>s</sub> = 3</th></tr></thead><tbody><tr><td>scaling function</td><td>38.7</td><td>58.1</td><td>60.4</td></tr><tr><td>wavelet #1</td><td>9.9</td><td>39.2</td><td>58.1</td></tr><tr><td>wavelet #2</td><td></td><td>11.2</td><td>39.2</td></tr><tr><td>wavelet #3</td><td></td><td></td><td>11.2</td></tr></tbody></table>

Considering for instance the frame with 3 levels of resolution, one may notice on Fig. 2 that the wavelets spatial support decreases with the wavelets refinement level. This is confirmed by radius data summarized in Table I, and consistent with the fact that, for a given frame, the wavelets frequency support increases with the refinement level.

Interestingly, the spatial support of the finest wavelet is consistent throughout all tested frames, in accordance with the way wavelet frames have been designed : starting with the finest wavelet to catch details, and completing the frame with a scaling function and possibly additional coarser wavelets.
1) *Simulation protocol:* Brain geometry, leadfield and baseline covariance matrices were imported from the "audiovisual" dataset, which contains EEG/MEG evoked data from one subject performing an audiovisual experiment. MNE-Python was also used to simulate a uniform brain activity over patches of two different sizes, generated using parcels from Destrieux atlas (74 × 2 parcels). Examples are given in Fig. 3. Time courses for source activity (see Fig. 2) were given an evoked response shape, still using MNE-python resources. A real spatial covariance matrix was imported from the MNE-python database to simulate baseline data, and added to simulated evoked sensor data. To reproduce the classical 1/f behavior of the noisy signal Power Spectral Density, an AR model was fitted on raw experimental data, and simulated evoked time series were filtered using the corresponding IIR filter. For each case, 320 epochs were simulated, and the trial average was used for the inverse problem.

2) *Compared solutions:* In the results section below, we compared a number of solutions described in Section II together with wavelet based solutions proposed in this paper. The variational reference solutions include MNE (quadratic penalization $f(\mathbf{S}) = \lambda_{\text{MNE}} \|\mathbf{S}\|_F^2$, closed form solution), MCE ($f(\mathbf{S}) = \lambda_{\text{MCE}} \|\mathbf{S}\|_1$, solved using the FISTA algorithm), sVB-SCCD ($f(\mathbf{S}) = \lambda_{\text{MCE}} (\gamma \|\mathbf{S}\|_1 + (1-\gamma) \|\mathbf{S}\|_{TV})$, solved using the ADMM algorithm), with $\gamma = 0.40$ (L1-TV1) and $\gamma = 0.60$ (L1-TV2). These are compared with (sg)w-MCE 3 ($g(\mathbf{A}) = \lambda_{\text{MCE}} \|\mathbf{A}\|_1$, 3 wavelet scales, solved using FISTA) and (sg)w-SBL 1/2/3 (see equations (14-15)), using 1 to 3 wavelet scales and Champagne algorithm (100 iteration steps).

**Remark 2 (Regularization of the noise spatial covariance matrix):** despite a whitening and dimension reduction process on the noise spatial covariance matrix, most inverse problem solvers require a regularization step of the covariance matrix before solving the inverse problem ; interestingly, that is not the case for solvers using Sparse Bayesian Learning (SBL) approaches, in which the cortical source variances are "learned" through an optimization algorithm.

<figure><img src="image_4.png" /><figcaption>Fig. 2. Time course signal designed to mimick an evoked response</figcaption></figure>

For the following model performance evaluation, all EEG and MEG channels were considered as input data for the inverse problem. Given two specific spatial extents (10 and 20 mm from each parcel center, see Fig. 3) for the simulated activity, noisy evoked-like response activities were generated

<figure><img src="image_5.png" /><figcaption>Fig. 3. Reference cortical activity # 1 (over 148) simulated from Destrieux's brain parcellation - extent 10 mm (upper left) and 20 mm (upper right) and SNR values boxplot (lower)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. RLE (upper) and entropy (lower) for the 10 mm case - quartiles $q_1$ (dashed), $q_2$ and $q_3$ (dashed) are represented for each solver</figcaption></figure>

for each of the 148 cortical regions of Destrieux's brain parcellation, with rather similar SNR, as displayed on Fig. 3.

Given this simulation protocol, all the considered solvers were compared and challenged with respect to the previously detailed metrics. To facilitate interpretation and comparison, let us start with all threshold-based metrics.

### Threshold-based metrics (without binarisation)

For all threshold-based metrics (RLE, Shannon entropy, nb. active sources and spatial extent), (sg)w-SBL outperforms significantly the other solvers, from low threshold values, as can be seen on Fig. 4 and 5 for the 10 mm case : (sg)w-SBL
<figure><img src="image_4.png" /><figcaption>Fig. 3. Reconstructed (normalized) sources amplitudes for auditory evoked potentials (left stimulation, raw data processed with MNE-Python pipeline). From top to bottom: Left column : MNE, eLORETA, L1-TV (1 and 2); Right column : SBL and w-SBL (1,2 and 3).</figcaption></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Reconstructed (normalized) sources amplitudes for median nerve evoked potentials (left wrist electric stimulation, raw data processed with BRAINSTORM pipeline). From top to bottom: Left column : MNE, eLORETA, L1-TV (1 and 2) in the left hemisphere (on purpose). Right column : SBL and w-SBL (1, 2 and 3) in the right hemisphere. For each solver, views that display best the estimated activity are chosen.</figcaption></figure>

<figure><img src="image_7.png" /></figure>

## B. Real open-access data

MNE-Python software package [13] and BRAINSTORM application [45] offer several open-access databases, such as auditory evoked potentials (pure tone stimulation) and median nerve stimulation (electric stimulation of left and right wrists). BRAINSTORM software is documented and freely available for download online under the GNU general public license (http://neuroimage.usc.edu/brainstorm). We would like to acknowledge the authors Elizabeth Bock, Esther Florin, Francois Tadel, and Sylvain Baillet for the median nerve dataset. For both protocols, we studied averaged evoked responses (MEG channels only), and compared a number of solutions described in Section II together with wavelet-SBL solutions proposed in this paper. The reference solutions include MNE, eLORETA, MCE (FISTA algorithm), SBL in the sources domain, sVB-SCCD (ADMM algorithm), with $\gamma = 0.3/0.4$ (respectively real/simulated data, tagged L1-TV(1)) and $\gamma = 0.7/0.6$ (respectively real/simulated data, tagged L1-TV(2)). These are compared with w-SBL 1, 2 and 3. (equations (12-13)), using 1 to 3 resolution levels in the wavelet frame and Champagne algorithm (100 iterations).

Fig. 3-4 display the (normalized) amplitudes of reconstructed sources (at peak-amplitude time) obtained with multiple solvers, for respectively a left auditory and a left wrist (i.e. median nerve) stimulation.

Black contour lines indicate respectively the "Aud-rh" label suggested by MNE-Python, and Brodmann parcel BA3b within the right primary somatosensory cortex. The additional red and green lines are level curves corresponding respectively to 1% and 10% source amplitude levels. For clarity, these lines are not represented for some solvers, as they spread across the whole hemisphere. We again stress that no amplitude thresholding is applied on reconstructed sources. We also stress that, for reference MNE and eLORETA solvers em-

bedded within MNE-Python, regularization (hyper)parameters are tuned according to MNE-Python own heuristics. In the auditory case (Fig. 3), most solvers are able to reconstruct higher sources amplitudes in the contralateral hemisphere (where the first sources are detected after stimulation). This is no longer the case for the left wrist stimulation (Fig. 4), for which we have chosen to represent the left hemisphere on the left column, where MNE, eLORETA, L1-TV (1 and 2) reconstruct significant activity, in addition to the expected and here relevant contralateral (right) hemisphere for SBL and w-SBL solvers (complete views are given in Supplementary material). While SBL and w-SBL do satisfy the contralateral property in this median-nerve stimulation, left hemisphere views in Fig. 4 highlight the inability of MNE, eLORETA, and L1-TV (1 and 2) to satisfy the contralateral property. Corresponding MNE/eLORETA/L1-TV results in the right (contralateral) hemisphere are given in Supplementary Material, where the obtained localizations are way more unsatisfactory than SBL and w-SBL. Note that the current results are here obtained with the median nerve evoked data processed following BRAINSTORM [45] detailed pipeline (https://neuroimage.usc.edu/brainstorm/Tutorials/MedianNerveCtf). The MNE-Python processing of the same median nerve raw data leads to slightly different evoked data (https://mne.tools/stable/auto_examples/datasets/brainstorm_data.html), and less unsatisfactory reconstructed sources for classical solvers. We here chose to display results based on the most challenging (BRAINSTORM-processed) evoked data.

In both cases, MNE and eLORETA solutions suffer from the usual spatial leakage, MNE also yielding depth bias towards superficial sources. Sparse L1-TV (1 and 2) solvers manage to mitigate spatial leakage in the auditory case (even if they miss MNE-Python suggested "Aud-rh" label), while, in the
extent 10 mm

<figure><img src="image_7.png" /><figcaption>Fig. 5. Nb. active sources (upper) and spatial extent (lower) - quartiles $q_1$ (dashed), $q_2$ and $q_3$ (dashed) are represented for each solver</figcaption></figure>

extent 20 mm

<figure><img src="image_8.png" /><figcaption>Fig. 7. Nb. active sources (upper) and spatial extent (lower) - quartiles $q_1$ (dashed), $q_2$ and $q_3$ (dashed) are represented for each solver</figcaption></figure>

extent 20 mm

<figure><img src="image_9.png" /><figcaption>Fig. 6. RLE (upper) and entropy (lower) - quartiles $q_1$ (dashed), $q_2$ and $q_3$ (dashed) are represented for each solver</figcaption></figure>

extent 20 mm

<figure><img src="image_10.png" /><figcaption>Fig. 8. Boxplots for geometrical metrics with amplitude-weighting : $\delta_{SD}$ (upper) and distance between 3D centers of mass (lower)</figcaption></figure>

minimizes (with its tight $q_1, q_2, q_3$ quartiles) the RLE metrics, and reaches the fastest the reference value for the other metrics. Notice, however, that (sg)w-SBL 1 deviates from the reference solution when amplitude thresholding exceeds 40%. Similar results are obtained for the 20 mm case (see Fig. 6-7), where (sg)w-SBL reaches the reference extended solution even quicker. As for w-MCE 3, which promotes sparsity in the synthesis wavelet coefficients domain, its outperforms L1-TV solvers with respect to RLE and spatial extent, and MCE with respect to the number of active sources.

## Geometrical metrics with amplitude weighting

Fig. 8 displays, for the 20 mm case, the boxplots obtained for two geometrical metrics with amplitude weighting : the maximum difference between reference and estimated SD vectors ($\delta_{SD}$) and Euclidean distance between 3D centers of mass. Again (sg)w-SBL clearly minimizes these metrics.

Furthermore, given the small distance between centers of mass obtained with (sg)w-SBL, the depth metrics is less affected by the skull curvature. w-MCE 3 solver comes in second position behind (sg)w-SBL, improving a bit MCE results.

Fig. 9 focuses on the depth metrics for the estimated centers of mass (here restricted to the cortex graph). Again, (sg)w-SBL exhibits the best performances, in both cases. Note that one retrieves the classical MNE depth bias towards superficial solutions. L1-TV (1 and 2) solvers also exhibit very unsatisfactory results, with an even stronger bias.

Fig. 10 displays the evolution of the difference between reference and estimated depth (for a c.o.m. restricted to the
left wrist stimulation, this ability is much more sensitive to the tuning of regularization parameters. SBL performs fairly well in terms of localization (peak-amplitude vertex spot on within the target "Aud-rh" label), but gives a highly focal estimate, completely missing the spatial extent.

For both auditory and median nerve datasets, the coupling of Spectral Graph Wavelets with Sparse Bayesian Learning clearly exhibits the best performances in terms of localization, especially with wavelet frames designed with 1 or 2 resolution levels. Additional levels seem to lead to an overestimation of the sources spatial extent. This is consistent with the wavelets widening spatial support summarized in Table I.

To complete this first visual and qualitative performance comparison, Table II summarizes, for left auditory and left wrist stimulations, the results on the absolute metrics (spatial dispersion, depth, entropy) previously introduced. This quantitative comparison also highlights the significantly better performances of w-SBL (especially 1 and 2) on these auditory and median nerve datasets. It also confirms some characteristics of reference solvers, among which the overestimation of spatial extent for MNE, and the underestimation of MCE and SBL. We again stress that these metrics results have been obtained without any thresholding on sources amplitudes. In the absence of ground truth data, the metrics reference values are computed with respect to the primary auditory cortex label "Aud-rh" suggested by MNE-Python, and the Brodmann BA3b parcel associated to the primary somatosensory cortex.

<table><caption>TABLE II
METRICS FOR SOURCE RECONSTRUCTION ON REAL DATA: LEFT AUDITORY (WHITE) AND LEFT WRIST STIMULATIONS (GRAY). FROM TOP TO BOTTOM: $\Delta_{SD}$(MM), CENTER OF MASS DEPTH (MM), ENTROPY. Ref takes the reference parcel as ground truth (1 inside, 0 outside)</caption><thead><tr><th>Ref.</th><th>MNE</th><th>MCE</th><th>SBL</th><th>L1-TV1</th><th>L1-TV2</th><th>wSBL1</th><th>wSBL2</th><th>wSBL3</th></tr></thead><tbody><tr><td>13</td><td>111</td><td>19</td><td>8</td><td>73</td><td>76</td><td><strong>14</strong></td><td><strong>18</strong></td><td>27</td></tr><tr><td>34</td><td>100</td><td>46</td><td>54</td><td>93</td><td>96</td><td><strong>28</strong></td><td><strong>30</strong></td><td><strong>32.2</strong></td></tr><tr><td>15</td><td>0.72</td><td>3.2</td><td>10.5</td><td>1.6</td><td>11.2</td><td>10.5</td><td><strong>17.3</strong></td><td><strong>15.5</strong></td></tr><tr><td>2.8</td><td>16.3</td><td><strong>2.1</strong></td><td><strong>3</strong></td><td>7.8</td><td>16.3</td><td><strong>2.5</strong></td><td>4</td><td>7.8</td></tr><tr><td>4.75</td><td>12.27</td><td>1.86</td><td>1.21</td><td>10.18</td><td>11.79</td><td><strong>5.17</strong></td><td><strong>6.57</strong></td><td><strong>7.19</strong></td></tr><tr><td>6.27</td><td>12.3</td><td>4.27</td><td>2.51</td><td>12.1</td><td>12.4</td><td><strong>5.82</strong></td><td><strong>7.52</strong></td><td><strong>8.35</strong></td></tr></tbody></table>

In Table II, the first two and last rows encapsulate information related to the spatial extent of the reconstructed sources. While w-SBL yields well-localized sources with connected spatial support, MNE, eLORETA and L1-TV (1 and 2) generate, in the median nerve case, disconnected patterns, which significantly impacts the value and interpretation of $\Delta_{SD}$.

Confirming the qualitative analysis, SBL is really good in terms of depth, thanks to its good localization of the sources center of mass. Nevertheless, its seemingly acceptable results, like MCE, with respect to $\Delta_{SD}$ should not be taken for granted: again, results are here affected by the hyper-sparsity and thus excessive amplitudes of the reconstructed sources. For such solvers, entropy appear to be a more robust metric.

MNE and L1-TV are way off with respect to both $\Delta_{SD}$ and entropy. Regarding depth, MNE exhibits a strong bias towards superficial sources (especially in the auditory case) in the qualitative analysis. Quantitative results in terms of depth are here significantly affected by the strong over-estimation of the sources spatial support (with even a disconnected support in the median nerve case), which undermines the proper estimation of the sources center of mass, and consequently its

depth. Regarding L1-TV, the regularization parameters may lead to highly different results (especially in the median nerve case) for sources center of mass and thus depth.

Confirming again the qualitative analysis, w-SBL solvers show the best performances with respect to all absolute metrics considered. Following the spatial support properties of the spectral graph wavelets (see Table I), w-SBL 1 and 2 seem better suited to reconstruct the considered auditory and median nerve sources, given the spatial spread of wavelets #1 and #2 embedded in the associated frames.

### C. Simulated data

In addition to the two previous open-access datasets included in the MNE-Python software package [13], we also evaluated w-SBL and several other solvers on simulated data, generated according to a MNE-Python pipeline.

1) **Simulation protocol:** Brain geometry, leadfield and baseline covariance matrices were imported from the "audiovisual" dataset, which contains EEG/MEG evoked data from one subject performing an audiovisual experiment. MNE-Python was also used to simulate a uniform brain activity over patches of two different sizes (10 and 20 mm from parcel center), generated using parcels from Destrieux atlas (74 × 2 parcels). An example of a 20mm patch can be seen in Fig. 12 (top row). Fig. 5 displays the time course for source activity, designed as an evoked response shape, still using MNE-Python resources. A real spatial covariance matrix was imported from the MNE-Python database to simulate baseline data, added to simulated evoked sensor data. To reproduce the 1/f behavior of the signal Power Spectral Density, an AR model was fitted on raw data, and simulated evoked time series were filtered using the corresponding filter. For each case, 320 epochs were simulated, and the trial average was used for the inverse problem.

<figure><img src="image_8.png" /><figcaption>Fig. 5. Time course signal designed to mimic an evoked response.</figcaption></figure>

2) **Compared solutions and results:** The approaches tested in Section V-B (except eLORETA and regular SBL) were compared on data simulated according to the above protocol.

For performance evaluation, all EEG/MEG channels were considered as input data for the inverse problem. Given two specific spatial extents, noisy evoked-like response activities were generated for each of the 148 cortical regions of Destrieux's brain parcellation, with rather similar SNR : ($q_{25\%} = 2.7, q_{50\%} = 3.6, q_{75\%} = 4.9$) for the 10mm case, and ($q_{25\%} = 2.6, q_{50\%} = 3.7, q_{75\%} = 5.4$) for the 20mm case.

Given this protocol, considered solvers were compared and challenged using the previously detailed metrics (except entropy, less relevant for extended patch-like activities, thus not shown here). Unlike previous real data, the ground truth is now known, which allows for the computation of dissimilarity metrics (in addition to previous absolute metrics) between simulated and reconstructed sources.
extent 20 mm

<figure><img src="image_11.png" /><figcaption>Fig. 9. Boxplots for geometrical metrics with amplitude-weighting : depth of centers of mass (restricted to the graph)</figcaption></figure>

extent 20 mm

<figure><img src="image_12.png" /><figcaption>Fig. 10. Depths of centers of mass (restricted to the graph) with respect to reference depth</figcaption></figure>

extent 20 mm

<figure><img src="image_13.png" /><figcaption>Fig. 11. Statistical scores with respect to threshold</figcaption></figure>

graph) with respect to the reference depth. While MNE and L1-TV solvers clearly exhibit a linearly evolving bias with the reference depth, (sg)w-SBL shows very good and stationary performances throughout all the tested samples and depths. For clarity, MCE and w-MCE 3, which did not exhibit good results in Fig. 9, were not included in Fig. 10.

## Binary classifier metrics

Now considering binary statistical maps built from reference and estimated energy maps, Fig. 11 displays the Precision, Recall and F1 score of the multiple tested solvers. (sg)w-SBL clearly exhibits, by far, and from low threshold values, the best performances with respect to Precision and F1 score. Both scores are rapidly increasing with the threshold value. The recall decay is also more contained with (sg)w-SBL.

## 1-Wasserstein distance between amplitudes maps

To assess solvers performances with respect to amplitude reconstruction, Fig. 12 displays, for the 20 mm case, boxplots for the 1-Wasserstein distance between normalized reference and estimated amplitudes maps. Once again, (sg)w-SBL minimizes by far this amplitude and geometry-based metrics. Note that L1-TV best results with respect to the 1-Wasserstein distance are better than (sg)w-SBL, but once again the huge results dispersion undermines the robustness of this solver.

## $l^2$-norm conservation between amplitudes maps

To further evaluate solvers quantitative performance with respect to amplitude, Fig. 13 displays the ratio of $l^2$ norms between (un-normalized) reference and estimated amplitudes maps averaged on a time segment of interest $A = \sqrt{\langle S_*^2 \rangle_{\text{soi}}}$.

<figure><img src="image_14.png" /><figcaption>Fig. 12. 1-Wasserstein distance boxplot - extent 20 mm</figcaption></figure>

Similarly to the Wasserstein metrics, (sg)w-SBL (especially with 3 wavelet scales) again outperforms the other solvers, with tight $q_1, q_2, q_3$ quartiles very close to the ideal unit value. w-MCE 3 then yields the second best results after (sg)w-SBL with respect to this purely amplitude-based metrics.

## The impact of SNR on L1-TV solvers and hyperparameters

Finally, in order to better comprehend the surprisingly unsatisfactory results obtained with L1-TV solvers, one may investigate the impact of SNR on the solvers performances. Fig. 14 displays SNR impact on the 1-Wasserstein distance between normalized amplitudes maps. Thus, given the considered heuristics for the hyperparameters involved in L1-TV solvers, higher SNR values ($\ge 4$) appear necessary to see results improve from the MNE solution, while (sg)w-SBL performances take the lead from low SNR values.

## Visualization of a reconstructed activity
a) Dissimilarity metrics: w-SBL 1, 2 and 3 minimize (with their tight $q_1$, $q_2$, $q_3$ quartiles) the RLE metrics on Fig. 6, for both spatial extents. The same is true for the 1-Wasserstein distance (see Fig. 7), which doesn't involve any threshold. As already noticed on real data, too many resolution levels in the wavelet frame tend to decrease performances for smaller spatial extents (see Fig. 7, upper figure).

<figure><img src="image_9.png" /><figcaption>RLE - extent 20 mm</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 6. Evaluation of various solvers on simulated data, for two different spatial extent values. RLE metric as a function of the threshold - quartiles $q_1$ (dashed), $q_2$ and $q_3$ (dashed) are represented for each solver.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>$W_1$ distance - extent 20 mm</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 7. Evaluation of various solvers on simulated data, for two different spatial extent values. 1-Wasserstein distance boxplot. The figure is split into two because MCE and w-MCE sometimes yield null solutions for which metrics couldn't be computed, making boxplots less representative.</figcaption></figure>

As for binary classifier metrics, w-SBL again exhibits, by far, and from low threshold values, the best performances with respect to Precision and F1 score (see Fig. 8). Both scores are rapidly increasing with the threshold value. The recall decay is also more contained with w-SBL. w-SBL 3 seems again more adequate for larger spatial extents.

<figure><img src="image_13.png" /><figcaption>Classifier metrics - extent 20 mm</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 8. Evaluation of various solvers on simulated data, for two different spatial extent values. Precision, Recall and F1 with respect to threshold.</figcaption></figure>

It is also interesting to quantify the impact of Signal to Noise Ratio (SNR) on performances. In this respect, 1-Wasserstein distance (Section IV-B) is an ideal scope. Fig. 9 shows that w-SBL solvers quickly take the lead among MNE and L1-TV solutions, from low SNR values, which we interpret as follows:
<figure><img src="image_15.png" /><figcaption>Fig. 13. $l^2$-norm ratio boxplot on amplitudes $A = \sqrt{\langle S_*^2 \rangle_{SOI}}$ - extent 20 mm</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14. 1-Wasserstein distance with respect to SNR - extent 20 mm</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 15. Reference (upper) and reconstructed (lower: w-SBL 3 then L1-TV2) activity - patch 128 - parietal and medial views - extent 20 mm</figcaption></figure>

Finally, Fig. 15 displays the visualization of a reference and reconstructed activity (time-dependent absolute amplitude), for patch #128 (SNR = 4.14, 20 mm extent case), which corresponds to the first quartile for the w-SBL 3 solver with respect to the 1-Wasserstein distance. For snapshots taken at peak-amplitude time, (sg)w-SBL recovers the true source location, with a good spatial extent and a time course close to the reference one. L1-TV solution deviates notably from the true location (from medial to parietal view), and is not able to reproduce the reference time course.

## C. Discussion

From these extensive simulations of extended (non-focal) sources, one may note that, as expected, (w)MNE tends to produce source estimates that are too diffuse, and require post-processing (not addressed here) to yield fairly reliable source localization. MCE produces source estimates that are extremely focal, and is therefore not suitable for extended sources. Surprizingly enough, L1-TV approaches such as sVB-SCCD turn out to be extremely unstable, and often yield source estimates that are too extended and superficial. The simulation protocol was expected to favor these approaches. Our interpretation is that the two hyper-parameters involved there are difficult to tune, and the heuristics used in this work is not sufficient (to our knowledge, no systematic approach has been proposed in the literature). Wavelet-based solvers are indeed able to reconstruct such sources, more accurately and robustly than the other tested approaches. More specifically, (sg)w-MCE ($l^1$ penalization on wavelet coefficients) yields fairly good estimates, that compare favourably to the state of

the art for most evaluation metrics tested here. It however requires choosing an hyperparameter value, for which the heuristics derived in this paper seems relevant. The sparse Bayesian learning (sgw-SBL) approach, which provides an estimate for the wavelet covariance matrix together with wavelet coefficients, turns out to outperform consistently and significantly other approaches, with respect to a set of metrics designed to capture various aspects and possible artifacts. Quite remarkably, the sgw-SBL approach tends to perform quite well in terms of depth bias, which is a major difficulty in the M/EEG inverse problem.

We also plan to study similar approaches that use (orthonormal) wavelet bases on the cortical surface, instead of overcomplete wavelet frames as we do here. In the sgw case, scaling coefficients (low spatial frequencies) turn out to be almost systematically set to zero or very small values by the solvers (sgw-MCE or sgw-SBL). While relevant in practice, we conjecture that this could be a side effect of the redundancy of the wavelet frame, which we would like to confirm. We also plan to study the problem of hyperparameters tuning, in particular for L1-TV models.

## VI. CONCLUSIONS

This paper proposed a new approach for the M/EEG inverse problem, in the framework of the distributed sources model, which exploits sparsity in the spectral graph wavelet (spatial) domain. Such wavelets forming an overcomplete system, wavelet expansion of the solution is not unique. We have chosen to focus on the so-called synthesis approach, where
while SBL framework allows for a systematic tuning and optimization of *N* different regularization (hyper)parameters (for the *N* different sources) given the current SNR level (estimated on the whole time window), sparse non-linear solvers such as L1-TV do not even have a well-accepted heuristic for tuning these (hyper)parameters. As a result, higher SNR values (≥ 4) appear necessary to see L1-TV results improve from the MNE solution. These higher SNR values classically allow prioritizing the data fitting term $||Z - GS||^2$ over the regularization term $f(S)$. Regarding w-SBL solvers, w-SBL 3 again exhibits better performances for the larger spatial extent. In the 10 mm extent case, w-SBL 1 and 2 are slightly better, again in accordance with Table I results.

<figure><img src="image_15.png" /><figcaption>Fig. 9. Evaluation of various solvers on simulated data, for two different spatial extent values. 1-Wasserstein distance with respect to SNR</figcaption></figure>

b) Absolute metrics: Starting with the center of mass depth bias, w-SBL solvers once again exhibit better performances (see Fig. 10), by properly and consistently estimating the center of mass depth, even for deep reference sources (and for both spatial extents). MNE and L1-TV solutions, on the other hand, exhibit the well-known linearly evolving depth bias phenomenon towards superficial sources. For clarity, MCE and w-MCE 3 solvers, which did not produce exploitable results in almost half the tested patches, were not included in Fig. 10.

To gain insight on the sources spatial extent, we now turn to the spatial dispersion-based "radius" $\Delta_{SD}$, and entropy. In both extent cases, w-SBL exhibits, by far, the best performances with respect to $\Delta_{SD}$ "radius" metric, with narrow boxplots, as can be seen on Fig. 11. We recall that spatial dispersion and $\Delta_{SD}$ do not require any amplitude-thresholding on the sources. As already noticed, additional levels of resolution in the wavelet frame tend to decrease performances for the smaller extent scenario, and appear more suited for wider sources. As previously outlined in the evaluation metrics Section IV, $\Delta_{SD}$ can only be interpreted as a radius with certain care, and it is thus not surprising to witness a bias between the pre-defined reference sources extents and $\Delta_{SD}$.

<figure><img src="image_16.png" /><figcaption>Fig. 10. Evaluation of various solvers on simulated data, for two different spatial extent values. Center of mass depths bias (restricted to the graph) as a function of the reference depth</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 11. Evaluation of various solvers on simulated data, for two different spatial extent values. $\Delta_{SD}$ radius metric boxplots</figcaption></figure>

c) Visualization of a reconstructed activity: Fig. 12 displays the visualization of a reference and reconstructed activity (time-dependent absolute amplitude), for patch #128 (SNR = 4.14, 20 mm extent case), corresponding to the first quartile for the w-SBL 3 solver with respect to the 1-Wasserstein distance. For snapshots taken at peak-amplitude time, w-SBL recovers the true source location, with a good spatial extent. L1-TV solution deviates notably from the true location (from medial to parietal view), and here exhibits a significant depth bias.
optimization is performed with respect to wavelet coefficients. In variational approaches, the unknown (solutions, or wavelet coefficients) are obtained as solutions of an optimization problem, for which hyperparameters have to be chosen. The sparse Bayesian learning paradigm was adapted to the considered wavelet setting, leading to an optimization problem where hyperparameters are also optimized. These approaches were compared with standard ones using numerical simulations of joint EEG/MEG data provided by the MNE-python software. These extensive simulations and the numerous and discussed evaluation metrics show that the combination of spatial wavelets with sparse Bayesian learning brings a real improvement for extended brain sources reconstruction, both in terms of localization and quantitative estimation.

In a forthcoming paper, we will exploit the proposed solutions on real data for specific applications, and investigate evaluation of M/EEG source localization using companion intracerebral data.

Being able to identify localization, depth and time course of brain activity is important in clinical applications for all diseases that require source localization from MEG or EEG measurements. Our results show that the proposed approach improves significantly the state of the art for the reconstruction of extended brain sources.

## REFERENCES

[1] S. Baillet, J. C. Mosher, and R. M. Leahy, "Electromagnetic brain mapping," *IEEE Signal Processing Magazine*, vol. 18, no. 6, pp. 14–30, 2001.

[2] S. Baillet, *Encyclopedia of Computational Neuroscience*. New York, NY: Springer New York, 2014, ch. Forward and Inverse Problems of MEG/EEG, pp. 1–8. [Online]. Available: https://doi.org/10.1007/978-1-4614-7320-6_529-1

[3] A. Sorrentino and M. Piana, "Inverse modeling for MEG/EEG data," *bioRxiv*, 2017. [Online]. Available: https://www.biorxiv.org/content/early/2017/05/08/135269

[4] W.-T. Chang et al., "Spatially sparse source cluster modeling by compressive neuromagnetic tomography," *NeuroImage*, vol. 53, no. 1, pp. 146–160, 2010. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1053811910007196

[5] D. K. Hammond, P. Vanderghynst, and R. Gribonval, "Wavelets on graphs via spectral graph theory," *Applied and Computational Harmonic Analysis*, vol. 30, no. 2, pp. 129–150, 2011. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S1063520310000552

[6] A. Gramfort et al., "MEG and EEG data analysis with MNE-Python," *Frontiers in Neuroscience*, vol. 7, p. 267, 2013. [Online]. Available: https://www.frontiersin.org/article/10.3389/fnins.2013.00267

[7] S. Mokhtari et al., "Sparse wavelet-based solutions for the M/EEG inverse problem," in *SampTA, International Conference on Sampling Theory and Applications*, 2023. [Online]. Available: https://arxiv.org/abs/2306.15262

[8] J. Mosher, P. Lewis, and R. Leahy, "Multiple dipole modeling and localization from spatio-temporal MEG data," *IEEE Transactions on Biomedical Engineering*, vol. 39, no. 6, pp. 541–557, June 1992.

[9] K. Jerbi et al., "Localization of realistic cortical activity in MEG using current multipoles," *NeuroImage*, vol. 22, no. 2, pp. 779–793, 2004. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S1053811904001028

[10] B. Van Veen et al., "Localization of brain electrical activity via linearly constrained minimum variance spatial filtering," *IEEE Transactions on Biomedical Engineering*, vol. 44, no. 9, pp. 867–880, 1997.

[11] K. Sekihara and S. S. Nagarajan, *Adaptive Spatial Filters for Electromagnetic Brain Imaging*, ser. Series in Biomedical Engineering, J. H. Nagel, Ed. Springer-Verlag Berlin Heidelberg, 2008.

[12] M. S. Hämäläinen and R. J. Ilmoniemi, "Interpreting magnetic fields of the brain: minimum norm estimates," *Medical and Biological Engineering & Computing*, vol. 32, no. 1, pp. 35–42, 1994.

[13] M. Fuchs et al., "Linear and nonlinear current density reconstructions," *Journal of clinical Neurophysiology*, vol. 16, no. 3, pp. 267–295, 1999.

[14] F.-H. Lin et al., "Assessing and improving the spatial accuracy in MEG source localization by depth-weighted minimum-norm estimates," *NeuroImage*, vol. 31, no. 1, pp. 160–171, 2006. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S1053811905024973

[15] R. D. Pascual-Marqui, "Discrete, 3d distributed, linear imaging methods of electric neuronal activity, part 1: exact, zero error localization," 2007.

[16] A. M. Dale et al., "Dynamic statistical parametric mapping: Combining fMRI and meg for high-resolution imaging of cortical activity," *Neuron*, vol. 26, no. 1, pp. 55–67, 2000. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0896627300811381

[17] R. Pascual-Marqui, "Standardized low resolution brain electromagnetic tomography (SLORETA): Technical details," *Methods and findings in experimental and clinical pharmacology*, vol. 24 Suppl D, pp. 5–12, 2002.

[18] H. Becker et al., "A performance study of various brain source imaging approaches," in *IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*, 2014, pp. 5869–5873. [Online].

[19] K. Matsuura and Y. Okabe, "Selective minimum-norm solution of the biomagnetic inverse problem," *IEEE Transactions on Biomedical Engineering*, vol. 42, no. 6, pp. 608–615, June 1995.

[20] K. Uutela, M. Hämäläinen, and E. Somersalo, "Visualization of magnetoencephalographic data using minimum current estimates," *NeuroImage*, vol. 10, no. 2, pp. 173–180, 1999.

[21] L. Ding, "Reconstructing cortical current density by exploring sparseness in the transform domain," *Physics in Medicine and Biology*, vol. 54, no. 9, pp. 2683–2697, apr 2009. [Online]. Available: https://doi.org/10.1088/0031-9155/54/9/006

[22] H. Becker et al., "Fast, variation-based methods for the analysis of extended brain sources," in *22nd European Signal Processing Conference (EUSIPCO)*, 2014, pp. 41–45.

[23] Y. Li et al., "s-smooth: Sparsity and smoothness enhanced EEG brain tomography," *Frontiers in Neuroscience*, vol. 10, p. 543, 2016. [Online]. Available: https://www.frontiersin.org/article/10.3389/fnins.2016.00543

[24] W. Ou, M. S. Hämäläinen, and P. Golland, "A distributed spatio-temporal EEG/MEG inverse solver," *NeuroImage*, vol. 44, no. 3, pp. 932–946, 2009. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S1053811908007155

[25] A. Gramfort, M. Kowalski, and M. Hämäläinen, "Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods," *Physics in Medicine and Biology*, vol. 57, no. 7, pp. 1937–1961, mar 2012. [Online]. Available: https://doi.org/10.1088/0031-9155/57/7/1937

[26] A. Gramfort et al., "Time-frequency mixed-norm estimates: Sparse M/EEG imaging with non-stationary source activations," *NeuroImage*, vol. 70, pp. 410–422, 2013. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S1053811912012372

[27] D. Wipf and S. Nagarajan, "A unified bayesian framework for MEG/EEG source imaging," *NeuroImage*, vol. 44, no. 3, pp. 947–966, Feb. 2009. [Online]. Available: https://pubmed.ncbi.nlm.nih.gov/18602278

[28] D. P. Wipf, B. D. Rao, and S. Nagarajan, "Latent variable bayesian models for promoting sparsity," *IEEE Trans. Inf. Theor.*, vol. 57, no. 9, p. 6236–6255, sep 2011. [Online]. Available: https://doi.org/10.1109/TIT.2011.2162174

[29] A. Hashemi et al., "Unification of sparse bayesian learning algorithms for electromagnetic brain imaging with the majorization minimization framework," *NeuroImage*, vol. 239, p. 118309, 2021. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1053811921005851

[30] D. R. Hunter and K. Lange, "A tutorial on MM algorithms," *The American Statistician*, vol. 58, no. 1, pp. 30–37, 2004. [Online]. Available: https://doi.org/10.1198/00031300428236

[31] C. Amblard, E. Lapalme, and J.-M. Lina, "Biomagnetic source detection by maximum entropy and graphical models," *IEEE Transactions on Biomedical Engineering*, vol. 51, no. 3, pp. 427–442, March 2004.

[32] J.-M. Lina et al., "Wavelet-based localization of oscillatory sources from magnetoencephalography data," *IEEE Transactions on Biomedical Engineering*, vol. 61, no. 8, pp. 2350–2364, 2014.

[33] M.-C. Roubaud et al., "Space-Time Extension of the MEM Approach for Electromagnetic Neuroimaging," in *IEEE International Workshop on Machine Learning for Signal Processing*, Aalborg, Denmark, Sep. 2018. [Online]. Available: https://hal.archives-ouvertes.fr/hal-01843557

[34] S. Kolouri et al., "Optimal mass transport: Signal processing and machine-learning applications," *IEEE Signal Processing Magazine*, vol. 34, no. 4, pp. 43–59, 2017.
<figure><img src="image_18.png" /><figcaption>Fig. 12. Reconstruction of simulated activity for a specific patch, with various solvers. From top to bottom: reference, w-SBL 3, L1-TV(2) solutions - patch 128 - parietal and medial views - extent 20 mm.</figcaption></figure>

## D. Discussion

From these extensive simulations of extended (non-focal) sources, one may note that, as expected, (w)MNE tends to produce source estimates that are too diffuse, and require post-processing (not addressed here) to yield fairly reliable source localization. MCE produces source estimates that are extremely focal, and is therefore not suitable for extended sources. Surprizingly enough, L1-TV approaches such as sVB-SCCD turn out to be extremely unstable and sensitive to regularization parameters, and often yield source estimates that are too extended and superficial. The simulation protocol (with patch-like uniform sources activities) was expected to favor these approaches. Our interpretation is that the two regularization (hyper)parameters involved there are difficult to tune, and the heuristics used in this work is not sufficient (to our knowledge, no systematic approach has been proposed in literature). Wavelet-based solvers are indeed able to reconstruct such sources, more accurately and robustly than the other tested approaches. More specifically, w-MCE (ℓ¹ penalization on wavelet coefficients) yields fairly good estimates, that compare favourably to the state of the art for most evaluation metrics tested here, but for only half the tested patches. w-MCE moreover requires the tuning of a regularization (hyper)parameter, for which the heuristics derived in this paper seems relevant. The sparse Bayesian learning (w-SBL) approach, which provides an estimate for the wavelet covariance matrix together with wavelet coefficients, turns out to outperform consistently and significantly (without any amplitude thresholding) other approaches, from low SNR values and with respect to a set of metrics designed to capture various aspects and possible artifacts. Quite remarkably, w-SBL performs quite well in terms of depth bias, which is a major difficulty in the M/EEG inverse problem.

With spectral graph wavelets, scaling coefficients (low spatial frequencies) turn out to be almost systematically set to zero or very small values by the solvers (w-MCE or w-SBL). While relevant in practice, we conjecture that this could be a side effect of the redundancy of the wavelet frame. In this respect, coupling SBL with wavelet bases instead of overcomplete frames will be a relevant comparison.

Besides comparison with wavelet bases, the impact of the graph wavelet parameters on the solution also deserves further investigations. This includes in particular the number of decomposition levels, but also parameters that were fixed in this study, such as the filter bandwidths, and the choice of the Laplacian matrix. Regarding the number of decomposition levels, simulations show that increasing the number $N_s$ of decomposition levels tends to increase the spatial extent of the solution (see e.g. Fig. 11). This would for example impact a choice of threshold if one wants to estimate spatial extent by thresholding the solution (Fig. 8 shows that the optimal F1-based threshold would depend on $N_s$). One may then consider that $N_s$ should be interpreted as a new hyper-parameter, since the spatial extent is not known in advance. Even though w-SBL consistently outperforms other approaches, whatever the value of $N_s$ for all considered metrics, this point deserves further investigations.

## VI. CONCLUSIONS

This paper proposed a new approach for the M/EEG inverse problem, in the framework of the distributed sources model, which exploits sparsity in the spectral graph wavelet (spatial) domain. Such wavelets forming an overcomplete system, wavelet expansion of the solution is not unique. We have chosen to focus on the so-called synthesis approach, where optimization is performed with respect to wavelet coefficients. In variational approaches, the unknown (solutions, or wavelet coefficients) are obtained as solutions of an optimization problem, for which regularization (hyper)parameters have to be chosen. The sparse Bayesian learning paradigm was adapted to the considered wavelet setting, leading to an optimization problem where regularization (hyper)parameters are also optimized. These approaches were compared with standard ones using real open-access datasets and numerical simulations of joint EEG/MEG data provided by the MNE-Python software. Results on a third real dataset are available as Supplementary Material. These extensive simulations and the numerous discussed evaluation metrics (with and without amplitude thresholding) show that the combination of spatial wavelets with sparse Bayesian learning brings a real improvement for extended brain sources reconstruction, both in terms of localization and quantitative estimation. A natural extension of the simulations would be to consider multiple sources, with partially correlated time courses, for which additional difficulties may be expected.

In a forthcoming paper, we will investigate evaluation of M/EEG source localization using companion intracerebral data. This requires the inclusion of additional wavelets, built from deeper 3D structures instead of surface, spectral graph wavelets being able to handle both within the same framework.
We are also currently working at a quantitative assessment of the impact of the wavelet design parameters (number of scales, Laplacian matrix, filter bandwidth). At this point, we do not believe we can estimate the optimal number of decomposition levels directly from data, but we can try to make the result less sensitive to that parameter. One may also imagine that in some applications, an expert could be interested to compare the estimated extent to prior knowledge, and therefore be able to make a choice, so that this apparent drawback may become valuable.

## REFERENCES

[1] S. Baillet, J. C. Mosher, and R. M. Leahy, "Electromagnetic brain mapping," *IEEE Signal Process. Mag.*, vol. 18, no. 6, pp. 14–30, 2001.

[2] S. Baillet, *Encyclopedia of Computational Neuroscience*. Springer, 2014, ch. Forward and Inverse Problems of MEG/EEG, pp. 1–8.

[3] A. Sorrentino and M. Piana, *Inverse Modeling for MEG/EEG Data*. Cham: Springer International Publishing, 2017, pp. 239–253.

[4] B. He, A. Sohrabpour, E. Brown, and Z. Liu, "Electrophysiological source imaging: a noninvasive window to brain dynamics," *Annu. Rev. Biomed. Eng.*, vol. 20, p. 171–196, 2018.

[5] K. Liao, M. Zhu, L. Ding, S. Valette, W. Zhang, and D. Dickens, "Sparse imaging of cortical electrical current densities via wavelet transforms," *Phys. Med. Biol.*, vol. 57, no. 21, pp. 6881–6901, oct 2012.

[6] L. Ding, "Reconstructing cortical current density by exploring sparseness in the transform domain," *Phys. Med. Biol.*, vol. 54, no. 9, pp. 2683–2697, apr 2009.

[7] M. Vetterli, J. Kovacevic, and V. K. Goyal, *Fourier and Wavelet Signal Processing*. Cambridge University Press, 2014.

[8] P. Berkes, R. Turner, and M. Sahani, "On sparsity and overcompleteness in image models," in *NeurIPS*, J. Platt, D. Koller, Y. Singer, and S. Roweis, Eds., vol. 20. Curran Associates, Inc., 2007.

[9] W.-T. Chang, A. Nummenmaa, J.-C. Hsieh, and F.-H. Lin, "Spatially sparse source cluster modeling by compressive neuromagnetic tomography," *NeuroImage*, vol. 53, no. 1, pp. 146–160, 2010.

[10] K. Liao, M. Zhu, and L. Ding, "A new wavelet transform to sparsely represent cortical current densities for EEG/MEG inverse problems," *Comput. Methods Programs Biomed.*, vol. 111, no. 2, pp. 376–388, 2013.

[11] M. Zhu, W. Zhang, D. Dickens, and L. Ding, "Reconstructing spatially extended brain sources via enforcing multiple transform sparseness," *Neuroimage*, vol. 86, p. 280–293, 2014.

[12] D. K. Hammond, P. Vandergheynst, and R. Gribonval, "Wavelets on graphs via spectral graph theory," *Appl. Comput. Harmon. Anal.*, vol. 30, no. 2, pp. 129 – 150, 2011.

[13] A. Gramfort, M. Luessi, E. Larson, D. Engemann, D. Strohmeier, C. Brodbeck, R. Goj, M. Jas, T. Brooks, L. Parkkonen, and M. Hämäläinen, "MEG and EEG data analysis with MNE-Python," *Frontiers in Neuroscience*, vol. 7, p. 267, 2013.

[14] S. Mokhtari, J.-M. Badier, C. G. Bénar, and B. Torrésani, "Sparse wavelet-based solutions for the M/EEG inverse problem," in *SampTA, International Conference on Sampling Theory and Applications*, 2023.

[15] J. Mosher, P. Lewis, and R. Leahy, "Multiple dipole modeling and localization from spatio-temporal MEG data," *IEEE Trans. Biomed. Eng.*, vol. 39, no. 6, pp. 541–557, June 1992.

[16] K. Jerbi, S. Baillet, J. Mosher, G. Nolte, L. Garnero, and R. Leahy, "Localization of realistic cortical activity in MEG using current multipoles," *NeuroImage*, vol. 22, no. 2, pp. 779 – 793, 2004.

[17] B. Van Veen, W. Van Drongelen, M. Yuchtman, and A. Suzuki, "Localization of brain electrical activity via linearly constrained minimum variance spatial filtering," *IEEE Trans. Biomed. Eng.*, vol. 44, no. 9, pp. 867–880, 1997.

[18] K. Sekihara and S. S. Nagarajan, *Adaptive Spatial Filters for Electromagnetic Brain Imaging*, ser. Series in Biomedical Engineering, J. H. Nagel, Ed. Springer-Verlag Berlin Heidelberg, 2008.

[19] H. Sah, A. Sohrabpour, M. Akçakaya, and B. He, "Electromagnetic brain source imaging by means of a robust minimum variance beamformer," *IEEE Trans. Biomed. Eng.*, vol. 65, p. 2365–2374, 2018.

[20] M. S. Hämäläinen and R. J. Ilmoniemi, "Interpreting magnetic fields of the brain: minimum norm estimates," *Med. Biol. Eng. Comput.*, vol. 32, no. 1, pp. 35–42, 1994.

[21] M. Fuchs, M. Wagner, T. Köhler, and H.-A. Wischmann, "Linear and nonlinear current density reconstructions," *J. Clin. Neurophysiol.*, vol. 16, no. 3, pp. 267–295, 1999.

[22] F.-H. Lin, T. Witzel, S. P. Ahlfors, S. M. Stufflebeam, J. W. Belliveau, and M. S. Hämäläinen, "Assessing and improving the spatial accuracy in MEG source localization by depth-weighted minimum-norm estimates," *NeuroImage*, vol. 31, no. 1, pp. 160 – 171, 2006.

[23] R. D. Pascual-Marqui, "Discrete, 3d distributed, linear imaging methods of electric neuronal activity. part 1: exact, zero error localization," 2007.

[24] K. Matsuura and Y. Okabe, "Selective minimum-norm solution of the biomagnetic inverse problem," *IEEE Trans. Biomed. Eng.*, vol. 42, no. 6, pp. 608–615, June 1995.

[25] K. Uutela, M. Hämäläinen, and E. Somersalo, "Visualization of magnetoencephalographic data using minimum current estimates," *NeuroImage*, vol. 10, no. 2, pp. 173–180, 1999.

[26] H. Becker, L. Albera, P. Comon, R. Gribonval, F. Wendling, and I. Merlet, "A performance study of various brain source imaging approaches," in *IEEE Int Conf Acoust Speech Signal Process (ICASSP)*, 2014, pp. 5869–5873.

[27] A. Sohrabpour and B. He, "Exploring the extent of source imaging: Recent advances in noninvasive electromagnetic brain imaging," *Curr. Opin. Biomed. Eng.*, vol. 18, p. 100277, 2021.

[28] N. von Ellenrieder, L. Beltrachini, C. Muravchik, and J. Gotman, "Extent of cortical generators visible on the scalp: effect of a subdural grid," *NeuroImage*, vol. 101, p. 787–795, 2014.

[29] Z.-L. Lü and S. Williamson, "Spatial extent of coherent sensory evoked cortical activity," *Exp Brain Res*, vol. 84, p. 411–416, 1991.

[30] H. Becker, L. Albera, P. Comon, R. Gribonval, and I. Merlet, "Fast, variation-based methods for the analysis of extended brain sources," in *22nd Eur. Signal Process. Conf (EUSIPCO)*, 2014, pp. 41–45.

[31] L. Ding and B. He, "Sparse source imaging in electroencephalography with accurate field modeling," *Hum. Brain Mapp.*, no. 29, pp. 1053–1067, 2008.

[32] D. Wipf and B. Rao, "Sparse bayesian learning for basis selection," *IEEE Trans. Signal Process.*, vol. 52, no. 8, pp. 2153–2164, Aug 2004.

[33] D. P. Wipf, J. P. Owen, H. T. Attias, K. Sekihara, and S. S. Nagarajan, "Robust bayesian estimation of the location, orientation, and time course of multiple correlated neural sources using MEG," *NeuroImage*, vol. 49, no. 1, pp. 641–655, 2010.

[34] A. Hashemi, C. Cai, G. Kutyniok, K.-R. Müller, S. S. Nagarajan, and S. Haufe, "Unification of sparse bayesian learning algorithms for electromagnetic brain imaging with the majorization minimization framework," *NeuroImage*, vol. 239, p. 118309, 2021.

[35] D. Wipf and S. Nagarajan, "A unified bayesian framework for MEG/EEG source imaging," *NeuroImage*, vol. 44, no. 3, pp. 947–966, Feb. 2009.

[36] D. R. Hunter and K. Lange, "A tutorial on MM algorithms," *The American Statistician*, vol. 58, no. 1, pp. 30–37, 2004.

[37] V. P. Oikonomou and I. Kompatsiaris, "A novel bayesian approach for EEG source localization," *Comput. Intell. Neurosci*, vol. 2020, p. 8837954, Oct. 2020.

[38] S. Atasoy, I. Donnelly, and J. Pearson, "Human brain networks function in connectome-specific harmonic waves." *Nat. Commun.*, vol. 7, p. 10340, Jan 2016.

[39] J. G. Samuelsson, N. Peled, F. Mamashli, J. Ahveninen, and M. S. Hämäläinen, "Spatial fidelity of MEG/EEG source estimates: A general evaluation approach," *NeuroImage*, vol. 224, p. 117430, 2021.

[40] M. Stenroos and O. Hauk, "Minimum-norm cortical source estimation in layered head models is robust against skull conductivity error," *NeuroImage*, vol. 81, pp. 265–272, 2013.

[41] A. Molins, S. Stufflebeam, E. Brown, and M. Hämäläinen, "Quantification of the benefit from integrating MEG and EEG data in minimum $\ell^2$-norm estimation," *NeuroImage*, vol. 42, no. 3, pp. 1069 – 1077, 2008.

[42] S. Kolouri, S. R. Park, M. Thorpe, D. Slepcev, and G. K. Rohde, "Optimal mass transport: Signal processing and machine-learning applications," *IEEE Signal Process. Mag.*, vol. 34, no. 4, pp. 43–59, 2017.

[43] G. Peyré and M. Cuturi, "Computational Optimal Transport," *Found. Trends Mach. Learn.*, vol. 11, no. 5-6, pp. 355–206, 2018.

[44] T. Saito and M. Rehmsmeier, "The precision-recall plot is more informative than the ROC plot when evaluating binary classifiers on imbalanced datasets," *PloS one*, vol. 10, no. 3, p. e0118432, 2015.

[45] F. Tadel, S. Baillet, J. C. Mosher, D. Pantazis, and R. M. Leahy, "Brainstorm: A user-friendly application for MEG/EEG analysis," *Comput. Intell. Neurosci*, vol. 2011, p. 879716, 2011.