Alexandre Gramfort
email: alexandre.gramfort@inria.fr

Matthieu Kowalski

Matti Hämäläinen

Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods

Keywords: Magnetoencephalography, Electroencephalography, inverse problem, convex optimization, mixed-norms, structured sparsity, functional brain imaging

Magneto-and electroencephalography (M/EEG) measure the electromagnetic fields produced by the neural electrical currents. Given a conductor model for the head, and the distribution of source currents in the brain, Maxwell's equations allow one to compute the ensuing M/EEG signals. Given the actual M/EEG measurements and the solution of this forward problem, one can localize, in space and in time, the brain regions than have produced the recorded data. However, due to the physics of the problem, the limited number of sensors compared to the number of possible source locations, and measurement noise, this inverse problem is ill-posed. Consequently, additional constraints are needed. Classical inverse solvers, often called Minimum Norm Estimates (MNE), promote source estimates with a small ℓ 2 norm. Here, we consider a more general class of priors based on mixed-norms. Such norms have the ability to structure the prior in order to incorporate some additional assumptions about the sources. We refer to such solvers as Mixed-Norm Estimates (MxNE). In the context of M/EEG, MxNE can promote spatially focal sources with smooth temporal estimates with a twolevel ℓ 1 /ℓ 2 mixed-norm, while a three-level mixed-norm can be used to promote spatially non-overlapping sources between different experimental conditions. In order to efficiently solve the optimization problems of MxNE, we introduce fast first-order iterative schemes that for the ℓ 1 /ℓ 2 norm give solutions in a few seconds making such a prior as convenient as the simple MNE. Furhermore, thanks to the convexity of the optimization problem, we can provide optimality conditions that guarantee global convergence. The utility of the methods is demonstrated both with simulations and experimental MEG data.

Introduction

Inverse problems are common in applied physics. They consist of estimating the parameters of a model from incomplete and noisy measurements. Examples are tomography which is a key technology in the field of medical imaging, or identifying the targets using sonars and radars. Blind source separation, which is an active topic of research in audio-processing, also falls in this category. In this contribution, we target the localization of the brain regions whose neural activations produce electromagnetic fields measured by Magnetoencephalography (MEG) and Electroencephalography (EEG), which we will refer to collectively as M/EEG. The sources of M/EEG are current generators classically modeled by current dipoles. Given a limited number of noisy measurements of the electromagnetic fields associated to neural activity, the task is to estimate the positions and amplitudes of the sources that have generated the signals. By solving this problem, M/EEG become noninvasive methods for functional brain imaging with a high temporal resolution.

Finding a solution to an inverse problem requires finding a good model for the observed data given the model parameters: this is called the forward problem. The task in the inverse problem is to infer the model parameters given the measurements. This is particularly challenging for an under-determined problem where the number of parameters to estimate is greater than the number of measurements. In such settings, several different source configurations can explain the experimental data and additional constraints are needed to provide a sound solution. In addition, the solution may be highly sensitive to noise in the measurements. Such problems are said to be ill-posed. Note that even over-determined problems can be ill-posed.

The linearity of Maxwell's equations implies that the signals measured by M/EEG sensors are linear combinations of the electromagnetic fields produced by all current generators. The linear operator, called gain matrix in the context of M/EEG, predicts the fields measured by the sensors due to a configuration of sources [START_REF] Mosher | EEG and MEG: Forward solutions for inverse methods[END_REF]. Computing the gain matrix accurately is particularly crucial for EEG, and involves complex numerical solvers [START_REF] Kybic | A common formalism for the integral formulations of the forward EEG problem[END_REF], Gramfort et al. 2010). In the M/EEG literature, solvers known as distributed inverse solvers essentially seek to invert the gain matrix. In practice, the distribution of estimated currents is defined over a discrete set of locations where are positioned current dipoles. The distribution is scalar valued when only their amplitudes are unknown, and vector valued when both amplitudes and orientations of the dipoles need to be estimated. The current generators are commonly assumed to lie on the cortex and are in practice fixed at the locations of the vertices of a cortical mesh [START_REF] Dale | Improved localization of cortical activity by combining EEG and MEG with MRI cortical surface reconstruction: a linear approach[END_REF]. However, the number of generators largely exceeds the number of M/EEG sensors. To tackle this problem, one needs to use a priori knowledge on the characteristics of a realistic source configuration.

The priors most commonly used in the M/EEG community are based on the ℓ 2 norm, leading to what is known as the Minimum Norm (MN) inverse solver [START_REF] Wang | Magnetic source images determined by a lead-field analysis: the unique minimum-norm least-squares estimation[END_REF][START_REF] Dale | Improved localization of cortical activity by combining EEG and MEG with MRI cortical surface reconstruction: a linear approach[END_REF][START_REF] Hämäläinen | Interpreting magnetic fields of the brain: minimum norm estimates[END_REF][START_REF] Pascual-Marqui R D | Low resolution electromagnetic tomography: A new method for localizing electrical activity of the brain[END_REF]. This MN inverse solver leads to a linear solution: i.e., Minimum Norm Estimates (MNE) are obtained by simple matrix multiplication [START_REF] Tikhonov | Solutions of Ill-Posed Problems[END_REF]. This makes the estimation extremely fast. However, ℓ 2 -based solvers suffer from several limitations. Among which is the smearing of the even focal activations, often leading to overestimation of the extents of the activated areas. Also, they require to use a two-step approach where the MNE are post-processed to obtain an interpretable picture of the spatio-temporal activation patterns [START_REF] Pantazis | Spatiotemporal localization of significant activation in MEG using permutation tests[END_REF][START_REF] Gramfort | Tracking cortical activity from M/EEG using graph-cuts with spatiotemporal constraints[END_REF]. To address these limitations many alternatives to MNE have been proposed.

In the mid 90's, Matsuura et al. [START_REF] Matsuura | Selective minimum-norm solution of the biomagnetic inverse problem[END_REF] proposed to regularize the amplitudes of the estimated sources with an ℓ 1 prior using an optimization procedure based on the simplex method. This approach was then later slightly modified by [START_REF] Uutela | Visualization of magnetoencephalographic data using minimum current estimates[END_REF], who called the ℓ 1 penalized solutions minimum-current estimates (MCEs). About the same time, Gorodnitsky et al. proposed to use Iterative Reweighted Least Squares (IRLS) with the FOCUSS algorithm [START_REF] Gorodnitsky | Neuromagnetic source imaging with focuss: a recursive weighted minimum norm algorithm[END_REF] to approximate the solution that would be obtained with an ℓ 0 prior. Subsequently, [START_REF] Phillips | MEG-based imaging of focal neuronal current sources' Medical Imaging[END_REF]) proposed a Bayesian approach based of Markov random fields (MRF) and solved with mean field annealing. All these approaches are motivated by the fact that realistic source configurations are likely to have only a limited number of active sites. For example, only a few brain regions are typically significantly activated by a given cognitive task. The source configuration is said to be spatially sparse. This assumption has proved to be relevant for clinical applications [START_REF] Huppertz | Cortical current density reconstruction of interictal epileptiform activity in temporal lobe epilepsy[END_REF]) and also justifies dipole fitting which is currently the most widely used method in clinical settings.

However, the above approaches suffer from significant limitations. As they promote a sparse solution independently at each time instant, they fail to recover the time courses of cortical sources. In order to go beyond these limitations, there has been a growing interest for methods that promote spatially sparse solutions while taking into the temporal dynamics of the data [START_REF] Phillips | An empirical bayesian solution to the source reconstruction problem in eeg[END_REF][START_REF] Friston | Multiple sparse priors for the M/EEG inverse problem[END_REF][START_REF] Wipf | A unified bayesian framework for MEG/EEG source imaging[END_REF][START_REF] Valdés-Sosa P A | EEG source imaging with spatio-temporal tomographic nonnegative independent component analysis[END_REF][START_REF] Haufe | Combining sparsity and rotational invariance in EEG/MEG source reconstruction[END_REF][START_REF] Ou | A distributed spatio-temporal EEG/MEG inverse solver[END_REF]. While the methods proposed in [START_REF] Phillips | An empirical bayesian solution to the source reconstruction problem in eeg[END_REF][START_REF] Friston | Multiple sparse priors for the M/EEG inverse problem[END_REF][START_REF] Wipf | A unified bayesian framework for MEG/EEG source imaging[END_REF]) are related to sparse Bayesian learning where the problem boils down to the maximization of a non-convex cost function called the model evidence, [START_REF] Haufe | Combining sparsity and rotational invariance in EEG/MEG source reconstruction[END_REF][START_REF] Ou | A distributed spatio-temporal EEG/MEG inverse solver[END_REF]) address the problem using a sparsity-inducing prior that mixes both ℓ 1 and ℓ 2 norms. A ℓ 1 prior is used to promote a spatially sparse solution and a smooth ℓ 2 prior is used either for orientations [START_REF] Haufe | Combining sparsity and rotational invariance in EEG/MEG source reconstruction[END_REF] or both time and orientations [START_REF] Ou | A distributed spatio-temporal EEG/MEG inverse solver[END_REF], leading to a convex optimization problem. A problem is convex when it consists in minimizing a convex function over a convex set [START_REF] Boyd | Convex Optimization[END_REF]. The main reason for the success of these solvers is the structured sparsity induced by the ℓ 21 mixed-norm. Figure 1 illustrates source estimates with a simple ℓ 1 norm compared to a structured prior with a ℓ 21 mixed-norm. The latter leads to a structured sparsity pattern while a simple ℓ 1 norm provides a scattered pattern that is not consistent with what is known about the sources. Here, the ℓ 21 prior guarantees that the active source sites will be the same over the entire time interval of interest. Furthermore, grouping the temporal coefficients with an ℓ 21 norm is a natural way to take into account the smooth temporal dynamics of M/EEG data. More generally, mixed-norm based priors offer a general way to take the structure of a problem into consideration. We call solutions obtained with such priors Mixed-Norm Estimates (MxNE). For an application to other brain imaging methods, see, for example [START_REF] Varoquaux | Advances in Neural Information Processing Systems Advances in Neural Information Processing Systems[END_REF], where a two-level mixed norm was employed for the identification of brain networks using functional Magnetic Resonance Images (MRI) data.

Despite this growing interest, the use of sparsity-inducing priors is still limited to a small group of researchers. One possible reason is that solvers proposed so far are slow when applied to the analysis of real datasets. Another explanation is that algorithms proposed so far are complex and difficult to implement. Indeed, while a basic minimum norm can be computed in a few hundreds of milliseconds, sparse inverse solvers as proposed in [START_REF] Haufe | Combining sparsity and rotational invariance in EEG/MEG source reconstruction[END_REF][START_REF] Ou | A distributed spatio-temporal EEG/MEG inverse solver[END_REF]) can take an hour to converge when realistic dimensions are used. A challenge is therefore to develop efficient optimization strategies that can solve the M/EEG inverse with such priors in a few seconds.

In the last few years, the machine learning and signal processing communities have devoted a lot of efforts into the improvement of the optimization methods that help to solve non-differentiable problems arising when considering sparse priors. One reason is that, under certain conditions, it has been proved that sparsity could enable the perfect resolution of ill-posed problems [START_REF] Donoho | Compressed sensing[END_REF][START_REF] Candès | The dantzig selector: statistical estimation when p is much larger than n[END_REF]. Among the list of algorithms that have been proposed are IRLS methods, similar to the FOCUSS algorithm, that consist in iteratively computing weighted MN solutions with weights updated after each iteration [START_REF] Li | A globally convergent method for lp problems[END_REF][START_REF] Daubechies | Iteratively re-weighted least squares minimization: Proof of faster than linear rate for sparse recovery[END_REF]. The LARS-LASSO algorithm [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF][START_REF] Efron | Least angle regression[END_REF], which is a variant of the homotopy method from Osborne [START_REF] Osborne | A new approach to variable selection in least squares problems[END_REF], is an extremely powerful method for solving the ℓ 1 problem. Simple coordinate descent methods [START_REF] Friedman | Pathwise coordinate optimization[END_REF] or blockwise coordinate descent, also called Block Coordinate Relaxation (BCR) [START_REF] Bruce | Block coordinate relaxation methods for nonparamatric signal denoising[END_REF], are also possible strategies. Alternatively, methods based on projected gradients and proximity operators have been proposed [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF][START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF], Nesterov 2007a[START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. Even if some MxNE can be obtained efficiently, e.g., with coordinate descent, the algorithms proposed in this contribution rely on proximal operators and gradient based methods as they provide a generic approach for all MxNE. They are also grounded on the current mathematical understanding and convergence properties of these solvers.

In this paper, we introduce efficient methods to compute mixed-norm estimates from M/EEG data. The three main contributions of this article are:

(i) We present the M/EEG inverse problem as a convex optimization problem and we explain how structured solutions can be promoted via appropriate priors based on mixed-norms. (ii) We present in detail optimization methods that outperform in terms of convergence speed previously proposed algorithms and derive optimality conditions to control the convergence of the algorithm.

(iii) We then give two examples of MxNE that are relevant for M/EEG using two and three-level mixed-norms, including application to real data.

The first section of the paper provides the necessary background and notation. Mixed-norm estimates with two or three-level mixed-norms are introduced. The second contains the algorithmic and mathematical details of the optimization methods. The third section provides experimental results on real MEG data, demonstrating the efficiency and relevance of the proposed methods.

Mixed-norm estimates (MxNE)

In this section, we introduce inverse problems with linear forward models and more specifically the M/EEG inverse problem. We then define formally the one, two and three-level mixed-norms, explaining their influence on the solutions when used as priors. We explain how a three-level mixed-norm can be used for functional mapping and detail how the ℓ 21 norm can be combined to obtain focal source estimates while promoting smooth time courses.

Framework and notation

Solving an inverse problem consists of estimating one or more unknown signals from observations, typically incomplete and noisy. When considering linear models, the observations, also called measurements, are linear combinations of the signals, also called sources. The linear relationship between the sources and the measurements, of this model, also called the forward model, is commonly derived from the physics of the problem.

Distributed source models in M/EEG use the individual anatomical information derived from MRI [START_REF] Dale | Improved localization of cortical activity by combining EEG and MEG with MRI cortical surface reconstruction: a linear approach[END_REF]. The putative source locations can be then restricted to the brain volume or to the cortical mantle. Due to the linearity of the forward problem, each source adds its contribution independently to the measured signal. We focus here on source models where one dipole with a known orientation is positioned at each location. Source estimates are the amplitudes of the dipoles. Such models are known as fixed orientation. The framework however holds also in the free orientation case where three dipoles share a same spatial location. In this case both amplitudes and orientations need to be estimated.

The measurements M ∈ R N ×T (N number of sensors and T number of time instants) are obtained by multiplying the source amplitudes X ∈ R S×T (S number of dipoles) by a forward operator G ∈ R N ×S , i.e., M = GX. In addition, the measurements are corrupted by an additive noise E:

M = GX + E .
In the context of M/EEG, N lies between 50 for EEG only and 400 for M/EEG combined measurements, while S lies between 5000 and 50000 depending on the precision of the source model considered.

A classical approach to estimate X given M consists in introducing a cost function F whose minimum provides the solution:

X * = arg min X F (X) = arg min X (f 1 (X) + λf 2 (X)) .
(1)

The cost function is composed of two terms:

• A data-fit term, f 1 , that quantifies how well the estimated sources match the measured data. This term takes into account the characteristics of the measurement noise. • A regularization term, a.k.a., penalty term or prior, denoted f 2 , that is used to introduce a priori knowledge on the solution. This term is mandatory to render the solution unique when considering ill-posed problems.

These two terms are balanced by the regularization parameter λ > 0. In the context of M/EEG, f 2 can be directly a function of the source amplitudes or introduce a regularization matrix like a spatial Laplacian D leading to a regularization term of the form f 2 (DX) [START_REF] Pascual-Marqui R D | Low resolution electromagnetic tomography: A new method for localizing electrical activity of the brain[END_REF]). This contribution focuses on cases where f 1 and f 2 are convex functions [START_REF] Boyd | Convex Optimization[END_REF]). As will be detailed later, the convexity of F is a key assumption that allows to obtain globally optimal solutions which are independent of the initialization of the solver. As will be discussed in Section 3, this assumption allows us to employ very efficient optimization procedures whose mathematical properties in terms of complexity and convergence rate are fully understood. Another benefit of convexity observed in practice is the increased stability of the solutions in the presence of noise.

In M/EEG, f 1 (X) is usually the squared ℓ 2 norm of the residual R = M -GX:

f 1 (X) = 1 2 M -GX 2 2 = 1 2 R 2 2 = 1 2 N n=1 T t=1 R 2 n,t . (2)
The smaller is the residual, the better the sources explain the data. The minimization (2) is equivalent to finding a maximum likelihood estimate under the assumption that the additive noise E is Gaussian, i.e., E ∼ N (0, I). In practice, the M/EEG noise is not white but one can estimate the noise-covariance matrix, which can be employed in whitening, either from empty-room data or from periods of actual brain signals void of data of interest [START_REF] Hansen | MEG: An Introduction to Methods[END_REF]. Note that what is called here ℓ 2 norm is in fact the Frobenius matrix norm, since, generally, T > 1. Note also that the factor 1 2 is included for convenience in the derivation of the optimization methods.

The problem therefore reads:

X * = arg min X 1 2 M -GX 2 2 + λf 2 (X) , λ ∈ R + . (3)
The variance of each of the uncorrelated whitened signals is unity, which can help to set the λ parameter. Indeed, assuming that E is Gaussian, the expected value of E 2 2 is N T . It suggests that λ should be chosen such that M -G X 2 2 ≈ N T . This is known as the discrepancy principle [START_REF] Morozov | On the solution of functional equations by the method of regularization[END_REF]). It will be used in Section 4.

We now proceed to discuss suitable priors f 2 for M/EEG.

The ℓ w;p norm

The most common choice for f 2 (X) is the squared ℓ 2 (Frobenius) norm of the sources amplitudes X. This lies in the category of ℓ p norms [START_REF] Matsuura | Selective minimum-norm solution of the biomagnetic inverse problem[END_REF][START_REF] Wagner | Adv Biomagn Res (BIOMAG)[END_REF][START_REF] Uutela | Visualization of magnetoencephalographic data using minimum current estimates[END_REF][START_REF] Gorodnitsky | Neuromagnetic source imaging with focuss: a recursive weighted minimum norm algorithm[END_REF] that work on a time-bytime basis. Throughout this paper, we are often interested in estimating X ∈ R S×T or X ∈ R S×KT when considering K > 1 datasets with a joint estimation. When considering ℓ w;p norms applied to such X, we consider X as a set of coefficients that can be seen as a vector x ∈ R P where P = ST or P = SKT .

Definition 1 (ℓ w;p norm) Let x = (x 1 , . . . , x P ) ∈ R P and w = (w 1 , . . . , w P ) ∈ R P +, some positive weights. Let p ≥ 1. Then the ℓ w;p norm of x is:

x w;p = P s=1 w s |x s | p 1/p
, which is known to be convex for p ≥ 1 and strictly convex for p > 1 [START_REF] Boyd | Convex Optimization[END_REF]).

The reason for introducing weights in the ℓ p norms is due to the fact the columns (G •s ) s of M/EEG forward operators are not normalized. The closer the dipole s is to the head surface, the bigger G •s 2 . This implies that a naïve inverse procedure would favor dipoles close to the head surface. In the M/EEG literature, this is known as the "depth bias" [START_REF] Lin | Assessing and improving the spatial accuracy in MEG source localization by depth-weighted minimumnorm estimates[END_REF]. Using a weighted norm is a way to address this problem.

2.2.1. The ℓ w;2 norm The squared ℓ 2 norm when used both for the data-fit and the penalty term f 2 is known as MNE in the M/EEG literature. The optimization problem reads:

X = arg min X 1 2 M -GX 2 2 + λ 2 X 2 w;2 . ( 4
)
This corresponds to a penalized maximum likelihood estimate assuming the sources are Gaussian and normally distributed, with a diagonal covariance matrix [START_REF] Wipf | A unified bayesian framework for MEG/EEG source imaging[END_REF]. By using such a prior, one spreads the energy of the solution over all the sources. In the context of M/EEG source localization, it leads to activation maps where every brain region has a non-zero amplitude (see Fig. 1-a) and where the extent of active regions is often over-estimated. Solvers based on ℓ w;2 penalty fail to recover high spatial frequencies. In order to avoid this, we can employ a prior that promotes spatially sparse solutions where the data will be explained by a few sources. Keeping f 2 convex, this is can be done with an ℓ 1 norm.

The ℓ w;1 norm

The ℓ 1 norm promotes sparse solutions, which is a strong hypothesis: the solution should only have a small number of non-zero coefficients.

While sparsity can be a valuable assumption in some applications, e.g., denoising [START_REF] Kowalski | Random models for sparse signals expansion on unions of basis with application to audio signals[END_REF][START_REF] Févotte | Sparse linear regression with structured priors and application to denoising of musical audio[END_REF][START_REF] Dupé | A proximal iteration for deconvolving poisson noisy images using sparse representations[END_REF], it can also lead to unrealistic solutions in other applications, e.g., blind source separation [START_REF] Bobin | Advances in Imaging and Electron Physics[END_REF], coding [START_REF] Daudet | International Conference Wavelet analysis and Applications[END_REF], and for M/EEG. Indeed, as illustrated in Fig. 1, an ℓ 1 prior should be used with some caution when performing M/EEG source imaging with temporally correlated data. Such a prior, which estimates the active sources independently at each time instant, will very likely fail to recover the smooth temporal dynamics of a realistic source. To address this limitation a solution recently proposed in the literature estimates the sources for all time instants jointly after introducing a coupling between the estimates [START_REF] Ou | A distributed spatio-temporal EEG/MEG inverse solver[END_REF]. This is achieved using a penalty based on a two-level mixed-norm.

Two-level mixed-norms

In order to define the two-level mixed-norm, we must consider a sequence indexed by a double index (g, m) ∈ N 2 such that (x) = (x g,m ) (g,m)∈N 2 . One can then consider the two canonical subsequences (x g,. ) = (x g,1 , x g,2 , . . .) for a fixed g, and (x .,m ) = (x 1,m , x 2,m , . . .) for a fixed m. This labeling convention introduces a grouping of the coefficients and will be utilized below.

Definition 2 (Two-level mixed-norms) Let x ∈ R P be indexed by a double index (g, m) such that x = (x g,m ).

Let p, q ≥ 1, and w ∈ R P +, * be a sequence of strictly positive weights labeled by a double index (g, m) ∈ N 2 . We call mixed-norm of x ∈ R P , the norm ℓ w;p,q defined by

x w;p,q =   g m w g,m |x g,m | p q/p   1/q
.

Cases p = +∞ and q = +∞ are obtained by replacing the corresponding norm by the supremum.

The two indices g and m can be interpreted as a hierarchy of the coefficients. The double indexing needed by the definition of mixed-norms allows to consider coefficients by groups. Coefficients are indeed distinguished between groups which are blind to each other, and the coefficients that belong to a same group are correlated. With the notation above, the g index can be seen as the "group index" and the m index as the "membership" index. Mixed-norms are then a practical way to induce explicitly a coupling between coefficients, instead of the independence hypothesis behind the ℓ p norms. Hence, mixed-norms allow to promote some structures that have been observed in real signals. Properties of such norms, convexity in particular, enable the use of efficient optimization strategies.

In order to illustrate this, let us consider the use of the ℓ w;21 norm in the problem (3). The ℓ w;21 norm defined on a matrix X ∈ R S×T , and with weights w depending only on the space index s (t indexes time), is given by:

X w;21 = s t w s X 2 s,t .
This corresponds to the sum of the ℓ 2 norm of the lines. As a consequence, an estimation of X given by the minimization of Eq. ( 3) is sparse through the lines, i.e., all the coefficients of a line of X are either jointly nonzero, or all set to zero (see Fig. 1-c). Such a behavior will become more explicit with the definition of the so called proximity operator, see Section 3.2. This approach, proposed earlier for M/EEG [START_REF] Ou | A distributed spatio-temporal EEG/MEG inverse solver[END_REF], avoids the irregular time series obtained with a simple ℓ 1 norm. Note that the general formulation in Definition 2 using (g, m) covers the case with sources having unconstrained orientations. In this case g indexes each spatial location which contains three dipoles [START_REF] Haufe | Combining sparsity and rotational invariance in EEG/MEG source reconstruction[END_REF][START_REF] Ou | A distributed spatio-temporal EEG/MEG inverse solver[END_REF].

The two-level mixed-norms were introduced during the 60's in [START_REF] Benedek | The space l p with mixed norm[END_REF]). These norms were then studied more formally in the context of Besov and modulation spaces [START_REF] Samarah | Local Fourier bases and modulation spaces[END_REF][START_REF] Feichtinger | Modulation spaces: Looking back and ahead[END_REF][START_REF] Rychkov | On restrictions and extensions of the besov and triebel-lizorkin spaces with respect to lipschitz domains[END_REF][START_REF] Grochenig | Non-linear approximation with local fourier bases[END_REF]. Also see [START_REF] Kowalski | Sparse regression using mixed norms[END_REF]) and [START_REF] Kowalski | Sparsity and persistence: mixed norms provide simple signals models with dependent coefficients[END_REF], who introduced the use of the ℓ 12 norm under the name Elitist-Lasso.

Three-level mixed-norms

In this section, we are interested by models where sources X can be indexed by three indices. In the context of M/EEG, these three indices can correspond to the spatial location, the experimental conditions, and the time. For example, for somatosensory data of Section 4, an experimental condition corresponds to the finger that is stimulated. Let us denote this new index by k. The sources, with elements indexed by (s, k, t), are denoted by X ∈ R S×KT (K concatenated datasets) or simply x ∈ R P with P = SKT . Using this notation we can define a three-level mixed norm.:

Definition 3 (Three-level mixed-norms) Let x ∈ R P be indexed by a triple index (s, k, t) such that x = (x s,k,t ). Let p, q, r ≥ 1 and w ∈ R P +, * a sequence of strictly positive weights. We call mixed norm of x the norm ℓ w;p,q,r defined by

x w;pqr =    S s=1   K k=1 T t=1 w s,k,t |x s,k,t | p q/p   r/q    1/r
. Cases p = +∞, q = +∞ and r = +∞ are obtained by replacing the corresponding norm by the supremum.

The inverse problem then reads:

X * = arg min X 1 2 M -GX 2 2 + λ r X r w;pqr , λ ∈ R + . (5)
For our application, we will use the ℓ w;212 mixed-norm. Note that the ℓ 2 , ℓ 12 and ℓ 21 norms are special cases of the latter norm. Indeed ℓ 2 is obtained by setting K = 1, ℓ 21 by setting S = 1 and ℓ 12 by setting T = 1. This suggests that an optimization procedure for the ℓ 212 norm readily works for both ℓ 12 and ℓ 21 norms. Also, it can be observed that ℓ w;221 is equivalent to ℓ w;21 after grouping conditions as well as time instants. By doing so one imposes the active sources to be common between all experimental conditions.

With the ℓ 1 norm to penalize the experimental conditions, while keeping the ℓ 2 norm on other indices, source estimates with non-zero activations for few conditions are promoted. By doing so, one penalizes the overlap between the active regions for the different conditions. With the somatosensory example, such a mixed-norm promotes activations where a given spatial location is active only for one, or at least few, experimental conditions. By definition, this norm corresponds to the a priori information that the stimulation of the different fingers leads to brain activations at different cortical locations, see Section 4.

Algorithms

The algorithms we employ are first-order methods that fit in the same category as the iterative thresholding procedures proposed in [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF]) for the ℓ 1 penalty. We extend them to problems where f 2 is a convex mixednorm [START_REF] Kowalski | Sparse regression using mixed norms[END_REF], which, as explained in Section 2, can take into account the specific characteristics of, e.g., M/EEG source localization. The properties of such algorithms are based on recent mathematical results [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF].

Let us first introduce the notion of proximity operator, a.k.a., proximal operator [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF]):

Definition 4 (Proximity operator) Let φ : R P → R be a proper convex function. The proximity operator associated to φ and λ > 0, denoted by prox λφ : R P → R P reads:

prox λφ (y) = arg min x∈R P 1 2 y -x 2 2 + λφ(x) . (6)

Iterative proximal gradient methods

In a nutshell, proximal gradient methods are a natural extension of gradient-based techniques when the objective function to minimize has an amenable non-smooth part. Such procedures based on iterative thresholding and more generally on projected gradients require that the cost function (1) meets the following hypothesis [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]):

• f 1 is a proper convex function whose gradient is Lipschitz continuous: ∃L ∈ R + such that ∇f 1 (x) -∇f 1 (y) ≤ L xy for all x and y in R P . L is called the Lipschitz constant. • f 2 is a proper convex function (not necessarily differentiable).

In our case, the gradient of the data-fit (2) is Lipschitz continuous. It reads:

∇f 1 (X) = -G T (M -GX) ,
and its Lipschitz constant is given by L = G T G (spectral norm which corresponds to the largest singular value).

The "simplest" iterative scheme to minimize (1) given in Algorithm 1 is called Iterative soft-thresholding and sometimes referred to as Forward-Backward or Landweber iterations [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF].

Algorithm 1: ISTA (Iterative shrinkage/thresholding algorithm) Initialization: Let X (0) ∈ R S×KT (for example 0). repeat

X (k+1) = prox µλ f2 X (k) + µG T (M -GX (k) ) , with 0 < µ < 2 L . until convergence;
The idea is to alternate the minimization over f 1 using a small gradient step and the computation of the proximal operator associated with f 2 . As the proximal operator can be seen as a generalized projection, this algorithm is a generalized iterative projected gradient method (see [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF] for a proof of convergence).

Unfortunately, this algorithm may converge rather slowly. It has been proved that its convergence rate is O(1/k), where k is the number of iterations

∃C > 0, ∀k such that F (X (k) ) -F (X * ) ≤ C k .
To improve the convergence speed, at least two accelerated projected gradient schemes whose convergence speed is O(1/k 2 ) have been proposed (Nesterov 2007b[START_REF] Weiss | Algorithmes rapides d'optimisation convexe[END_REF][START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. The FISTA (Fast Iterative shrinkage/thresholding algorithm) algorithm [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]) is one of them. It's a small modification of ISTA that takes into account the previous descent direction. It is a two-step approach. Note that a classical example of a multi-step approach is the conjugate gradient algorithm used to solve positive definite linear systems. More details on these approaches can be found in [START_REF] Tseng | Approximation accuracy, gradient methods, and error bound for structured convex optimization[END_REF].

Algorithm 2: FISTA Initialization:

X (0) ∈ R S×KT , Z (1) = X (0) , τ (1) = 1, k = 1, 0 < µ < 1 L repeat X (k) = prox µλf2 Z (k) + µG T (M -GZ (k) ) τ (k+1) = 1+ √ 1+4τ (k) 2 2 Z (k+1) = X (k) + τ (k) -1 τ (k+1) (X (k) -X (k-1) ) until convergence;
In order to tackle the optimization problem (3), one just needs to know how to compute prox µλf2 where f 2 is a mixed-norm presented in Section 1.

Proximity operators corresponding to mixed-norms

The following proposition gives details of the proximity operators associated with the mixed-norms presented in Section 1. It corresponds to the solutions of (6) when φ is a mixed-norm.

Proposition 1 (Proximal operators for MxNE) Let x ∈ R P and y ∈ R P . Let w ∈ R * P + be a vector of weights.

ℓ 2 norm Let x be indexed by s. The proximity operator associated to the squared ℓ 2 norm is given by x = prox λ . 2 w,2 (y) where x reads coordinate by coordinate:

x s = y s 1 + λw s .
ℓ 1 norm Let x be indexed by s. The proximity operator associated to the ℓ 1 norm is given by x = prox λ . w,1 (y) where x reads coordinate by coordinate:

x s = y s |y s | (|y s | -λw s ) + .
The function (•) + is defined as (a) + = max(a, 0) and we use the convention 0/0 = 0. The proximity operator for the ℓ 1 norm is known as "soft-thresholding".

ℓ 21 norm Let x be indexed by (s, t). Let us consider a vector of weights used to weight each group. The proximity operator associated to the ℓ 21 norm is given by x = prox λ . w,21 (y) where x reads for each coordinate:

x s,t = y s,t 1 - λ √ w s y s 2 + ,
where y s is the vector formed by the coefficients indexed by s.

ℓ 12 norm Let x be indexed by (s, k). Let us consider w a vector of weights used to weight each group. Let r s,k be defined such that r s,k def = y s,k /w s,k . For each s, let the indexing denoted by k ′ s be defined such that ∀k ′ s , r s,k ′ s +1 ≤ r s,k ′ s . Let the index K s be such that:

λ Ks k ′ s =1 w 2 s,k ′ s r s,k ′ s -r s,Ks < r s,Ks ≤ λ Ks+1 k ′ s =1 w 2 s,k ′ s r s,k ′ s -r s,Ks
The proximal operator x = prox λ . w,12 (y) is given coordinate by coordinate:

x s,k = y s,k |y s,k |   |y s,k | - λ 1 + λK ws Ks k ′ s =1 y s,k ′ s   + ,
where

K ws = Ks k ′ s =1 w 2 s,k ′ s .
This proposition shows the effect of such proximal operators on their inputs. For ℓ 2 , it is a simple weighting. The associated cost function being differentiable, a proof is obtained simply by computing the derivative with respect to each x s . For ℓ 1 it is a thresholding of all the coefficients independently (see [START_REF] Donoho | De-noising by soft-thresholding[END_REF] for a proof). As a result, some coefficients are set to zero which reflects the sparsity obtained with such a penalty. With the ℓ 21 norm, a group is globally set to zero depending on its norm. A coefficient is non-zero only if the norm of the group it belongs to is large enough. If groups are formed by rows then the ℓ 21 prior promotes a row structured sparsity pattern as illustrated in Fig. 1-c. For completeness, a derivation of this proximal operator is given in Appendix A. Note that such results have been previously obtained like in [START_REF] Kowalski | Sparse regression using mixed norms[END_REF]). However, the later work does not address the weighted case.

From Proposition 1, the proximal operators associated to any mixed-norm combining ℓ 1 and ℓ 2 norms can be derived. This is in particular the case of the ℓ 212 norm for which the proximal operator is given by the following proposition.

Proposition 2 (Proximal operator associated to the ℓ w;212 norm) Let y ∈ R P be indexed by (s, k, t). Let w ∈ R P be a vector of positive weights such that ∀t, w s,k,t = w s,k . Let us define [y s,k ] def = w s,k t y 2 s,k,t and r s,k def = [y s,k ]/w s,k . For each s, let the indexing denoted by k ′ s be defined such that ∀k ′ s , r s,k ′ s +1 ≤ r s,k ′ s . Let the index K s be such that:

λ Ks k ′ s =1 w s,k ′ s r s,k ′ s -r s,Ks < r s,Ks ≤ λ Ks+1 k ′ s =1 w s,k ′ s r s,k ′ s -r s,Ks .
Then, x = prox λ . 2 w;212 (y) is given, for each coordinate (s, k, t), by:

x s,k,t = y s,k,t      1 - λ √ w s,k 1 + λK ws Ks k ′ s =1 [y s,k ′ s ] y s,k 2      + ,
where K ws = Ks ks=1 w s,ks . A proof of this proposition is given in Appendix B. Having established the minimization procedures for MxNE, we need next to test for the convergence and the optimality of the current iterate X (k) .

Optimality conditions and stopping criterion

When the cost function F to be optimized is smooth, a natural optimality criterion is obtained by checking whether the gradient is small: ∇F (X (k) ) < ε. Unfortunately, this approach does not apply to the non-differentiable cost-functions involving ℓ 1 norms.

An answer for convex problems more generally consists of computing, if possible, a duality gap. For a subset of these problems the Slater's conditions apply and, consequently, the gap at the optimum proves to be zero [START_REF] Boyd | Convex Optimization[END_REF]. Computing the gap starts by deriving a dual formulation of the original problem, also called the primal problem. The duality gap is defined as the difference between the minimum of the primal cost function F p and the maximum of the dual cost function F d . For a value of X (k) of the primal variable at iteration k, if one can exhibit a dual variable Y (k) , the duality gap η (k) is defined as:

η (k) = F p (X (k) ) -F d (Y (k) ) ≥ 0
At the optimum (corresponding to X * ), if the Y (k) associated to X (k) is properly chosen, η (k) is 0. By exhibiting a pair (X (k) , Y (k) ) one can guarantee that:

F p (X (k) ) -F p (X * ) ≤ F p (X (k) ) -F d (Y (k) ) .
A good stopping criteria is therefore given by η (k) < ε. The solution meeting this condition is said to be ε-optimal. The challenge in practice is to find an expression for F d and to be able to associate a "good" Y to a given X for the problems of the form (3).

We now give a general answer for the class of problems detailed in this contribution.

The solution is derived from the Fenchel-Rockafellar duality theorem [START_REF] Rockafellar | Convex Analysis[END_REF]) which leads to the following dual cost function:

F d (Y ) = - 1 2 Y 2 2 + Tr(Y T M ) -λf * 2 (G T Y /λ) (7)
where Tr stands for the trace of a matrix and f * 2 is the Fenchel conjugate of f 2 defined by:

f * 2 (Z) def = sup X∈R P Tr(Z T X) -f 2 (X) .
In Appendix C we provide the Fenchel-Rockafellar duality theorem in order to obtain this result.

The Fenchel conjugates of mixed-norms and squared mixed-norms, which remain to be given in (7), can be computed thanks to the following proposition: Proposition 3 (Fenchel conjugate of a mixed-norm) (i) The Fenchel conjugate of norm u p1,...,pn is the indicator function of the dual norm:

v → i v * p 1 ,...,pn = 0 if v p ′ 1 ,...,p ′ n ≤ 1 +∞ if not where ∀i, p ′ i is such that 1 pi + 1 p ′ i = 1
(ii) The Fenchel conjugate of the function u → 1 2 u 2 p1,...,pn is the function:

v → 1 2 v 2 p ′ 1 ,...,p ′ n
Moreover, the Karush-Khun-Tucker (KKT) conditions of the Fenchel-Rockafellar duality theorem (see Appendix C) give a natural mapping from the primal space to the dual space:

Y (k) = M -GX (k) ,
which then needs to be modified in order to satisfy the constraint of f * 2 . When f * 2 is an indicator function, it corresponds to a projection on the associated convex set. It is then possible to use the dual gap as stopping criterion. Algorithm 3 summarizes the computation of the dual gap, in cases of ℓ 21 and ℓ 212 penalty function. An 1D illustration with the ℓ 1 is provided in Fig. 2.

Algorithm 3: Duality gap for ℓ 212 or ℓ 21 Entries: X (k) Mapping to the dual space:

Y (k) = M -GX (k) Compute f * 2 (G T Y (k) ): if f 2 = ℓ 21 then
Project dual variable on the constraint if necessary:

Y (k) = Y (k) / max( G T Y (k) 2∞ /λ, 1) f * 2 (G T Y (k) /λ) = 0 if f 2 = ℓ 212 then f * 2 (G T Y (k) /λ) = 1 2 G T Y (k) /λ 2 2,∞,2
Compute duality gap: From a numerical point of view, every solution is ε-optimal for a particular value of ε. The duality gap observed at the end of the computation is for example limited by machine precision. Also, Algorithm 3 shows that η (k) depends on the scaling of the data. Therefore, in practice the duality gap is meaningful if the input data have been scaled or normalized in a certain way. This is guaranteed with M/EEG data by the pre-whitening step. Our experimental results show that for whitened data a duality gap lower than 10 -5 does not produce distinguishable solutions.

η (k) = 1 2 M -GX k 2 2 + λf 2 (X k ) + 1 2 Y (k) 2 2 -Tr(Y T M ) + λf * 2 (G T Y /λ)
3.4. An active set method to improve the convergence speed using the ℓ 21 norm Like the ℓ 1 norm, the ℓ 21 norm leads to sparse solutions; only a few sources will have non-zero activations. Knowing this, we can accelerate the computation with an active set strategy which will restrict computations to sources likely to have non-zero activations. This amounts to solving problems with only a subset of columns of G. Let us call Γ the set of sources considered in the sub-problem and X * Γ the associated estimated sources. By computing the duality gap associated to the full problem for a value of X, such that X restricted to Γ is equal to X * Γ and where the rest of X is filled with 0, one can test the quality of the solution for the full problem. If this solution is not good enough according to the stopping criteria, one needs to add to Γ sources that are likely to be active. Such sources violate the KKT optimality conditions [START_REF] Boyd | Convex Optimization[END_REF]). These conditions are specific to the penalty considered. With a ℓ w;21 penalty, denoting by W the diagonal matrix of weights, the KKT optimality conditions impose the following constraint on X:

W -1 G T (M -GX) 2∞ ≤ λ . ( 8
)
The indices of the sources that need to be added to the active set at a next iteration, are given by the indices of the rows of W -1 G T (M -GX) that do not meet this constraint. Intuitively it says that the sources that should be added to the active set are the ones whose forward fields correlate the most with the current residual. Such an active set strategy is known as forward as the size of the problem keeps increasing.

In practice, one can simply add to Γ the source that violates the most the latter constraint. This can however be rather slow if the active set contains hundreds of variables. That would mean running FISTA hundreds of times. A natural idea consists in adding groups of sources, i.e., the set of sources that violate the most the constraint. When no more source violates the KKT constraint, the solver has converged to an optimal solution. The number of sources that should be added to the active set at each iteration is however application specific. For M/EEG, we have found that adding blocks of 10 new variables is a good trade-off. For an optimal solution containing at the most about a hundred active sources a solution is obtained in practice by running the solver ten times at the most on very small problems. Note that the procedure do guarantee the optimality of the solution at the end as the active set can only grow meaning that the solver will end up solving the original full problem [START_REF] Roth | ICML '08: Proceedings of the 25th international conference on Machine learning[END_REF].

Using the active set strategy, the solution corresponding to an experimental data set (with about 300 channels, 200 time samples and 10000 sources) can be obtained in a few seconds. This means that the proposed optimization strategy makes the use of the ℓ 21 penalty computationally trackable in practical M/EEG applications, which was not the case using the methods relying on second order cone solvers proposed in [START_REF] Haufe | Combining sparsity and rotational invariance in EEG/MEG source reconstruction[END_REF][START_REF] Ou | A distributed spatio-temporal EEG/MEG inverse solver[END_REF]).

Convergence results on simulated data

In order to illustrate the convergence rate of the algorithms detailed above in a realistic experimental setting, we performed a simulation using a real MEG gain matrix (151 sensors and 5,000 sources). The implementation used is written in Matlab and involves only linear algebra and standard vectorized operations. Fig. 3 shows the size of the duality gap as a function of computation time using ISTA, FISTA, or FISTA with the active set (AS-FISTA) approach on a problem with an ℓ 21 prior. One can observe that FISTA actually converges much faster than ISTA and that AS-FISTA clearly outperforms both of them. For comparison we also ran on the same configuration of dipoles a SOCP (Second Order Cone Program) solver using the CVX Matlab toolbox (http://cvxr.com/cvx/) as in [START_REF] Ou | A distributed spatio-temporal EEG/MEG inverse solver[END_REF]). Employing only 1,000 dipoles out of 5,000, the computation took 86 s. Such a solver relies of the inversion of a Hessian matrix whose size in S × S. It's cost per iteration is therefore O(S 3 ), i.e., cubic in the number of dipoles, and it also requires to store in memory of matrix of size S × S which may be prohibitive. To tackle the realistic problem used for Fig. 3, it suggests that besides the problem of storing a 5000 × 5000 matrix in memory, computation should approximately be multiplied by 125 which corresponds to more than 2 hours of computation.

Figure 3. Convergence of ISTA, FISTA and AS-FISTA using a . w;21 penalty with a real MEG lead field (151 sensors and 5,000 sources) and synthetic measurements. It can be observed that ISTA can be slow to converge compared to FISTA and that the active set strategy speeds significantly the convergence.

Simulations and MEG results

The following section first presents results with the ℓ 21 norm applied to M/EEG data, and then some simulation and experimental results obtained with the ℓ 212 prior. We show that our solver applied with an ℓ 21 norm provides accurate results in a few seconds on a real auditory M/EEG dataset and that the ℓ 212 norm can improve the accuracy of reconstructions when performing functional mapping of the somatosensory cortex.

MxNE with the ℓ 21 norm

The data used to illustrate the performance of the ℓ 21 MxNE consists of MEG and EEG combined recordings of the evoked response to left-ear auditory pure-tone stimulus ‡. Data were acquired using a 306-channel MEG Neuromag Vectorview system with 60 EEG electrodes, sampled at 600 samples/s. The signals were low-pass filtered at 40 Hz. The noise covariance matrix was estimated from the 200 ms of recordings preceding each stimulus. The source space consisted of 8192 dipoles on the cortex with orientations constrained to be normal to the cortical mantle. Two channels with technical artifacts were ignored from the analysis resulting in a lead field matrix of size 364 × 8192. We averaged 55 epochs and the sources were estimated between 0 and 400 ms resulting in 241 time samples, hence M ∈ R 364×241 .

The results are presented in Fig. 4. The computation time using AS-FISTA for the entire source estimation was 20 s on a laptop (4 GB of RAM and 2.8 GHz CPU). At the optimum the cost function values was 2073.2 and the duality gap about 10 -5 corresponding to a change in fifth significant digit of the cost function. The estimated sources are located in both contralateral and ipsilateral primary auditory cortices (cA and iA). A first early component is observed in cA between 30 and 50 ms with a peak around 90 s in cA and later at 100 ms in iA.

In many respects the source estimates look similar to standard multi-dipole fittings results. However, a few remarks should be made about the present results. First, when working with constrained orientations, the signs of the estimated wave forms are dependent on the orientation used for the dipole, i.e., the wall of a sulcus on which the dipoles are located. This is a fundamental problem of M/EEG source imaging well known from the classical MNE. Also, the cluster of 3 active dipoles in iA illustrates a natural behavior of convex priors. These 3 dipoles have very similar forward fields making them very hard to disambiguate with M/EEG. The stability of the ℓ 21 convex prior produces this cluster of dipoles while a non-convex prior, e.g., ℓ 2p with p < 1, would certainly pick any one of these dipoles. 1 overlaps with the active region of experimental condition 2, which may in practice be unrealistic considering the underlying physiology. The primary somatosensory cortex (S1) [START_REF] Penfield | The Cerebral Cortex of Man: A Clinical Study of Localization of Function[END_REF] and the primary visual cortex (V1) [START_REF] Wandell | Visual field maps in human cortex[END_REF]) are examples of brain regions with such known organizations. The different body parts are mapped to locations at S1 and a location in the visual field maps to an area at V1 This later is known as the retinotopic organization of V1 [START_REF] Wandell | Visual field maps in human cortex[END_REF]. Recent work, such as [START_REF] Sharon | The advantage of combining MEG and EEG: Comparison to fMRI in focally stimulated visual cortex[END_REF][START_REF] Hagler | Source estimates for MEG/EEG visual evoked responses constrained by multiple, retinotopically-mapped stimulus locations[END_REF], emphasize the interest for advanced methods able to reveal such functional organizations noninvasively. However, while [START_REF] Hagler | Source estimates for MEG/EEG visual evoked responses constrained by multiple, retinotopically-mapped stimulus locations[END_REF] propose to use functional MRI data to improve MNE results, we propose the ℓ 212 prior which is a principled way of taking into account the spatial properties of such regions in order to obtain better functional mappings of the brain using only M/EEG data.

To illustrate the fact that an ℓ 2 prior tends to smear the activations and therefore to overestimate the extent of the active regions, we first performed a simulation study.

4.2.2.

Simulation study We generated a synthetic dataset that mimics the organization of the primary somatosensory cortex (S1) [START_REF] Penfield | The Cerebral Cortex of Man: A Clinical Study of Localization of Function[END_REF]. Three non-overlapping cortical regions with a similar area (cf. Fig. 5a), that could correspond to the localization of three right hand fingers were assumed and were used to generate synthetic measurements corrupted with an additive Gaussian noise. The amplitude of activation for the most lateral region (colored in red in Fig. 5), that could correspond to the thumb, was set to be two times as big as the amplitudes of the two other regions. An inverse source estimate was then computed with a standard ℓ w;2 (4) norm, an ℓ w;1 norm, and the ℓ w;212 mixed-norm (5). Each source in the three simulated active regions was then assigned a label corresponding to the condition for which its estimated amplitude was the largest. Quantification of performance was done for multiple values of signal-to-noise ratio (SNR) by counting the percentage of dipoles that have been incorrectly labeled. The SNR is defined here as 20 times the log of the ratio between the norm of the signal and the norm of the added noise. The results are also presented in Fig. 6. It can be observed that the ℓ w;212 always produces the best result and that the ℓ 1 norm is more strongly affected by the decrease in SNR. In order to have a fair comparison between all methods, the λ was set in each case to have M -GX * 2 equal to the norm of the added noise, always known in simulations. The depth bias was corrected in the ℓ w;212 norm case by setting w s,k = w s = G •s 2 2 . This amounts to scaling the columns of G. The depth bias correction was applied the same way for ℓ 1 and ℓ 2 solutions.

The results are illustrated in Fig. 5b and 5c on a region of interest (ROI) around the left primary somatosensory cortex. It can be observed that in the ℓ w;2 norm result the extend of the most lateral region is overestimated while the result obtained with the ℓ w;212 mixed-norm is clearly more accurate.

MEG data

We also analyzed somatosensory data acquired using a CTF Systems Inc. Omega 151 system at a 1250 Hz sampling rate. The somatosensory stimulus was an electrical square-wave pulse delivered randomly to the thumb, index, middle, and little finger of the right hand of a healthy right-handed subject. The evoked response was computed by averaging 400 repetitions of the stimulation of each finger. The triangulation over which cortical activations have been estimated contained a high number of vertices (about 55,000). The forward solution was computed using the spherically symmetric head model [START_REF] Sarvas | Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem[END_REF]. The sphere was positioned to match the shape on the inner surface of the skull near the primary somatosensory cortex (central sulcus). Prior to source estimation, data were whitened using a noise covariance matrix estimated during baseline periods.

The source estimates during the period between 42 and 46 ms are presented in Fig. 7. For each type of prior, the regularization parameter λ was set in order for X * to satisfy M -GX * 2 2 ≈ N KT , knowing that after whitening the data, N KT is a good estimate of the noise variance.

During the time interval of interest the measured magnetic fields indicate the currents are directed into the cortex at S1 which is known to lie on the posterior bank of the central sulcus. Therefore, we next ignored the positive activations located on the anterior bank. Each dipole with negative amplitude was then assigned a label between 1 and 4 based on its maximum amplitude across the 4 conditions. For each condition, equivalently each label, the biggest connected component of dipoles with the same label was kept. Each of the 4 estimated components, corresponding to the 4 right hand fingers are presented in Fig. 7. Solutions using both ℓ w;2 and ℓ w;212 norms are presented. The solution with ℓ w;1 is not represented as the sparsity promoted do not produce a continuous mapping which make results difficult to compare.

With ℓ w;212 the well known organization of the primary somatosensory cortex [START_REF] Penfield | The Cerebral Cortex of Man: A Clinical Study of Localization of Function[END_REF]) is successfully recovered with regions of similar size for each finger. While with ℓ w;2 , the component corresponding to the index finger is overestimated leading to an incorrect localization of the area corresponding to the thumb. Note that some activation does remain at the right location for the thumb using the ℓ w;2 norm. However, it is not strong enough to match with the biggest connected component represented here. These results demonstrate that an alternative to the standard ℓ 2 priors, can improve the localization of cortical activations by offering the possibility to use a prior between different conditions. By solving the inverse problem of multiple conditions simultaneously and by using a mixed-norm that sets an ℓ 1 prior between each condition, our method penalizes current estimates with an overlap between the corresponding active regions. When such a hypothesis holds, the localization of the neural activity becomes more accurate with increasing number of conditions recorded and included in the analysis.

Discussion

In the present paper, we capitalize on advanced numerical methods to tackle multiple convex optimization problems present in many applications such as functional brain imaging using M/EEG. Our paper provides a unifying view of many solvers previously proposed in the M/EEG literature and is to our knowledge the first demonstration that the M/EEG inverse problem can be solved in a few seconds using non-ℓ 2 priors. Rapid computations are essential in functional brain imaging, since interactive exploratory analysis is often needed.

This work relates to the distributed solvers based on sparse Bayesian regression that have been recently proposed [START_REF] Nummenmaa | Automatic relevance determination based hierarchical Bayesian MEG inversion in practice[END_REF][START_REF] Friston | Multiple sparse priors for the M/EEG inverse problem[END_REF][START_REF] Wipf | A unified bayesian framework for MEG/EEG source imaging[END_REF]). These solvers are not explicitly derived from cost-functions like (3) and they lead to non-convex optimization problems not covered by the algorithms detailed above. Note also that formulating the inverse problem as the minimization of a cost function does not guarantee convexity. For example, Valdes-Sosa et al. propose in (Valdés-Sosa et al. 2009) to estimate the source activations as a linear combination of a small number of spatiotemporal maps. Here again, sparsity is a key assumption of the method, however, the minimized cost function is not jointly convex in space and time. The consequence of the non-convexity for all these methods, is that the optimality of the solutions cannot be guaranteed and that solutions depend on the initialization of the algorithm. Our formulation of MxNE does not suffer from these shortcomings.

Another benefit of MxNE is the diversity of a priori knowledge that can be taken into account. With the same mathematical foundations and very similar implementations, the ℓ 21 norm can be used to promote sparse source estimates with smooth temporal activations, an ℓ 221 norm can furthermore impose a common set of active dipoles between conditions, and the ℓ 212 norm can be used for more accurate functional mapping. From the neuroscience perspective, the ℓ 21 prior models the a priori assumption that active brain regions should be consistent during a time interval. This assumption is adapted to some datasets like the auditory data presented here, however it is likely to be wrong for a long time interval during which active sources are changing. The ℓ 212 prior is motivated by its ability to explicitly model the functional specificity of brain regions thus leading to more insights on neural circuitry [START_REF] Chklovskii | Maps in the brain: what can we learn from them?[END_REF]. The latter application is also to our knowledge the first attempt to improve the M/EEG inverse problem by using multiple datasets jointly.

Finally, an important point is that there is no prior that fits all needs. That is why MxNE does not refer to a particular prior but to a class of solvers that use mixednorms to better constrain the M/EEG inverse problem. This implies that depending on the assumptions made about the underlying sources, a particular mixed-norm can be more relevant than others.

Conclusion

In this article, we have shown how mixed-norms can be used to promote structure on inverse estimates in order to take into account some a priori knowledge. In the case of M/EEG, the a priori knowledge is based on the understanding of both neurophysiology and biophysics.

This contribution provides principled first-order optimization strategies which are simple to implement and fast, especially when an active set approach is used. Furthermore, the speed of convergence of these algorithms is well understood thanks to a theoretical analysis. All these algorithms rely on proximity operators, which are in practice special thresholding operators. Moreover, we were able to construct practical criteria to stop the optimization process while guaranteeing the optimality of the solutions obtained.

The utility of mixed-norms is demonstrated with both synthetic and experimental MEG data. Our results match existing knowledge about auditory evoked responses and lead to a more accurate mapping of the somatosensory homunculus than the unstructured standard methods. We can apply this Theorem to the functional (3) with g(X) = -1 2 M -X 2 2 and f (X) = λf 2 (X). Then, one just have to compute the conjugate of f . By definition of the dual, given here for a concave function, we have

g * (Y ) def = inf X Tr(Y T X) + 1 2 M -X 2 2 = -sup X -Tr(Y T X) - 1 2 M -X 2 2
then, by the change of variable Y = M -X we have

g * (Y ) = Tr(Y T M ) -sup Z Tr(Y T Z) - 1 2 Z 2 2 .
Moreover, we know by Proposition 3 that the Fenchel conjugate of a squared norm is the squared dual norm. Then, we have

g * (Y ) = Tr(Y T M ) - 1 2 Y 2 2 .
For the Fenchel conjugate of f (X) = λf 2 (X), one can apply simple calculus rules given in [START_REF] Boyd | Convex Optimization[END_REF], which give f * (Y ) = λf * 2 (Y /λ). Finally, the Fenchel-Rockafellar dual function of (3) is given by

Tr(Y T M ) - 1 2 Y 2 2 -λf * 2 (G T Y /λ) .
Furthermore, one can check that the dual variable Y = M -GX verifies the KKT condition

- 1 2 M -GX 2 2 + Tr(Y T M ) - 1 2 Y 2 2 = GX, Y .

NonFigure 1 .

 1 Figure 1. (a), (b) and (c) show in white the non-zero in the estimated source amplitudes obtained with the three norms. The non-zero coefficients are shown in white. While ℓ 2 yields only non-zero coefficients (all sources have a non-zero amplitude), ℓ 21 promotes non-zero coefficients with a row structure (only a few sources have non-zero amplitude over the entire time interval of interest).

Figure 2 .

 2 Figure 2. Duality gap illustration with a 1D ℓ 1 problem (λ = 2 and M = 2). It can be observed that the minimum of Fp is equal to the maximum of F d , i.e., that the duality gap is 0.

Figure 4 .

 4 Figure 4. ℓ w;21 estimates on auditory M/EEG data. Estimation leads to 4 active dipoles in both contralateral and ipsilateral auditory cortices (cA and iA).

4. 2 .

 2 MxNE with the ℓ 212 prior: Functional mapping 4.2.1. Motivation During an M/EEG experiment, a subject is generally asked to perform different cognitive tasks or to respond to various stimuli. Without an adapted prior, it may occur that the estimated active cortical region in experimental condition ‡ Condition 1 of the sample data provided by the MNE software.

Figure 5 .

 5 Figure 5. Simulation results on the primary somatosensory cortex (S1) (SNR = 20dB). Neighboring active regions reproduce the organization of S1. It illustrates that the ℓ w;212 prior improves over a simple ℓ w;2 prior.

Figure 6 .

 6 Figure 6.Evaluation of ℓ w;2 , ℓ w;212 and ℓ w;1 estimates on synthetic somatosensory data. The error represents the percentage of wrongly labeled dipoles.

Figure 7 .

 7 Figure 7. Labeling results of the left primary somatosensory (S1) cortex in MEG using both ℓ w;2 and ℓ w;212 priors. Source estimated during the period between 42 and 46 ms. The ℓ w;212 leads to a more coherent estimate of the functional organization of S1.

.

  Note that this result gives also the proof of the proximal operator associated to the Elitist-Lasso in Proposition 1.Appendix C. Proof of equation (7) Theorem 1 (Fenchel-Rockafellar duality[START_REF] Rockafellar | Convex Analysis[END_REF])) Let f : R M ∪ {+∞} → R be a convex function and g : R N ∪ {+∞} → R a concave function. Let G be a linear operator mapping vectors ofR M to R N . Then inf X∈R M {f (X)g(GX)} = sup Y ∈R N {g * (Y )f * (G T Y )}where f * (resp. g * ) is the Fenchel conjugate associated to f (resp. g), and G T the adjoint operator of G.Moreover, the Karush-Kuhn-Tucker (KKT) conditions read:f (X) + f * (G T u) = X, G T Y , g(GX) + g * (Y ) = GX, Y .

  We can differentiate the functional, when x k w;2 = 0, to obtain the variational system:|x k,t | = |y k,t |λ √ w k |x k,t | x k Appendix B. Proof of Proposition 2 (ℓ 212 proximal operator)We can differentiate the functional, when x s,k,t = 0, to obtain the variational system:|x s,k,t | = |y s,k,t |λ √ w s,k |x s,k,t | x s,k s w;21 = y s,k 2 ⇒ x s,k 2 + λ √ w s,k x s w;21 = y s,k 2 . (B.3)This last equality is true for all k. Then for k and l satisfying x s,k 2 > 0 and x s,l 2 > 0, we have:⇔ x s,k 2 = y s,k 2 -λK ws ( x s,k 2y s,k 2 )λ √ w s,k l/ x s,l 2 >0 √ w s,l y s,l 2 ⇔ x s,k 2 = y s,k 2 -, where K ws = k/ x s,k 2>0 w s,k . We therefore have for all s: If we reorder the √ w s,k y s,k 2 and introduce K s defined in the proposition, it leads to: | x s,k 2 x s,k 2 + λ √ w s,k x s w;21 Using (B.3), we get: |x s,k,t | = |y s,k,t | x s,k 2 y s,k 2 . By injecting the result (B.6) in this equation we get: |x * s,k,t | = |y s,k,t | y s,k 2 -

	1 2 y s,k 2 -x w;21 = Proof We are looking for: with x * = arg min x x s,k 2 = with [y s ] = Ks k=1 √ w s,k y s,k 2 . Let us rewrite (B.2): y -x 2 2 + λ x w;21 , k t w k x 2 k,t . y -x 2 2 + λ x 2 λ + √ w s,k 1 + λK ws [y s ] . w;212 , with x 2 w;212 = s   k t w s,k x 2 s,k,t 1/2   2 -1 . 2 sign(x s,k,t ) = sign(y s,k,t ) which gives: |x s,k,t | 1 + λ √ w s,k x s,k -1 2 By summing over t, we get: ⇒ x 2 s,k,t 1 + λ √ w s,k x s,k -1 2 x s w;21 s,k,t 2 = y 2 x s w;21 = |y s,k,t | |x s,k,t | = |y s,k,t | 1 + λ √ w s,k x s,k -1 2 x s w;21 = |y s,k,t λ x s w;21 √ w s,k 1+λKw s [y s ] + y s,k 2 λ √ w s,k = |y s,k,t | 1 -1 + λK ws	(B.6) (B.1) (B.2)
	x s,k 2 1 + λ √ w s,k x s,k	-1 2
									-1 2
	x s,k 2 √ w s,k	=	x s,l 2 √ w s,l	+	y s,k 2 √ w s,k	-	y s,l 2 √ w s,l	.	(B.4)
	Furthermore, we have that:
	x s w;21 =
						 +
						
					ν x 2 k,ν
								+
	|x k,t | = = = |y k,t | 1 -   |y k,t | -  |y k,t | -λ λ √ w s,k 1+λKw s √ w k |x k,t | λ √ w k |x k,t | ν y 2 k,ν x 2 k,t y 2 k,t  + √ w s,l y s,l 2    l/ x s,l 2>0  |x k,t | |y k,t | y k 2 λ √ w k y k 2 + , x s,k 2 =   y s,k 2 -λ √ w s,k 1 + λK ws k/ x s,k 2 >0	y s,l 2 √ w s,l -w s,k y s,k 2 √ w s,k y s,k 2  
	which is the desired result.

k

x k w;2 = sign(x k,t ) = sign(y k,t ) which gives:

|x k,t | 1 + λ √ w k x k -1 2 = |y k,t | As 1 + λ √ w k x k -1

2 does not depend on t, it implies that for all t and ν:

|x k,t | = |y k,t | |y k,ν | |x k,ν | .

By injecting this last result in Eq. (A.1) we obtain

|x k,t | =   |y k,t | -λ √ w k |x k,t | x l/ x s,l 2>0 √ w s,l x s,l 2 .

which implies that:

x s,k 2 = y s,k 2λ √ w s,k l/ x s,l 2>0 √ w s,l x s,l 2 . (B.5)

By injecting (B.4) in (B.5) we get:

x s,k 2 = y s,k 2λ √ w s,k l/ x s,l 2 >0 w s,l

x s,k 2 √ w s,k + + .

Acknowledgments

The authors acknowledge support from the ANR ViMAGINE ANR-08-BLAN-0250-02, the National Center for Research Resources P41 RR014075-11, and the NIH National Institute of Biomedical Imaging and Bioengineering 5R01EB009048. The authors also wish to thank the anonymous reviewers for their careful comments which helped improve the presentation and the readability of this article.

Appendix A. Derivation of the (ℓ 21 proximal operator)

We are looking for:

x * = arg min x


Gil Moretto

Nicolas Epchtein

Maud Langlois

Isabelle Vauglin

An Off-Axis Telescope Concept for Antarctic Astronomy

Keywords: infrared cases for Antarctica, off-axis telescope, low scattered light, low emissivity

The Antarctic plateau offers exceptional atmospheric and environmental conditions for astronomical observations over a wide range of wavelengths and uniquely favorable to infrared astronomy. Exceptional low sky brightness throughout the near-and mid-infrared and a telescope facility complying with the highest possible dynamic range for photometry, angular resolution and the wide-field leads to the possibility of a modest-sized 2m off-axis telescope achieving comparable sensitivity to that of a larger ground-based 8-10m class telescope or a same sized space-based ones.

INTRODUCTION

We address here the scientific relevancy and feasibility of creating an astronomical observatory for the optical/infrared spectral range on the Antarctic Plateau/Dome C. The properties of the atmosphere above the Antarctic Plateau are known to be unique on the Earth surface. Atmospheric turbulence is concentrated in a thin layer of a few tens of meters, the sky opacity, particularly in the infrared, is considerably reduced and the thermal infrared sky background radiation is lower by a factor of 10 to 20 in the 2-3 µm windows. The combination of these advantages makes a 2m class telescope comparable to an 8m one in a conventional site in terms of sensitivity. The Antarctic Plateau is thus likely to become in the coming decades a new major platform for advanced astronomical observations. In the short/medium term, it is expected that synoptic infrared imaging and spectro-imaging surveys exploring the temporal dimension, from Antarctica could play an essential role in the context of the future large ground based and space projects such as E-ELT, JWST, LSST, EUCLID, GAIA. A medium/large aperture telescope on the Antarctic Plateau have the potential to undertake tasks previously thought to be possible only in space, for example the imaging and crude spectroscopy of Earth-size extra-solar planets [START_REF] Angel | Concept for a second Giant Magellan Telescope (GMT) in Antarctica[END_REF]. Three astrophysical key-domains for an optical/IR facility had been identified, for instance in the ARENA vision (2010), namely (1) the search and characterization of Type Ia supernovae and the investigation of the distant Universe, (2) the stellar populations and evolution in the local group of galaxies, (3) the identification and characterization of extra-solar planets. We present here a preliminary new concept of telescope able to achieve the three above mentioned key-domains. Such concept offers a wide range of cutting edge technological investigations that would definitely overcome the scope of Antarctic astronomy.

DOME C SITE QUALITY

The Dome C unique quality has been reported for day-and night-time astronomical observations and presented by [START_REF] Aristidi | Dome C site testing: surface layer, free atmosphere seeing, and isoplanatic angle statistics[END_REF] and Gredel (2009). At a glance here below are the main results of atmospheric measurements carried out, so far: Precipitable water vapour (PWV) -Monthly average ranges from 0.72(+/-0.20) mm in December to 0.26(+/-0.1) mm in March/April period. Observations in the 200 µm window opens at a transmission level of about 20% during 25% of the time. Observations at 350 µm and 450 µm are possible all year round (see, [START_REF] Tremblin | Site testing for submillimetre astronomy at Dome C, Antarctica[END_REF]. Note that low PWV also results in an improved transmission in the NIR and MIR windows and in an extended wavelength coverage (notably in the important 3µm window).

Sky Brightness -Beyond 2.3µm, the dominant factor is the emission of the molecules (mainly H 2 O) and aerosols at the temperature of the atmosphere (-50 o C to -90 o C). Thanks to these lower temperatures and the lower amount of PWV, the sky is dramatically darker than in the other observatories by at least 3 magnitudes [START_REF] Phillips | The near-infrared sky emission at the South Pole in Winter[END_REF][START_REF] Walden | First Measurements of the infrared sky brightness at Dome C Antarctica[END_REF]. In the visible, sky is dark at solar elevation of -12° resulting in a total of 2506 dark hours per year. The fraction of clear nights is of the order of 80-85% during winter.

Temperature and Wind Profiles -Quasi-periodic ground temperature oscillations occur during summer and winter months. In winter dT/dt is up to 30 o C/week., in summer dT/dt > 10 o C/hour. Wind speed is very low [START_REF] Aristidi | An analysis of temperatures and wind speeds above Dome, Antarctica[END_REF]. Ground layer wind profiles do not show any strong diurnal variations.

Seeing -A summary of the seeing measurements at Dome C has been reported by [START_REF] Aristidi | Dome C site testing: surface layer, free atmosphere seeing, and isoplanatic angle statistics[END_REF] and [START_REF] Ziad | First statistics of the turbulence outer scale at Dome C[END_REF]. Data obtained with three DIMMs set up at 3m, 8m, and 20m height above the ground level show that the seeing is of poor quality on the ground because of a turbulent boundary layer of mean height ~ 23 m in winter. However, the free atmosphere seeing above the boundary layer is indeed very good with a median value of ~ 0.36 arcsec [START_REF] Aristidi | Dome C site testing: surface layer, free atmosphere seeing, and isoplanatic angle statistics[END_REF]. The turbulent layer below the boundary layer can be corrected very efficiently using adaptive optics allowing to achieve 0.40 arcsec over 60 degrees FOV in H and K bands (Le Roux et al., 2007).

AN OFF-AXIS TELESCOPE FOR DOME C

Astrophysical key programs as identified above (1)+( 2 Diffraction limited images of faint objects should be soon obtained with existing 8-m class ground based telescopes, equipped with laser guide stars adaptive optics systems. Ultra deep images, particularly in the 2-5μm region where the natural zodiacal background is very low, should be obtained by the James Webb space telescope, operating at cryogenic temperature. Observations to complement the JWST will be possible from the ground, they will require a very large telescope with the lowest possible thermal background and excellent emissivity and scattered light performance. An off-axis telescope provides an inherently low scattered light design because there are no obstructions in the beam. There are a minimal number of scattered light sources. All mirrors can be robustly supported and articulated because of the easy access allowed by this design. All warm components, sources of the telescope self-thermal emission, will be out of optical beam minimizing its emissivity.

Here we refer to an off-axis telescope configuration as an off-axis, or decentered section of a concentric, on-axis optical configuration in which the primary and auxiliary mirrors are an off-axis piece of a concentric parent mirror. Optically speaking an off-axis telescope is not an asymmetric system, it is a decentered system preserving its bi-lateral symmetry. This means that the system preserves many of the optical performances, tolerances and sensitivities characteristics of the parent concentric system -a deal for system opto-mechanical alignment. In fact, since only a small section of the parent optical system is illuminated, some geometrical performance characteristics of the daughter system supersede the parent optics.

A three-mirror decentered Paul-Baker [START_REF] Paul | Systèmes Connecteurs pour Réflecteurs Astronomique[END_REF][START_REF] Baker | On Improving the Effectiveness of Large Telescopes[END_REF] system is our guideline concept reference, where some variation on the shape of the three mirrors is allowed to achieve good image quality across a wide field of view. One of the main constraints during such optimizations is the necessity to have compact designs, making it possible to be installed on top of a ~23m high (or higher) tower to overcome the boundary layer developed some 2m above ice and up to 23m (Gredel 2009). Above that, the free atmospheric seeing reach its best value during the winter, see A baseline design for Dome C relies on a 2.5 meter unobstructed aperture M1 which should produce a F/8 system optimized over an 1x1 Deg FOV, hereafter dubbed the 2.5@C. We generate a Paul-Baker configuration by adding auxiliary convex secondary mirror (M2) and a concave tertiary (M3). This configuration produces excellent correction across wide FOV's but typically (concentric configuration) also suffers from a lot of obscuration because the large size of the secondary and the tertiary mirrors. Also a concentric configuration produces a trapped focal plane between M2 and M3. This obscuration is avoided if we consider a decentered (off-axis) optical configuration.

The decision on plate scale (PS=10.31Arcsec/mm) was a compromise between an optimal sampling for high resolution and a compact design. An effective focal length (EFL) of 20m produces a support structure long of V1V2 = V1FP < 7m, short enough for our proposes. As shown in Figure 3 we consider a typical detector pixel size of 15µm, which for an EFL=20m results in a sampling of FWHM=2.35 pixels. This is the Nyquist criterion, the standard deviation. Some discussion has been flown confirming that a sampling of FWHM=3.5 pixels or more is the optimal for high resolution. The optical performance at the edges of 1x1 Deg 2 FOV is show in Figure 4. The PSF spots at corneredges (F5, F6, F8, F9) of 1X1 Deg 2 FOV are diffraction limited at 1000nm wavelength and elsewhere they are diffraction limited at 550nm or even better. Note that the blur in this system is only weakly dependent on the off-axis angle and the telescope will be entirely seeing limited.

The optical performance and opto-mechanical characteristic of such concept show its versatility and optimal solution to minimize the light scattering and emissivity and capitalizing the Antarctica/Dome C exceptional seeing and sky transparence on near and thermal-IR.

In terms of feasibility there are several myths about off-axis telescopes. The assessment of new technologies count on the off-axis telescopes already built, for example SOLARC A and NST B . Both act as pathfinders for future off-axis telescope projects to improve the maturity of evolving technologies; enabling more thoroughly understand how to overcome the technical and engineering challenges specific to this kind of design. In particular, resolving issues of alignment, heat dissipation and guiding off of asymmetrical phenomena.

A preliminary opto-mechanical sensitivity for the 6.5m off-axis NPT C was done by [START_REF] Moretto | Optical Performance of the 6.5-m Off-Axis New Planetary Telescope[END_REF]. The NPT opto-mechanical system required accurate positioning to 0.1 mm and 20-40 arcsec. Even a relatively floppy mechanical mount with a flexure of 1 mm in combination with commercial hexapod secondary mount mechanisms and a simple metrology system should be able to maintain telescope collimation. As the 2.5@C the NPT mirrors (M1, M2, and M3) do not physically contain their vertices (so there is no intersection between them and the parent primary mirror optical axis), it is straightforward to define an optical and mechanical fiducial surface near the focal plane that is referenced to the off-axis primary mirror for deterministic or active metrology. In view of the 2.5@C's M1, M2, M3 and FP vertices co-linearity (it was not the case for the NPT), it should make easier its opto-mechanical sensitivities and its alignment and collimation.

The possibility of building a large off-axis optical telescope depends critically on the feasibility of the optics (see Moretto, Kuhn and Goode, 2012, at this conference). Steward Observatory Mirror Lab (SOML) (Arizona, USA) had just finish to polish the largest off-axis segment mirror so far; an 8.4m in diameter off-axis paraboloid mirror. The most critical parts of polishing and tests have been demonstrated in the polishing and measurement of a 1.7m off-axis prototype and the primary mirror for the New Solar Telescope (NST) (see [START_REF] Martin | Production of 8.4m segments for the Giant Magellan Telescope[END_REF], at this conference). Sagem (France), through its Reosc department has also started polishing and testing prototypes segments for the ELT-M1 (see [START_REF] Rodolfo | Prototypes segment polishing and testing for the ELT M1[END_REF]

this conference)

A SOLAR C -Scatter-free Observatory for Limb Active Regions and Coronae is a 0.5m off-axis coronagraphic reflecting telescope located on Haleakala adjacent to the Mees Solar Observatory. Institute for Astronomy, University of Hawaii. B NST -New Solar Telescope is a 1.6m off-axis New Solar Telescope at the BBSO/NJIT. A collaboration with the University of Hawaii and the Korea Astronomy & Space Science Institute. See Goode et al (2012) at this conference. C NPT -New Planetary Telescope, a 6.5m off-axis telescope design study, encompassing many of the mechanical, optical fabrication issues, and facility concerns, has been undertaken by a NPT study group that was organized by the Institute for Astronomy and the Infrared Telescope Facility (IRTF) during 1999. To support the NASA mission, the primary objectives of the NPT are studies of Kuiper Belt objects, near-Earth objects, and direct support of space missions such as Cassini.

CONCLUSIONS

Many areas of modern astrophysics are not flux-limited but are rather dynamic-range limited. Simply collecting more photons will not solve the problem -for these topics, we don't need bigger telescopes, we need better telescopes set up in exceptional sites.

A better telescope concept, based on an off-axis optical design that allows exceptional dynamic range for photometry and high angular resolution + wide field imaging and capitalizes the exceptional atmospheric and environmental Antarctic conditions for astronomical observations over the optical and thermal infrared wavelengths is proposed.

The Antarctic Plateau and, in particular the French-Italian station Concordia, is indeed likely to provide in the future the best site on Earth to set up a medium size infrared telescope facility thank to both, exceptional turbulence at low height and minimal sky brightness conditions. To be costefficient and to take maximum benefit from this unique environment, the instrument and focal instrumentation must be fully optimized. For this reason, we propose an off-axis telescope concept that minimizes stray light due to the diffraction by mirror supports, and instrumental thermal selfemissivity. Since most of the turbulence is concentrated in a thin layer, the telescope must be set up above the boundary layer, i.e., at about 30 m above the ground at Dome C (perhaps at lower elevation in even better sites) and makes use of Adaptive optics device to alleviate the residual turbulence.

)+( 3 )

 3 for an optical/IR facility call for three driver cases (a) the highest possible dynamic range for photometry, (b) angular resolution and (c) the wide-field imaging in the optical and the thermal infrared. A telescope facility complying with the three driver cases (a)+(b)+(c) and capitalizing the Dome C outstanding atmospheric and environmental performance over the optical and thermal infrared wavelengths, assert for an off-axis telescope concept.

Figure 1 -

 1 Figure 1 -(A) The French/Italian Astro Concordia base at Dome C -showing in the backline the staff headquarters (two cylindrical structures) and in the frontline the seeing measurements instruments. The data for winter and summer seeing measurements are presented respectively by (B) and (C), from Aristidi et al. (2009).

Figure 1 (

 1 [START_REF] Aristidi | Dome C site testing: surface layer, free atmosphere seeing, and isoplanatic angle statistics[END_REF].

Figure 2 -

 2 Figure2-The 2.5@C an off-axis telescope for Dome C. A 2.5 meter unobstructed aperture M1 yielding in a F/8 system optimized over a 1x1 Deg FOV. AA' is the parent mirrors M1, M2 and M3 optical axis. M1 and M3 vertexes V1 and V3 are coincident as well as for M2 (V2) and FP. Note that the primary mirror M1, the secondary M2 and the tertiary M3 are not tilted mirrors but are a decentered piece of concentric system, preserving its bi-lateral symmetry.

Figure 3 -Figure 4 -

 34 Figure 3 -The EFL versus sampling for several commercially ready pixel sizes and considering a seeing of 0.36 Arcsec for the Dome C.

Proc. of SPIE Vol. 8444, 84445E-1 Proc. of SPIE Vol. 8444, 84445E-2

Proc. of SPIE Vol. 8444, 84445E-6

Last Picture -A naked eye view of our galaxy from Dome C base -© Jonathan Zaccaria (www.collectifpolaire.org).

Proc. of SPIE Vol. 8444, 84445E-7

Proc. of SPIE Vol. 8444, 84445E-8