Gaël Varoquaux

R Cameron Craddock

Learning and comparing functional connectomes across subjects

Keywords: Functional connectivity, connectome, group study, effective connectivity, fMRI, resting-state

Functional connectomes capture brain interactions via synchronized fluctuations in the functional magnetic resonance imaging signal. If measured during rest, they map the intrinsic functional architecture of the brain. With task-driven experiments they represent integration mechanisms between specialized brain areas. Analyzing their variability across subjects and conditions can reveal markers of brain pathologies and mechanisms underlying cognition. Methods of estimating functional connectomes from the imaging signal have undergone rapid developments and the literature is full of diverse strategies for comparing them. This review aims to clarify links across functional-connectivity methods as well as to expose different steps to perform a group study of functional connectomes.

Introduction

Functional connectivity reveals the synchronization of distant neural systems via correlations in neurophysiological measures of brain activity [START_REF] Biswal | Functional connectivity in the motor cortex of resting human brain using echo-planar MRI[END_REF][START_REF] Friston | Functional connectivity: the principal-component analysis of large (PET) data sets[END_REF]. Given that highlevel function emerges from the interaction of specialized units [START_REF] Tononi | Reentry and the problem of integrating multiple cortical areas: simulation of dynamic integration in the visual system[END_REF], functional connectivity is an essential part of the description of brain function, that complements the localizationist picture emerging from the systematic mapping of regions recruited in tasks [START_REF] Sporns | Organization, development and function of complex brain networks[END_REF]. However, while there exists a well-defined standard analysis framework for activation mapping that enables statistically-controlled comparisons across subjects [START_REF] Friston | Statistical parametric maps in functional imaging: A general linear approach[END_REF], group-level analysis of functional connectivity still face many open methodological challenges. Deriving a picture of a single subject's functional connectivity is by itself not straightforward, as the brain comprises a myriad of interacting subsystems and its connectivity must be decomposed into simplified and synthetic representations. An important view of brain connectivity is that of distributed functional networks depicted by their spatial maps [START_REF] Fox | The human brain is intrinsically organized into dynamic, anticorrelated functional networks[END_REF]. Another no less important and complementary view is that of connections linking localized functional modules depicted as a graph [START_REF] Bullmore | Complex brain networks: graph theoretical analysis of structural and functional systems[END_REF]. This representation of brain connectivity is often called the functional connectome [START_REF] Sporns | The human connectome: a structural description of the human brain[END_REF] and is the focus of intense worldwide research efforts as it holds promises of new insights in cognition and pathologies [START_REF] Biswal | Toward discovery science of human brain function[END_REF][START_REF] Fox | Clinical applications of resting state functional connectivity[END_REF][START_REF] Greicius | Resting-state functional connectivity in neuropsychiatric disorders[END_REF].

The purpose of this paper is to review methodological progress in the estimation of functional connectomes from blood oxygenation level dependent (BOLD) based functional magnetic resonance imaging (fMRI) data and their comparisons across individuals. It does not attempt to be exhaustive, as the field is wide and moving rapidly, but details specific tools and guidelines that, in the experience of the authors, lead to controlled and powerful inter-subject comparisons. The paper is focused on functional connectomes in contrast to structural connectomes, as the inference of functional connectivity requires important statistical modeling considerations that are vastly different from the complications involved with estimating structural connectivity. While the notion of functional connectomics is often associated with the study of resting state [START_REF] Biswal | Toward discovery science of human brain function[END_REF], the methods presented in this paper are also relevant for taskbased studies. On the other hand, the paper has a focus on fMRI; although the core concepts presented can be applied to magnetoencephalography (MEG) or electroencephalography (EEG) [START_REF] Stam | Functional connectivity patterns of human magnetoencephalographic recordings: a "small-world" network?[END_REF], additional specific problems such as source reconstruction must be considered [START_REF] Schoffelen | Source connectivity analysis with MEG and EEG[END_REF].

"Functional connectivity" is defined as a measure of synchronization in brain signals [START_REF] Friston | Functional and effective connectivity in neuroimaging: a synthesis[END_REF]. More generally, it is interesting as a window on underlying synchrony on neural processes [START_REF] Lee | A report of the functional connectivity workshop, dusseldorf[END_REF]. By "functional connectome", here we specifically denote a graph representing functional interactions in the brain, where the term "graph" is taken in its mathematical sense: a set of nodes connected together by edges. Graph nodes (brain regions) correspond to spatially-contiguous and functionally-coherent patches of gray matter and edges describe long-range synchronizations between nodes that are putatively subtended by large fiber pathways [START_REF] Marrelec | Exploring large-scale brain networks in functional MRI[END_REF]. A graph can be weighted or not, and is completely equivalent to its adjacency matrix, a symmetric matrix tabulating the connection weights between each pair of nodes. Functional-connectivity graphs are used to represent evoked activity, as in task-response studies [START_REF] Mcintosh | Towards a network theory of cognition[END_REF], as well as ongoing activity, present in the absence of specific tasks or in the background during task and often studied in so-called resting state experiments [START_REF] Raichle | Two views of brain function[END_REF]. Another important notion that arises from the study of distributed modes of brain function is that of specialized functional networks1 [START_REF] Fox | The human brain is intrinsically organized into dynamic, anticorrelated functional networks[END_REF]. With our definition of the functional connectome, functional networks are not directly building blocks of the connectome but appear as a consequence of the graphical structure [START_REF] Varoquaux | Brain covariance selection: better individual functional connectivity models using population prior[END_REF][START_REF] Varoquaux | Markov models for fMRI correlation structure: is brain functional connectivity small world, or decomposable into networks?[END_REF].

The paper is organized as follows. First we discuss estimation of functional connectomes. This part, akin to a first-level analysis in standard activation mapping methodology, is not in itself a group-level operation, but it is a critical step for inter-subject comparison. In a following section, we discuss several strategies for comparing connectomes across subjects. Finally we discuss the links between the representation of brain connectivity as graphs of functional connectivity and more complex models, such as effective-connectivity models.

Estimating functional connectomes

Here we discuss the inference of connectomes from functional brain imaging data. We start with preprocessing considerations, followed by the choice of nodes i.e. regions, signal extraction, and the estimation of graphs.

Preprocessing considerations

In addition to standard preprocessing performed for task-based analysis (slice-timing correction, realignment, spatial normalization, and possibly smoothing), connectivity-based analysis require additional denoising to separate intrinsic activity from confounding signals. This process involves regressing time series capturing sources of structured noise from the fMRI data. Physiological noise due to cardiac and respiration are two important noise signals [START_REF] Birn | Separating respiratory-variation-related fluctuations from neuronalactivity-related fluctuations in fMRI[END_REF][START_REF] Birn | The respiration response function: the temporal dynamics of fMRI signal fluctuations related to changes in respiration[END_REF][START_REF] Hu | Retrospective estimation and correction of physiological fluctuation in functional MRI[END_REF][START_REF] Lund | Non-white noise in fMRI: does modelling have an impact?[END_REF] that are difficult to control for and as a result are not commonly regressed out. Instead the mean signal from white matter (WM) and cerebrospinal fluid (CSF) are used as surrogates to measure these sources of noise as well as other scanner induced signal fluctuations [START_REF] Fox | The human brain is intrinsically organized into dynamic, anticorrelated functional networks[END_REF][START_REF] Lund | Non-white noise in fMRI: does modelling have an impact?[END_REF]. More complex models account for spatial variation in noise by incorporating voxel-specific regressors of neighboring WM (ANATICOR [START_REF] Jo | Mapping sources of correlation in resting state FMRI, with artifact detection and removal[END_REF]) or the top components from a principal components analysis of highvariance signals (CompCor [START_REF] Behzadi | A component based noise correction method (compcor) for bold and perfusion based fMRI[END_REF]). Head motion induced signal fluctuations are accounted for by incorporating movement parameters [START_REF] Fox | The human brain is intrinsically organized into dynamic, anticorrelated functional networks[END_REF][START_REF] Friston | Movement-related effects in fMRI time-series[END_REF][START_REF] Lund | Non-white noise in fMRI: does modelling have an impact?[END_REF]. The global mean time series has been proposed as an additional noise regressor that appears to improve the spatial specificity of connectivity results [START_REF] Fox | The human brain is intrinsically organized into dynamic, anticorrelated functional networks[END_REF][START_REF] Fox | The global signal and observed anticorrelated resting state brain networks[END_REF]. This practice has become controversial since the global signal regression introduces negative correlations [START_REF] Chang | Effects of model-based physiological noise correction on default mode network anti-correlations and correlations[END_REF][START_REF] Murphy | The impact of global signal regression on resting state correlations: are anti-correlated networks introduced?[END_REF][START_REF] Saad | Trouble at rest: how correlation patterns and group differences become distorted after global signal regression[END_REF]. Removing these sources of nuisance in addition to linear trends results in more contrasted correlation matrices that improve the delineation of functional structures (fig. 2).

Filtering to remove high frequencies is often performed, based on the initial observation that fluctuations implicated in resting-state functional connectivity are predominately slower than 0.1 Hz [START_REF] Biswal | Functional connectivity in the motor cortex of resting human brain using echo-planar MRI[END_REF][START_REF] Cordes | Frequencies contributing to functional connectivity in the cerebral cortex in "restingstate" data[END_REF]. While high-pass and low-pass filtering decrease the impact of some confounds, recent studies have shown that connectivity is present across the full spectrum of observed frequencies [START_REF] Smith | Temporally-independent functional modes of spontaneous brain activity[END_REF][START_REF] Van Oort | Resting state networks are characterized by high frequency BOLD fluctuations[END_REF]. Regressing out a good choice of confound signals is more specific than frequency filtering, and in our experience gives more contrasted correlation matrices2 . In addition, the recent developments of very rapid acquisition protocols prevent aliasing of the physiological noise with the neural signal and give access to more specific noise confounds than traditional low-TR sequences [START_REF] Boyacioglu | Generalized iNverse imaging (GIN): Ultrafast fMRI with physiological noise correction[END_REF].

It is important to keep in mind that the proposed correction strategies are approximate and not definitive techniques. This has become particularly apparent for head motion with reports that micromovements on the scale of ≤ 0.2 mm can induce artefactual group-level findings even when motion is accounted for in preprocessing [START_REF] Power | Spurious but systematic correlations in functional connectivity mri networks arise from subject motion[END_REF][START_REF] Satterthwaite | An improved framework for confound regression and filtering for control of motion artifact in the preprocessing of resting-state functional connectivity data[END_REF][START_REF] Van Dijk | The influence of head motion on intrinsic functional connectivity MRI[END_REF]. Special care must be taken to adequately control for residual impact of head motion in the group model [START_REF] Satterthwaite | An improved framework for confound regression and filtering for control of motion artifact in the preprocessing of resting-state functional connectivity data[END_REF][START_REF] Van Dijk | The influence of head motion on intrinsic functional connectivity MRI[END_REF].

Defining regions

The choice of regions of interests (ROIs) that define the nodes of the graphs can be very important both in the estimation of connectomes and for group comparison [START_REF] Wang | Parcellationdependent small-world brain functional networks: A restingstate fMRI study[END_REF]. Unsurprisingly, simulations have shown that extracting signal from ROIs that did not match functional units would lead to erronous graph estimation [START_REF] Smith | Network modelling methods for fMRI[END_REF]. Different strategies to define suitable ROIs coexist. While dense parcellation approaches cover a large fraction of the brain [START_REF] Achard | A resilient, low-frequency, small-world human brain functional network with highly connected association cortical hubs[END_REF][START_REF] Bellec | Identification of large-scale networks in the brain using fMRI[END_REF][START_REF] Craddock | A whole brain fMRI atlas generated via spatially constrained spectral clustering[END_REF][START_REF] Varoquaux | Brain covariance selection: better individual functional connectivity models using population prior[END_REF][START_REF] Wang | Parcellationdependent small-world brain functional networks: A restingstate fMRI study[END_REF], this coverage can be traded off to focus on some specific regions, in favor of increased functional specificity and thus better differentiation across networks [START_REF] Dosenbach | A core system for the implementation of task sets[END_REF][START_REF] Greicius | Functional connectivity in the resting brain: a network analysis of the default mode hypothesis[END_REF][START_REF] Varoquaux | Detection of brain functional-connectivity difference in post-stroke patients using group-level covariance modeling[END_REF]. In addition, while ROIs are most often defined as a hard selection of voxels, it is also possible to use a soft definition, attributing weights as with probabilistic atlases, or spatial maps of functional networks extracted from techniques such as independent component analysis (ICA) [START_REF] Kiviniemi | Functional segmentation of the brain cortex using high model order group PICA[END_REF][START_REF] Smith | Temporally-independent functional modes of spontaneous brain activity[END_REF].

Regions from atlases. Atlases can be used to define fullbrain parcellations. Popular choices are the Automatic Anatomic Labeling (AAL) atlas [START_REF] Tzourio-Mazoyer | Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain[END_REF], which benefits from an SPM toolbox, or the ubiquitous Talaraich-Tournoux atlas [START_REF] Talairach | Co-planar Stereotaxic Atlas of the Human Brain: 3-dimensional Proportional System[END_REF]. However, these atlases suffer from major shortcomings; namely i) they were defined on a single subject and thus do not reflect inter-subject variability, and ii) they focus on labeling large anatomical structures and do not match functional layout -for instance only two regions describe the medial part of the frontal lobe in the AAL atlas. Multi-subject probabilistic altases such as the Harvard-Oxford atlas distributed with FSL [START_REF] Smith | Advances in functional and structural MR image analysis and implementation as FSL[END_REF] or the sulci-based structural atlas used in [START_REF] Varoquaux | Brain covariance selection: better individual functional connectivity models using population prior[END_REF] mitigate the first problem, and the high number of regions defined using sulci also somewhat circumvent the second problem (see fig. 1).

Defining regions from the literature. Regions can be defined from previous studies, informally or with systematic meta-analysis. This strategy is used to define the main resting-state networks, such as the default mode network, but may also be useful to study connectivity in taskspecific networks [START_REF] Biswal | Functional connectivity in the motor cortex of resting human brain using echo-planar MRI[END_REF][START_REF] Dosenbach | A core system for the implementation of task sets[END_REF][START_REF] Grillon | Hyperfrontality and hypoconnectivity during refreshing in schizophrenia[END_REF][START_REF] Rissman | Measuring functional connectivity during distinct stages of a cognitive task[END_REF]. The common practice is to place balls of a given radius, 5 or 10 mm, centered at the coordinates of interest. Given that functional networks are tightly interleaved in some parts of the cortex, such as the parietal lobe, care must be taken not to define too many regions that would overlap and lead to mixing of the signal. FMRI-based function definition. Defining regions directly from the fMRI signal brings many benefits. First, it can capture subject-specific functional information. Second, it adapts to the signal at hand and its limitations, such as image distortions or vascular and movement artifacts that are isolated in ICA-like approaches. Lastly, incorporating functional information into regional definition will result in more homogenous regions that better represent connectivity present at the voxel level than anatomicallydefined atlases such as AAL or Harvard-Oxford [START_REF] Craddock | A whole brain fMRI atlas generated via spatially constrained spectral clustering[END_REF]. The simplest approach to define task-specific regions is to use activation maps derived from standard GLM-based analysis in a task-driven study (see for instance [START_REF] Poldrack | Handbook of functional MRI data analysis[END_REF]). Regions are extracted by thresholding the maps, or using balls around the activation peaks. For resting-state studies, unsupervised multivariate analysis techniques are necessary. Clustering approaches extract full-brain parcellations [START_REF] Bellec | Multi-level bootstrap analysis of stable clusters in resting-state fMRI[END_REF][START_REF] Craddock | A whole brain fMRI atlas generated via spatially constrained spectral clustering[END_REF][START_REF] Thirion | Dealing with the shortcomings of spatial normalization: Multi-subject parcellation of fMRI datasets[END_REF][START_REF] Yeo | The organization of the human cerebral cortex estimated by intrinsic functional connectivity[END_REF], and have been shown to segment well-known functional structures from rest data. Alternatively, decomposition methods, such as ICA [START_REF] Beckmann | Probabilistic independent component analysis for functional magnetic resonance imaging[END_REF], can unmix linear combinations of multiple effects and separate out partially-overlapping spatial maps that capture functional networks or confounding effects, as for instance with the presence of vascular structure in functional networks. At high model order, ICA maps define a functional parcellation [START_REF] Kiviniemi | Functional segmentation of the brain cortex using high model order group PICA[END_REF]. Extracting regions from these maps requires additional effort as they can display fragmented spatial features and structured background noise, but incorporating sparsity and spatial constraints in the decomposition techniques leads to contrasted maps that outline many different structures [START_REF] Varoquaux | Markov models for fMRI correlation structure: is brain functional connectivity small world, or decomposable into networks?[END_REF] (see fig. 1). [START_REF] Tzourio-Mazoyer | Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain[END_REF], the Harvard-Oxford atlas, the sulci atlas used in [START_REF] Varoquaux | Brain covariance selection: better individual functional connectivity models using population prior[END_REF], regions extracted by Ncuts [START_REF] Craddock | A whole brain fMRI atlas generated via spatially constrained spectral clustering[END_REF], the resting-state networks extracted in [START_REF] Smith | Correspondence of the brain's functional architecture during activation and rest[END_REF] by ICA, and in [START_REF] Varoquaux | Multi-subject dictionary learning to segment an atlas of brain spontaneous activity[END_REF] by sparse dictionary learning.

Optimal number of regions. Defining an optimal number of regions to use for whole-brain connectivity analysis bears careful consideration. On one hand we desire a sufficiently large number of regions to guarantee that they are functionally homogeneous regions and adequately represent the connectivity information present in the data. On the other hand too many regions will render statistical inference challenging, result in an explosion in computational complexity, and interfere with the interpretability of observed connections. For functional parcellation, crossvalidation methods can be employed to estimate an optimal number of regions based on homogeneity, the ability to reproduce connectivity information present at the voxel scale, and the ability to obtain the same parcellations from independent data [START_REF] Blumensath | Resting-state fMRI single subject cortical parcellation based on region growing[END_REF][START_REF] Craddock | A whole brain fMRI atlas generated via spatially constrained spectral clustering[END_REF]. In general these metrics do not result in an obvious peak at a "best" number of regions, but instead offer a range over which the number of regions can be chosen based on the needs of the analysis at hand. Finally, it is important to keep in mind that there is no universally better parcellation and associated number of regions. From a practical standpoint, these choices will depend on the task at hand, and more fundamentally, a good description of brain function should cover multiple scales. Given that it is not clear that an optimal parcellation can be identified from the sample size of a typical study, randomized parcellation, as used in structural connectomes [START_REF] Zalesky | Whole-brain anatomical networks: does the choice of nodes matter?[END_REF] or activation mapping [START_REF] Varoquaux | Small-sample brain mapping: sparse recovery on spatially correlated designs with randomization and clustering[END_REF], may also be considered.

Estimating connections

The concept of functional connectivity has been called elusive [START_REF] Horwitz | The elusive concept of brain connectivity[END_REF]: it has many mathematical instantiations al- though in essence they all strive to extract simple statistics from functional imaging in order to characterize synchrony and communication between large ensembles of neurons.

Here we choose to focus on second order statistics that can be related to Gaussian models, the simplest of which being the correlation matrix of the signals of the different ROIs.

Signal extraction. Given a set of graph nodes, the next step is to extract a representative time series for each node.

To study intrinsic activity, e.g. with rest data, signal extraction can be achieved by either averaging the fMRI time series across the voxels in a region, or by taking the first eigenvariate from a principle components analysis of the time series [START_REF] Friston | A critique of functional localisers[END_REF]. Comparisons of these methods has shown that the eigenvariate method is more sensitive to function inhomogeneity [START_REF] Craddock | A whole brain fMRI atlas generated via spatially constrained spectral clustering[END_REF] and exhibits worse test-retest reliability than averaging time series [START_REF] Zuo | Mean or SVD? A test-retest reliability perspective on seed timeseries generation in RSFC[END_REF]. In addition, improved specificity to BOLD signal can be enforced by using only signal in voxels near gray-matter tissues. For this purpose, we suggest summarizing the signal in an ROI by a mean of the different voxels weighted by the subjectspecific gray matter probabilistic segmentation, as output by e.g. SPM's segmentation tool [START_REF] Ashburner | Unified segmentation[END_REF] or FSL's FAST program [START_REF] Zhang | Segmentation of brain MR images through a hidden Markov random field model and the expectation-maximization algorithm[END_REF].

Studying connectivity from evoked activity with taskdriven studies requires disambiguating task-specific connectivity effects from intrinsic connectivity mediated by shared neuromodulatory/task inputs, anatomical pathways, etc. In this regard, it can be beneficial to run a GLM-based first-level analysis, enforcing specificity of the measure extracted to the task. With slow event-related designs, task-specific functional connectivity can be captured in trial-to-trial fluctuations in the BOLD response, estimated using a GLM analysis with one regressor per trial [START_REF] Grillon | Hyperfrontality and hypoconnectivity during refreshing in schizophrenia[END_REF][START_REF] Mennes | Inter-individual differences in resting-state functional connectivity predict task-induced BOLD activity[END_REF][START_REF] Rissman | Measuring functional connectivity during distinct stages of a cognitive task[END_REF]. This approach, known as beta-series regression, has been adapted for rapid event-related designs, using multiple GLMs to optimize deconvolution of each trial [START_REF] Mumford | Deconvolving bold activation in event-related designs for multivoxel pattern classification analyses[END_REF].

Correlation and partial correlations. Given ROIs defining the nodes of the functional-connectome graph, one needs to estimate the corresponding edges connecting them. Functional connectivity between the ROIs can be measured by computing the correlation matrix of the extracted signals. An important and often neglected point is that the sample correlation matrix, i.e. the correlation matrix obtained by plugging the observed signal in the correlation matrix formula, is not the population correlation matrix, i.e. the correlation matrix of the data-generating process. If the number of measurements was infinite, the two would coincide, however if this number is not large compared to the number of connections (that scales as the square of the number of ROIs), the sample correlation matrix is a poor estimate of the underlying population correlation matrix. In other words, the sample correlation matrix captures a lot of sampling noise, intrinsic randomness that arises in the estimation of correlations from short time series. Conclusions drawn from the sample correlation matrix can easily reflect this estimation error. Varoquaux et al. [START_REF] Varoquaux | Brain covariance selection: better individual functional connectivity models using population prior[END_REF] and Smith et al. [START_REF] Smith | Network modelling methods for fMRI[END_REF] have shown respectively on rest fMRI and on realistic simulations that a good choice of correlation matrix estimator could recover the connectivity structure, where the sample correlation matrix would fail. In general, the choice of a better estimate depends on the settings and the end goals [START_REF] Varoquaux | Detection of brain functional-connectivity difference in post-stroke patients using group-level covariance modeling[END_REF][START_REF] Varoquaux | Markov models for fMRI correlation structure: is brain functional connectivity small world, or decomposable into networks?[END_REF], however the Ledoit-Wolf shrinkage estimate [START_REF] Ledoit | A well-conditioned estimator for largedimensional covariance matrices[END_REF] is a simple, computationally-efficient, and parameter-less alternative that performs uniformly better than the sample correlation matrix [START_REF] Varoquaux | Brain covariance selection: better individual functional connectivity models using population prior[END_REF][START_REF] Varoquaux | Markov models for fMRI correlation structure: is brain functional connectivity small world, or decomposable into networks?[END_REF] and should always be preferred.

For the problem of recovering the functionalconnectivity structure, i.e. finding which region is connected to which, sparse inverse covariance estimators have been found to be efficient [START_REF] Ryali | Estimation of functional connectivity in fMRI data using stability selection-based sparse partial correlation with elastic net penalty[END_REF][START_REF] Smith | Network modelling methods for fMRI[END_REF][START_REF] Varoquaux | Brain covariance selection: better individual functional connectivity models using population prior[END_REF]. The intuition for relying on inverse covariance rather than correlation stems from that fact that standard correlation (marginal correlation) between two variables a and b also capture the effects of other variables: strong correlation of a and b with a third variable c will induce a correlation between a and b. On the opposite, the inverse covariance3 matrix (also called precision matrix ) captures partial correlations, removing the effect of other variables [START_REF] Marrelec | Partial correlation for functional brain interactivity investigation in functional MRI[END_REF]. In the small sample limit, this removal is challenging from the statistical standpoint. This is why an assumption of sparsity, i.e. that only few variables need to be considered at a time, is important to estimate a good inverse covariance. Various estimation strategies exist for sparse inverse covariance, and have an impact on the resulting networks [START_REF] Varoquaux | Brain covariance selection: better individual functional connectivity models using population prior[END_REF][START_REF] Varoquaux | Markov models for fMRI correlation structure: is brain functional connectivity small world, or decomposable into networks?[END_REF]. The Graph Lasso (ℓ 1 -penalized maximum-likelihood estimator) [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF] is in general a good approach for structure recovery. In group studies, the ℓ 21 estimator [START_REF] Honorio | Simultaneous and group-sparse multitask learning of gaussian graphical models[END_REF][START_REF] Varoquaux | Brain covariance selection: better individual functional connectivity models using population prior[END_REF] is useful to impose a common sparsity structure across different subjects and achieve better recovery of this common structure. Simply put, these approaches are necessary because estimation noise creates a background structure (see fig. 3); however, unlike in a univariate situation, the parameters are not independent, and the spurious background connections degrade the estimation of the actual connections. The sparse estimators make a compromise between imposing simpler models, i.e. with less connections, and providing a good fit to the data. This compromise is set via a regularization parameter which controls the sparsity of the estimate. A good procedure to choose this parameter is via cross-validation [START_REF] Varoquaux | Brain covariance selection: better individual functional connectivity models using population prior[END_REF].

Network structure extracted. The correlation matrices and inverse-covariance matrices that we extract contain a lot of information on the functional structure of the brain. First, the correlation matrix (fig. 2) shows blocks of synchronized regions that can be interpreted as large-scale functional networks, such as the default mode network. Note that the split in networks is not straightforward. Different ordering of the nodes will reveal different networks. Indeed, because of the presence of hubs and interleaved networks, the picture in terms of segregated networks is not sufficient to explain full-brain connectivity [START_REF] Varoquaux | Markov models for fMRI correlation structure: is brain functional connectivity small world, or decomposable into networks?[END_REF]. Connectivity matrices, correlation matrices and inverse-covariance matrices, can be represented as graphs: nodes connected by weighted edges (inserts on fig. 2 and fig. 3). The inverse-covariance matrix, which captures partial correlations, appears then as extracting a backbone or core of the graph. While such structure has been used as a way to summarize anatomical brain connectivity graphs [START_REF] Hagmann | Mapping the structural core of human cerebral cortex[END_REF], here it has a clear-cut meaning with regard to the BOLD signal: it gives the conditional independence structure between regions [START_REF] Varoquaux | Markov models for fMRI correlation structure: is brain functional connectivity small world, or decomposable into networks?[END_REF]. In other words, regions a and b are not connected if the signal that they have in common can be explained by a third region c. In this light, the choice of nuisance regressors to remove confounding common signal is less critical with partial correlations than with correlations. Indeed, while with correlation matrices regressing out the global mean has a drastic effect (fig. 2 upper right and lower right), on inverse covariance it only changes the resulting matrices very slightly (fig. 4).

There have been debates on whether to regress out certain signals, such as the global mean, as it induces negative correlations [START_REF] Chang | Effects of model-based physiological noise correction on default mode network anti-correlations and correlations[END_REF][START_REF] Fox | The global signal and observed anticorrelated resting state brain networks[END_REF][START_REF] Murphy | The impact of global signal regression on resting state correlations: are anti-correlated networks introduced?[END_REF], and these may seem surprising: one network appears as having opposite fluctuations to another. However, correlation between two signals only takes its meaning with the definition of a baseline. A simple picture to explain anti-correlations between two regions is the presence of a third region, mediating the interactions. Using this third region as a baseline would amount to estimate partial correlations in the whole system. Using inverse-covariance matrices or partial correlations to understand brain connectivity makes the interpretation in terms of interactions between brain regions easier and more robust to the choice of confounds.

Comparing connectivity

We now turn to the problem of comparing functional connectivity across subjects or across conditions.

Detecting changing connections

First, we focus on detecting where the connectivity matrices estimated in the previous section differ.

Mass-univariate approaches. The most natural approach is to apply a linear model to each coefficient of the connectivity matrices [START_REF] Grillon | Hyperfrontality and hypoconnectivity during refreshing in schizophrenia[END_REF][START_REF] Lewis | Learning sculpts the spontaneous activity of the resting human brain[END_REF]. This approach is similar to the second-level analysis used in mass-univariate brain mapping, and gives rise to many of the well-known techniques used in such a context, such as the definition of a secondlevel design, with possibly the inclusion of confounding effects, and statistical tests (T tests or F tests) on contrast vectors. Importantly, in order to work with Gaussiandistributed variables, it is necessary to apply a Fisher Z transform 4 to the correlations. Note that in these settings, the Ledoit-Wolf estimator [START_REF] Ledoit | A well-conditioned estimator for largedimensional covariance matrices[END_REF] is often a good choice to estimate the correlation matrix, as it is parameter-free and gives good estimation performance without imposing any restrictions on the data. For hypothesis testing, correcting for multiple comparison can severely limit statistical power, as the number of tests performed scales as the square of the number of regions used. Controlling for the false discovery rate (FDR) mitigates this problem. Alternatively, as the assumptions underlying the Benjamini-Hochberg procedure [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF] for the FDR can easily be broken, non-parametric permutations-based tests give reliable approaches. In particular, the max-T procedure [START_REF] Ge | Resampling-based multiple testing for microarray data analysis[END_REF][START_REF] Nichols | Nonparametric permutation tests for functional neuroimaging: a primer with examples[END_REF] is interesting to avoid the drastic Bonferroni correction when controlling for multiple comparison in family-wise error rate.

Accounting for distributed variability. A specific challenge of connectivity analysis is that the connectivity strength between different regions tends to covary. For instance, with resting-state data, functional networks comprising many nodes can appear as more or less connected across subjects (see for instance fig. 5, showing variability in a control population at rest). In other words, non-specific variability is distributed across the connectivity graph, and it is structured by the graph itself. This observation brings the natural question of whether secondlevel analysis should be performed on correlation matrices, inverse-covariance matrices, or another parametrization that would disentangle effects and give unstructured (white) residuals. While inverse-covariance matrices show less distributed fluctuations than correlation matrices, they capture a lot of background noise, as partial correlations are intrinsically harder to estimate. Preliminary work [START_REF] Varoquaux | Detection of brain functional-connectivity difference in post-stroke patients using group-level covariance modeling[END_REF] suggests performing statistical tests on residuals of a parametrization intermediate between correlation matrix and inverse covariance matrix, as it can decouple effects and noise.

Taking a different stance on distributed variability, the "network-based statistics" approach [START_REF] Zalesky | Network-based statistic: Identifying differences in brain networks[END_REF] draws from the hypothesis that if, in a second-level analysis, an effect is detected on a connection that lies in a network of strongly connected nodes, a large sub-network is likely to carry an 4 See http://en.wikipedia.org/wiki/Fisher_transformation or [START_REF] Anderson | An introduction to multivariate statistical analysis[END_REF] e. Z score on residuals [114] Figure 5: Inter subject variability. Note that this is variability occurring in a healthy population at rest, in other words it is non specific variabilitya: single-subject correlation matrices for different subjectsb: Corresponding Z-score (effect / standard deviation) of the difference between a subject and the remaining othersc: singlesubject inverse-covariance matricesd: Corresponding Z-score for the inverse-covariance matricese: Corresponding Z-score for the subject residuals, as defined in [START_REF] Varoquaux | Detection of brain functional-connectivity difference in post-stroke patients using group-level covariance modeling[END_REF].

effect. Thus, they adapt cluster-level inference to connectivity analysis, in order to mitigate the curse of multiple comparisons.

Comparing network summary statistics

Both the multiple comparison issue and the networklevel distributed variability are a plague to edge-level comparison of connectomes. A possible strategy to circumvent these difficulties is to perform comparisons and statistical testing at the level of the network, rather than the individual connection.

Network integration. Marrelec et al. [START_REF] Marrelec | Regions, systems, and the brain: hierarchical measures of functional integration in fMRI[END_REF] introduce the use of entropy and mutual information as a measure of network-level functional integration5 . Gaussian entropy can be seen as a simple metric to generalize correlation or variance to multiple nodes (see [START_REF] Anderson | An introduction to multivariate statistical analysis[END_REF] §2.5.2 and §7.5). Indeed, let us consider 3 nodes: a, b and c. Their correlation structure is captured by three correlation coefficients: ρ ab , ρ bc and ρ ac . Summarizing these by their mean, as might seem natural, discards the relationship between the signals, while using the integration metric, defined as the Gaussian entropy, tells us how much two signals can be combined to form the third (see fig. 6). Cross-entropy -or mutual information- [START_REF] Marrelec | Regions, systems, and the brain: hierarchical measures of functional integration in fMRI[END_REF] measures the amount of crosstalk between two systems in a similar way as Gaussian entropy is used to measure the integration of a brain system. The functional-connectivity structure, or its representation in the form of a correlation matrix, can thus be characterized via the integration and cross-talk of some of its sub-systems. This approach gives a simplified representation with a small number of metrics that can be compared across subjects.

2.422

Graph-topological metrics. Functional connectivity graphs have been found to display specific topological properties6 that are characteristic of small-world networks [START_REF] Achard | A resilient, low-frequency, small-world human brain functional network with highly connected association cortical hubs[END_REF][START_REF] Bullmore | Complex brain networks: graph theoretical analysis of structural and functional systems[END_REF][START_REF] Salvador | Neurophysiological architecture of functional magnetic resonance images of human brain[END_REF][START_REF] Stam | Functional connectivity patterns of human magnetoencephalographic recordings: a "small-world" network?[END_REF]. These networks display excellent transport properties: although they have a relatively small number of connections, any two regions of the brain are well connected. Another interesting consequence of their specific topology is the resilience it gives the system to attacks such as resulting from brain lesions [START_REF] Achard | A resilient, low-frequency, small-world human brain functional network with highly connected association cortical hubs[END_REF]. This overall structure of functional-connectivity graphs can be summarized by a few metrics, such as the average path length between any two nodes, the local clustering coefficients, or the node degree centrality [START_REF] Rubinov | Complex network measures of brain connectivity: uses and interpretations[END_REF]. Given that pathologies without a localized focus, such as schizophrenia, are thought to have a global impact on brain connectivity [START_REF] Bassett | Hierarchical organization of human cortical networks in health and schizophrenia[END_REF][START_REF] Liu | Disrupted small-world networks in schizophrenia[END_REF], the graph-topological metrics are promising markers to perform inter-subject comparison. Such an approach is appealing as it is not subject to multiple comparison issues. However, it has been criticized as giving a fairly unspecific characterization of the brain and being fragile to noise [START_REF] Ioannides | Dynamic functional connectivity[END_REF]. Another caveat is that these properties are not specific to brain function: correlation matrices display smallworld properties such as local clustering by construction. Indeed, if two nodes are strongly correlated to a third, they are highly likely to be correlated to each other [START_REF] Zalesky | On the use of correlation as a measure of network connectivity[END_REF]. This observation highlights the need for well defined nullhypothesis [START_REF] Rubinov | Weight-conserving characterization of complex functional brain networks[END_REF][START_REF] Zalesky | On the use of correlation as a measure of network connectivity[END_REF], but also for controlled recovery of brain functional connectivity going beyond empirical correlation matrices, as discussed in the previous section.

Predictive Modeling

Predictive modeling is concerned with learning (or fitting) a model that is capable of predicting information from unseen data [START_REF] Pereira | Machine learning classifiers and fmri: a tutorial overview[END_REF]. In the context of connectomes, predictive modeling can extract connectivity-based biomarkers of disease diagnosis, prognosis, or other phenotypic outcomes [START_REF] Craddock | Disease state prediction from resting state functional connectivity[END_REF][START_REF] Dosenbach | Prediction of individual brain maturity using fmri[END_REF]. The accuracy of a predictive model provides a measure of the amount of information present in the connectome about the phenotypic measure being evaluated [START_REF] Kjems | The quantitative evaluation of functional neuroimaging experiments: Mutual information learning curves[END_REF][START_REF] Kriegeskorte | Information-based functional brain mapping[END_REF]. When combined with reproducibility, prediction accuracy provides a metric for evaluating experimental trade-offs for data acquisition, preprocessing, and analysis [START_REF] Laconte | The evaluation of preprocessing choices in single-subject bold fMRI using npairs performance metrics[END_REF][START_REF] Strother | The quantitative evaluation of functional neuroimaging experiments: the NPAIRS data analysis framework[END_REF]. Multivariate predictive models are attractive in connectomics because they are sensitive to dependencies between features and avoid the need to correct for multiple comparisons since the significance of an entire pattern is evaluated using a single statistical test. Additionally, modern predictive modeling techniques draw from the statistical learning literature, which specifically addresses high dimensional datasets with few observations. Predictive modeling has been successfully applied to identify connectome-based biomarkers of Alzheimer's disease [START_REF] Stonnington | Predicting clinical scores from magnetic resonance scans in alzheimer's disease[END_REF], depression [START_REF] Craddock | Disease state prediction from resting state functional connectivity[END_REF][START_REF] Zeng | Identifying major depression using whole-brain functional connectivity: a multivariate pattern analysis[END_REF], schizophrenia [START_REF] Cecchi | Discriminative network models of schizophrenia[END_REF][START_REF] Shen | Discriminative analysis of resting-state functional connectivity patterns of schizophrenia using low dimensional embedding of fMRI[END_REF], autism [START_REF] Anderson | Functional connectivity magnetic resonance imaging classification of autism[END_REF], ADHD [START_REF] Zhu | Fisher discriminative analysis of resting-state brain function for attention-deficit/hyperactivity disorder[END_REF], aging [START_REF] Dosenbach | Prediction of individual brain maturity using fmri[END_REF], as well as to classify mental operations [START_REF] Richiardi | Decoding brain states from fMRI connectivity graphs[END_REF][START_REF] Shirer | Decoding subject-driven cognitive states with whole-brain connectivity patterns[END_REF]. The growing interest for applying predictive modeling to connectivity analysis was highlighted by the ADHD200 Global Competition, in which the object was to identify a connectivity-based biomarker of ADHD [108]. Recent work has illustrated the utility of predictive modeling for deriving connectivity models at the individual level [START_REF] Chu | Measuring the consistency of global functional connectivity using kernel regression methods[END_REF].

Technically, predictive modeling is a supervised machine learning problem where a target to be predicted -e.g. age, disease state, cognitive state-is available for each observation of the data. In the context of comparing connectomes, features used in the predictive model correspond to bivariate measures of connectivity [START_REF] Dosenbach | Prediction of individual brain maturity using fmri[END_REF][START_REF] Richiardi | Decoding brain states from fMRI connectivity graphs[END_REF][START_REF] Shirer | Decoding subject-driven cognitive states with whole-brain connectivity patterns[END_REF], or any of the previously discussed graph summary metrics [START_REF] Cecchi | Discriminative network models of schizophrenia[END_REF][START_REF] Ekman | Predicting errors from reconfiguration patterns in human brain networks[END_REF]. The quality of a predictive model is determined by its prediction accuracy (or generalization ability) which is measured using one or more iterations of cross-validation. Cross-validation iteratively subdivides available data into a subset used for training the classifier and a dataset for evaluating classifier performance 7 [START_REF] Pereira | Machine learning classifiers and fmri: a tutorial overview[END_REF]. The significance of achieved prediction accuracy can be assessed using permutation tests [START_REF] Golland | Permutation tests for classification: Towards statistical significance in image-based studies[END_REF]. Predictive modeling approaches typically require the specification of several parameters, which may be chosen based on domain specific knowledge or requirements [START_REF] Cherkassky | Learning from Data: Concepts, Theory, and Methods[END_REF], determined using an analytical approach [START_REF] Cherkassky | Practical selection of SVM parameters and noise estimation for svm regression[END_REF], or optimized using a second-level cross validation procedure [START_REF] Friedman | The elements of statistical learning[END_REF]. 7 Several strategies exist for performing cross-validation and the commonly used approach of using only a single observation for testing (leave-one-out cross-validation) results in highly variable estimates of prediction accuracy [START_REF] Friedman | The elements of statistical learning[END_REF]. Alternative approaches such as (5 or 10)fold cross-validation, or 0.632+ bootstrap should be preferred [START_REF] Friedman | The elements of statistical learning[END_REF].

Although predictive modeling techniques are well suited for measuring whether information exists in the connectome about a phenotypic variable, they do not directly identify the connections that are most relevant to the prediction. This limitation can be somewhat mitigated by relying on previously described sparse inverse covariance estimation techniques to minimize the number of connections. Additionally, feature selection [START_REF] Guyon | An introduction to variable and feature selection[END_REF] can be performed by filtering out features based on their statistical relationship with the variable of interest [START_REF] Craddock | Disease state prediction from resting state functional connectivity[END_REF][START_REF] Shen | Discriminative analysis of resting-state functional connectivity patterns of schizophrenia using low dimensional embedding of fMRI[END_REF]. Importantly, it must be performed within cross-validation to avoid biasing estimates of prediction accuracy. The interpretation of connections used in predictive models and their relationship with a phenotypic outcome is difficult and requires insight into the mechanisms underlying the modeling approach. For linear models, the weights of the model are similar to the weights of an ordinary linear regression. If features are appropriately standardized prior to training, the magnitude of the weights can be interpreted to reflect the relative importance of the feature to the model, for instance the corresponding edge in the connectome. However, interpreting how the connections differ between classes or relate to a phenotypic variable can be more complicated given the multivariate nature of their involvement. Indeed, the inclusion or exclusion of a connection in the model can induce a change of the sign of another model weight [START_REF] Craddock | Disease state prediction from resting state functional connectivity[END_REF]. It is perhaps most reasonable to adapt a conservative interpretation in which predictive modeling is used to identify candidate connections that are later tested in follow-up experiments better suited to elucidating their relationship to the variable of interest.

To conclude on predictive modeling with a practical note for connectome comparison, we would like to stress that while machine learning algorithms are powerful tools, they work best if they are provided with discriminanting and noiseless features. In other words, as with all other connectome comparison methods, optimizing firstlevel analysis -subject-level connectome extraction-is paramount.

Beyond correlation, effective connectivity?

All the approaches that we have presented in this review are based on second-order statistics of the signal, in other words correlation analysis. Traditionally, these are defined as functional connectivity, defined as "temporal correlations between remote neurophysiological events" [START_REF] Friston | Functional and effective connectivity in neuroimaging: a synthesis[END_REF], and opposed to effective connectivity, i.e. "the influence one neural system exerts over another" [START_REF] Friston | Functional and effective connectivity in neuroimaging: a synthesis[END_REF]. To conclude this review, we would like to bridge the gap between these concepts, which in our eyes should be seen as a continuum rather than an opposition (this opinion is also expressed in [START_REF] Mcintosh | Moving between functional and effective connectivity[END_REF]).

A first step to move from purely descriptive statistics to interaction models with functional connectivity analysis is to consider a correlation matrix as a Gaussian graphical model, i.e. a well-defined probabilistic model that de-scribes observed correlations in terms of an independence structure and conditional relations [START_REF] Lauritzen | Graphical models[END_REF][START_REF] Varoquaux | Markov models for fMRI correlation structure: is brain functional connectivity small world, or decomposable into networks?[END_REF]. In such settings, the inverse covariance graph or the partial correlations are a measure of influence from one node to another, albeit undirected. Inferring directionality in a Gaussian model is impossible. Linear structural equation models (SEMs) [START_REF] Mcintosh | Structural equation modeling and its application to network analysis in functional brain imaging[END_REF] rely on a similar model that consists in specifying a candidate directed graphical structure. This structure constraints the covariance matrix of the signals and can thus be tested on observed data. In fact some forms of SEMs are known as "covariance structure models". There is thus a strong formal link between correlation analysis in the framework of graphical models and SEMs: the former is undirected but fully exploratory, as it does not require the specification of candidate structure, while the latter is directed but confirmatory. This link has been exploited to specify candidate structures for SEMs using partial correlations [START_REF] Marrelec | Using partial correlation to enhance structural equation modeling of functional MRI data[END_REF]. More complex models, such as dynamical causal models (DCMs) [START_REF] Friston | Dynamic causal modelling[END_REF] or Granger causality [START_REF] Goebel | Investigating directed cortical interactions in time-resolved fmri data using vector autoregressive modeling and granger causality mapping[END_REF] require additional hypotheses such as non-linear couplings or time lags.

Most importantly, more complex models can only be used to model interactions between a small number of nodes. This is not only due to a computational difficulty, but also to fundamental roadblocks in statistics: the complexity of the model must match the richness of the data. While injecting prior information can help model estimation, the more informative this prior is, the more fragile the inference becomes. The ongoing debate on the impact of hemodynamic lag on Granger-causality inference [START_REF] Smith | The danger of systematic bias in group-level FMRI-lag-based causality estimation[END_REF] is an example of such fragility. Note that although most of the theory underpinning correlation analysis (Gaussian graphical models) is based on a Gaussian assumption, the core results are robust to violations of this assumption [START_REF] Ravikumar | Highdimensional covariance estimation by minimizing 1-penalized log-determinant divergence[END_REF].

It is tempting to favor more neurobiologically-inspired models that give descriptions close to our knowledge of the brain basic mechanisms, however, as George Box famously said, "all models are wrong; some models are useful". Depending on the question and the data at hand, a trade-off should be chosen between complex models based on a biophysical description, and simple phenomenological models such as correlation matrices. In particular to model interactions between a large number of regions, as in full-brain analysis, and learn a large connectome, simple models are to be preferred. For more hypothesis-driven studies, such as the analysis of the mechanisms underlying a specific task, more complex models can be preferred, if rich data is available. Automatic choice of model is a difficult problem, however, cross-validation (as used in [START_REF] Craddock | A whole brain fMRI atlas generated via spatially constrained spectral clustering[END_REF][START_REF] Strother | Evaluating fMRI preprocessing pipelines[END_REF][START_REF] Varoquaux | Brain covariance selection: better individual functional connectivity models using population prior[END_REF]) is a useful tool. The central principle of cross-validation is to test a model on different data than the data used to fit the model. Models too complex for the data available will fit noise in the data, and thus generalize poorly. The main benefit of cross-validation is that it is a non-parametric method which does not rely strongly on modeling assump-tions 8 .

Conclusion

Horwitz el al. [START_REF] Horwitz | Network analysis of brain cognitive function using metabolic and blood flow data[END_REF] claimed almost 20 years ago that "the crucial concept needed for network analysis is covariance". In our eyes, this still holds today. Estimation functional connectomes relies largely on fitting covariance models. Their comparison requires understanding how these covariances vary and finding metrics to capture this variability. The additional secret ingredient may be using confounds regressors in all statistical steps. A good choice of a small number of relevant regions facilitates connectome comparison. However, such a choice cannot yet be fully factored out via methods and must rely on neuroscientific expertise.

Methodological challenges to functional-connectomebased group studies arise from the dimensionality and the variability of the connectome. With the current tools, inter-subject comparison of connectomes comprising many nodes is limited by the difficulty of estimating high-dimensional covariance matrices and the loss of statistical power due to multiple comparisons. Better algorithms integrating powerful a priori information are required to push the limits of covariance estimation. Better characterization of inter-subject variability of connectomes [START_REF] Kelly | Characterizing variation in the functional connectome: promise and pitfalls[END_REF] will help choosing parameterizations and invariants to avoid testing each edge for a difference, as this strategy inevitably leads to a needle in a haystack problem.

Reviewing methodological options to learn and compare connectomes highlights that there is currently no unique solution, but a spectrum of related methods and analytical strategies. More empirical results are required to guide the choices. However this diversity is probably unavoidable: a diffuse disease like schizophrenia will not lead to the same connectome modifications as a focal lesion. In statistical learning, "no free lunch" theorems [START_REF] Wolpert | The lack of a priori distinctions between learning algorithms[END_REF] tell us that no strategy can perform uniformly better in all situations. In practice, the key to a successful analysis is to understand well the assumptions and interpretation of each option, in order to match the method to the question. Similarly, the idealized notion of an unique functional connectome to describe connections in brain function is probably an utopia, and various connectomes should be considered in different settings, such as the study of varying 8 This is to be contrasted to Bayesian model comparison, which will give well-controlled results only if the true generative model is in the list of models compared. [START_REF] Friston | Ten ironic rules for non-statistical reviewers[END_REF] argues that, based on the Neyman-Pearson lemma, cross-validation is less powerful than likelihood ratio tests using the full dataset. However, it is important to keep in mind that these approaches only test for self-consistence, as the Neyman-Pearson lemma is established under the hypothesis that the model used to define the test is indeed the data-generating process [START_REF] Neyman | On the problem of the most efficient tests of statistical hypotheses[END_REF], while in practice it is often the case that this model gives poor fits to the data [START_REF] Lohmann | Critical comments on dynamic causal modelling[END_REF]. Applying test procedures on different data than that used to fit the model, as in cross-validation, is much more resilient to modeling errors. phenotypic conditions, or that of on-going activity versus activity related to specific tasks.

Figure 1 :

 1 Figure1: Different full-brain parcellations: the AAL atlas[START_REF] Tzourio-Mazoyer | Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain[END_REF], the Harvard-Oxford atlas, the sulci atlas used in[START_REF] Varoquaux | Brain covariance selection: better individual functional connectivity models using population prior[END_REF], regions extracted by Ncuts[START_REF] Craddock | A whole brain fMRI atlas generated via spatially constrained spectral clustering[END_REF], the resting-state networks extracted in[START_REF] Smith | Correspondence of the brain's functional architecture during activation and rest[END_REF] by ICA, and in[START_REF] Varoquaux | Multi-subject dictionary learning to segment an atlas of brain spontaneous activity[END_REF] by sparse dictionary learning.

Figure 2 :

 2 Figure 2: Correlation matrices of rest time-series extracted from the 39 main regions of the Varoquaux 2011 [115] parcellation with different choices of confound regressors -Left: regressing out CompCor signals, as well as white matter and CSF average signals and movement parameters. The insert shows the connections restricted to a few major nodes. -Upper right: regressing out only movement parameters. -Lower right: regressing out movement parameters and global signal mean. No frequency filtering was applied here. When no confounding brain signals are regressed, all regions are heavily correlated. Regressing out common signal, in the form of well-identified confounds or a global mean, teases out the structure.

R

  Po st Te m p R DL PF C R Pa r R Fr on t po l L Pa r L DL PF C L Fr on t po l L IP S R IP S L An t IP S R An t IP S M ot or L Au d R Au d L ST S R ST S L In s R In s Ci ng V AC C D AC C R A In s Ba sa l Br oc a R Pa rs O p Su p Fr on t S L TP J R TP J Ce re b L LO C R LO C Vi s St ria te O

Figure 3 :

 3 Figure 3:Different inverse-covariance matrices estimates corresponding to fig.2-Left: group-sparse estimate using the ℓ 21 estimator[START_REF] Varoquaux | Brain covariance selection: better individual functional connectivity models using population prior[END_REF]. The insert shows the connections restricted to a few major nodes. -Upper right: non-sparse estimate: inverse of the sample correlation matrix. -Lower right: sparse estimate using the Graph Lasso[START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF].

Figure 4 :

 4 Figure 4: Inverse-covariance matrices for different choice of confound regressors -Left: regressing out only movement parameters -Right: removal of the global mean, instead of the white matter, CSF, and CompCor time courses.

  section 4.2.3 for mathematical arguments.

Figure 6 :

 6 Figure 6: Two different correlation matrices with the same average correlation, but with very different integration values. Indeed, the matrix on the left was chosen to represent three signals a, b and c as different from each other as possible, given ρ ab + ρ bc + ρac = 1.35; it thus has a small integration value. On the opposite, for the matrix on the right, signal b can almost be fully recovered by combining signals a and c; the matrix thus has a large integration value.

In neuroimaging, the term network is sometimes used to denote a graph of brain function. To disambiguate the notion of segregated spatial mode[START_REF] Fox | The human brain is intrinsically organized into dynamic, anticorrelated functional networks[END_REF] from that of connectivity graphs, we will purposely restrict its usage in this paper.

Note that naive use of filtering can induce spurious correlations [26].

Covariance and correlation matrices differ simply by the fact that a covariance matrix captures the amplitude of a signal, via its variance, while a correlation matrix is computed on standardized (zero mean, unit variance) signals.

See[START_REF] Varoquaux | Brain covariance selection: better individual functional connectivity models using population prior[END_REF] for simplified formulas for network integration and mutual information.

In the neuroscience world, these descriptions are grouped under the terms of "graph-theoretical approaches", however graph theory is an entire division of mathematics and computer science that is concerned with much more than topology of random graphs.

Acknowledgments

GV acknowledges funding from the NiConnect grant and the Dynamic Diaschisis project DEQ20100318254 from Fondation pour la Recherche Médicale, as well as many insightful discussions with Andreas Kleinschmidt on on-going activity and Bertrand Thirion on statistical data processing. RCC would like to acknowledge support by a NARSAD Young Investigator Grant from the Brain & Behavior Research Foundation. The authors would like to thank the anonymous reviewers for their suggestions, which improved the manuscript.


Y Maday
email: maday@ann.jussieu.fr

O Mula

A generalized empirical interpolation method : application of reduced basis techniques to data assimilation

Introduction

The representation of some physical or mechanical quantities, representing a scalar or vectorial function that depends on space, time or both, can be elaborated through at least two -possiblycomplementary approaches: the first one, called explicit hereafter, is based on the measurement of some instances of the quantity of interest that consists in getting its value at some points from which, by interpolation or extrapolation, the quantity is approximated in other points than where the measurements have been performed. The second approach, called implicit hereafter, is more elaborated. It is based on a model, constructed by expertise, that implicitly characterizes the quantity as a solution to some problem fed with input data. The model can e.g. be a parameter dependent partial differential equation, the simulation of which allows to get an approximation of the quantity of interest, and, actually, many more outputs than the sole value of the quantity of interest. This second approach, when available, is more attractive since it allows to have a better understanding of the working behavior of the phenomenon that is under consideration. In turn, it facilitates optimization, control or decision making.

Nevertheless for still a large number of problems, the numerical simulation of this model is indeed possible -though far too expensive to be performed in a reasonable enough time. The combined efforts of numerical analysts, specialists of algorithms and computer scientists, together with the increase of the performances of the computers allow to increase every days the domains of application where numerical simulation can be used, to such an extent that it is possible now to rigorously adapt the approximation, degrade the models, degrade the simulation, or both in an intelligent way without sacrificing the quality of the approximation where it is required.

Among the various ways to reduce the problem's complexity stand approaches that use the smallness of the Kolmogorov n-width [START_REF] Kolmogoroff | Über die beste Annäherung von Funktionen einer gegebenen Funktionenklasse[END_REF] of the manifold of all solutions considered when the parameters varies continuously in some range. This idea, combined with the Galerkin method is at the basis of the reduced basis method and the Proper Orthogonal Decomposition (POD) methods to solve parameter dependent partial differential equations. These approximation methods allow to build the solution to the model associated to some parameter as a linear combination of some precomputed solutions associated to some well chosen parameters. The precomputations can be lengthy but are performed off-line, the online computation has a very small complexity, based on the smallness of the Kolmogorov n-width. We refer to [START_REF] Patera | Reduced Basis Approximation and A Posteriori Error Estimation for Parametrized Partial Differential Equations, Version 1.0[END_REF] [START_REF] Prud'homme | Reliable realtime solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF] for an introduction to these approaches.

Another possibility, rooted on the same idea, is the empirical interpolation method (EIM) that allows, from values of the quantity at some interpolating points, to build a linear combination of again preliminary fully determined quantities associated to few well chosen instances of the parameter. The linear combination is determined in such a way that it takes the same values at the interpolating points as the quantity we want to represent. This concept generalizes the classical -e.g. polynomial or radial basis -interpolation procedure and is recalled in the next section. The main difference is that the interpolating function are not a priori known but depend on the quantity we want to represent.

In this paper we first aim at generalizing further this EIM concept by replacing the pointwise evaluations of the quantity by more general measures, mathematically defined as linear forms defined on a superspace of the manifold of appropriate functions. We consider that this generalization, named Generalized Empirical Interpolation Method (GEIM), represents already an improvement with respect to classical interpolation reconstructions.

Bouncing on this GEIM, we propose a coupled approach based on the domain decomposition of the computational domain into two parts : one small domain Ω 1 where the Kolmogorov nwidth of the manifold is not small and where the parametrized PDE will be simulated and the other subdomain Ω 2 , much larger but with a small Kolmogorov n-width because for instance the solution is driven over Ω 2 by the behavior of the solution over Ω 1 . The idea is then to first construct (an approximation of) the solution from the measurements using the GEIM. In turn this reconstruction, up to the interface between Ω 1 and Ω 2 , provides the necessary boundary conditions for solving the model over Ω 1 . This is not the first attempt to use the small Kolmogorov width for another aim than the POD or reduced basis technique which are both based on a Galerkin approach. In [START_REF] Chakir | A two-grid finite-element/reduced basis scheme for the approximation of the solution of parametric dependent PDE[END_REF] e.g. the smallness of the Kolmogorov width is used to post-process a coarse finite element approximation and get an improved accuracy.

The problems we want to address with this coupled approach, stem from, e.g., actual industrial process or operations that work on a day-to-day basis; they can be observed with experimental sensors that provide sound data and are able to characterize part of their working behavior. We think that the numerical simulation and data mining approaches for analyzing real life systems are not enough merged in order to (i) complement their strength and (ii) cope for their weaknesses. This paper is a contribution in this direction.

In the last section, we evoke the problem of uncertainty and noises in the acquisition of the data, since indeed, the data are most often polluted by noises. Due to this, statistical data acquisition methods are used to filter out the source signals so that an improved knowledge is accessible. In many cases though, and this is more and more the case now, the data are far too numerous to all be taken into account, most of them are thus neglected because people do not know how to analyze them, in particular when the measures that are recorded are not directly related to some directly understandable quantity.

Generalized Empirical Interpolation Method

The rationale of all our approach relies on the possibility to approximately represent a given set, portion of a regular manifold (here the set of solution to some PDE), as a linear combination of very few computable elements. This is linked to the notion of n-width following Kolmogorov [START_REF] Kolmogoroff | Über die beste Annäherung von Funktionen einer gegebenen Funktionenklasse[END_REF]:

Definition 1 Let F be a subset of some Banach space X and Y n be a generic n-dimensional subspace of X . The angle between F and Y n is

E(F ; Y n ) := sup x∈F inf y∈Yn x -y X . The Kolmogorov n-width of F in X is given by d n (F, X ) := inf{E(F ; Y n ) : Y n a n-dimensional subspace of X} = inf Yn sup x∈F inf y∈Yn x -y X . (1)
The n-width of F thus measures to what extent the set F can be approximated by an ndimensional subspace of X . We assume from now on that F and X are composed of functions defined over a domain Ω ⊂ I R d , where d = 1, 2, 3 and that F is a compact set of X .

Recall of the Empirical Interpolation Method

We begin by describing the construction of the empirical interpolation method ( [START_REF] Barrault | An "empirical interpolation" method: Application to efficient reduced-basis discretization of partial differential equations[END_REF], [START_REF] Grepl | Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations[END_REF], [START_REF] Maday | A general multipurpose interpolation procedure: the magic points[END_REF]) that allows us to define simultaneously the set of generating functions recursively chosen in F together with the associated interpolation points. It is based on a greedy selection procedure as outlined in [START_REF] Nguyen | Certified real-time solution of parametrized partial differential equations[END_REF][START_REF] Prud'homme | Reliable realtime solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF][START_REF] Veroy | A posteriori error bounds for reduced-basis approximation of parametrized noncoercive and nonlinear elliptic partial differential equations[END_REF]. With M being some given large number, we assume that the dimension of the vectorial space spanned by F : span (F ) is of dimension ≥ M.

The first generating function is

ϕ 1 = arg max ϕ∈F ϕ( • ) L ∞ (Ω)
, the associated interpolation point satisfies x 1 = arg max x∈Ω |ϕ 1 (x)|, we then set q 1 = ϕ 1 (•)/ϕ 1 (x 1 ) and B 1 11 = 1. We now construct, by induction, the nested sets of interpolation points Ξ M = {x 1 , . . . , x M }, 1 ≤ M ≤ M max , and the nested sets of basis functions {q 1 , . . . , q M }, where M max ≤ M is some given upper bound fixed a priori. For M = 2, . . . , M max , we first solve the interpolation problem : Find

I M-1 [ϕ(•)] = M-1 j=1 α M-1,j [ϕ]q j , (2)
such that

I M-1 [ϕ(•)](x i ) = ϕ(x i ), i = 1, . . . , M -1 , (3)
that allows to define the α M-1,j [ϕ], 1 ≤ j ≤ M, as it can be proven indeed that the (M -1) × (M -1) matrix of running entry q j (x i ) is invertible, actually it is lower triangular with unity diagonal.

We then set

∀ϕ ∈ F, ε M-1 (ϕ) = ϕ -I M-1 [ϕ] L ∞ (Ω) , (4)
and define

ϕ M = arg max ϕ∈F ε M-1 (ϕ) , (5)
and

x M = arg max x∈Ω |ϕ M (x) -J M-1 [ϕ M ](x)| , (6)
we finally set r

M (x) = ϕ M (x) -J M-1 [ϕ M (x)], q M = r M /r M (x M ) and B M ij = q j (x i ), 1 ≤ i, j ≤ M .
The Lagrangian functions -that can be used to build the interpolation operator

I M in X M = span {ϕ i , 1 ≤ i ≤ M } = span {q i , 1 ≤ i ≤ M } over the set of points Ξ M = {x i , 1 ≤ i ≤ M } - verify for any given M , I M [u( • )] = M i=1 u(x i )h M i ( • ), where h M i ( • ) = M j=1 q j ( • )[B M ] -1 ji (note indeed that h M i (x j ) = δ ij ).
The error analysis of the interpolation procedure classically involves the Lebesgue constant

Λ M = sup x∈Ω M i=1 |h M i (x)|.
Lemma 2 For any ϕ ∈ F , the interpolation error satisfies

ϕ -I M [ϕ] L ∞ (Ω) ≤ (1 + Λ M ) inf ψM ∈XM ϕ -ψ M L ∞ (Ω) . (7)
The last term in the right hand side of the above inequality is known as the best fit of ϕ by elements in X M .

The generalization

Let us assume now that we do not have access to the values of ϕ ∈ F at points in Ω easily, but, on the contrary, that we have a dictionary of linear forms σ ∈ Σ -assumed to be continuous in some sense, e.g. in L 2 (Ω) with norm 1 -the application of which over each ϕ ∈ F is easy. Our extension consists in defining φ1 , φ2 ,. . . , φM and a family of associated linear forms σ 1 , σ 2 ,. . . , σ M such that the following generalized interpolation process (our GEIM) is well defined :

J M [ϕ] = M j=1 β j φj , such that ∀i = 1, . . . , M, σ i (J M [ϕ]) = σ i (ϕ) (8)
Note that the GEIM reduces to the EIM when the dictionary is composed of dirac masses, defined in the dual space of C 0 (Ω).

As explained in the introduction, our generalization is motivated by the fact that, in practice, measurements provides outputs from function ϕ that are some averages -or some momentsof ϕ over the actual size of the mechanical device that takes the measurement.

Among the questions raised by GEIM:

• is there an optimal selection for the linear forms σ i within the dictionary Σ ?

• is there a constructive optimal selection for the functions φi ?

• given a set of linearly independent functions { φi } i∈[1,M] and a set of continuous linear forms {σ i } i∈[1,M] , does the interpolant (in the sense of ( 8)) exist?

• is the interpolant unique?

• how does the interpolation process compares with other approximations (in particular orthogonal projections)?

• Under what hypothesis can we expect the GEIM approximation to converge rapidly to ϕ?

In what follows, we provide answers to these questions either with rigorous proofs or with numerical evidences.

The construction of the generalized interpolation functions and linear forms is done recursively, following the same procedure as in the previous subsection, based on a greedy approach, both for the construction of the interpolation linear forms φi and the associated forms selected in the dictionary Σ : The first interpolating function is, e.g.:

φ1 = arg sup ϕ∈F ϕ L 2 (Ω) ,
the first interpolating linear form is:

σ 1 = arg sup σ∈Σ |σ(ϕ 1 )|.
We then define the first basis function as: q1 = φ1 σ 1 ( φ1 )

. The second interpolating function is:

φ2 = arg sup ϕ∈F ϕ -σ 1 (ϕ)q 1 L 2 (Ω) .
The second interpolating linear form is:

σ 2 = arg sup σ∈L(X ) |σ( φ2 -σ 1 ( φ2 )q 1 )|,
and the second basis function is defined as:

q2 = φ2 -σ 1 ( φ2 )q 1 σ 2 ( φ2 -σ 1 ( φ2 )q 1 )
,

and we proceed by induction : assuming that we have built the set of interpolating functions {q 1 , q2 , . . . , qM-1 } and the set of associated interpolating linear forms {σ 1 , σ 2 , . . . , σ M-1 }, for M > 2, we first solve the interpolation problem : find

{ α M-1 j (ϕ)} j such that ∀i = 1, . . . , M -1, σ i (ϕ) = M-1 j=1 α M-1 j (ϕ)σ i (q j ),
and then compute:

J M-1 [ϕ] = M-1 j=1 α M-1 j (ϕ)q j
We then evaluate

∀ϕ ∈ F, ε M (ϕ) = ϕ -J M-1 [ϕ] L 2 (Ω) ,
and define: φM = arg sup

ϕ∈F ε M-1 (ϕ)
and:

σ M = arg sup σ∈Σ |σ( φM -J M-1 [ φM ])| The next basis function is then qM = φM -J M-1 [ φM ] σ M ( φM -J M-1 [ φM ])
.

We finally define the matrix B M such that BM ij = σ i (q j ), and set

X M ≡ span{q j , j ∈ [1, M ]} = span{ φj , j ∈ [1, M ]}.
It can be proven as in in [START_REF] Nguyen | Certified real-time solution of parametrized partial differential equations[END_REF][START_REF] Prud'homme | Reliable realtime solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF][START_REF] Veroy | A posteriori error bounds for reduced-basis approximation of parametrized noncoercive and nonlinear elliptic partial differential equations[END_REF].

Lemma 3 For any M ≤ M max , the set {q j , j ∈ [1, M ]} is linearly independent and XM is of dimension M .The matrix B M is lower triangular with unity diagonal (hence invertible) with other entries ∈ [-1, 1]. The generalized empirical interpolation procedure is well-posed in L 2 (Ω).

In order to quantify the error of the interpolation procedure, like in the standard interpolation procedure, we introduce the Lebesgue constant in the L 2 norm:

Λ M = sup ϕ∈F J M [ϕ] L 2 (Ω) ϕ L 2 (Ω) i.e.
the L 2 -norm of J M . A similar result as in the previous subsection holds Lemma 4 ∀ϕ ∈ F , the interpolation error satisfies:

ϕ -J M [ϕ] L 2 (Ω) ≤ (1 + Λ M ) inf ψM ∈ XM ϕ -ψ M L 2 (Ω)
A (very pessimistic) upper-bound for Λ M is:

Λ M ≤ 2 M-1 max i∈[1,M] q i L 2 (Ω)
Proof. The first part is standard and relies on the fact that, for any ψ ∈ XN then

J M (ψ M ) = ψ M . It follows that ∀ψ M ∈ XM , ϕ -J M [ϕ] L 2 (Ω) = [ϕ -ψ M ] -J M [ϕ -ψ M ] L 2 (Ω) ≤ (1 + Λ M ) ϕ -ψ M L 2 (Ω)
Let us now consider a given ϕ ∈ F and its interpolant

J M [ϕ] = M i=1 α M i (ϕ)q i in dimension M .
The constants α M i (ϕ) come from the generalized interpolation problem:

∀j ∈ [1, M ], σ j (ϕ) = j-1 i=1
α M i (ϕ)σ j (q i ) + α M j (ϕ)j(ψ). We infer the recurrence relation for the constants:

∀j ∈ [1, M ], α M j (ϕ) = σ j (ψ) - j-1 i=1 α i (ψ)σ j (q i ).
Based on the properties of the entries in matrix BM stated in lemma 3, we can obtain, by recurrence, an upper bound for each α

M j (ϕ): ∀j ∈ [1, M ], | α M j (ϕ)| ≤ 2 j-1 ϕ L 2 (Ω) . Then, ∀ϕ ∈ F, ∀M ≤ M max : J M (ϕ) L 2 (Ω) ≤ M i=1 2 j-1 q i L 2 (Ω) ϕ L 2 (Ω) . Therefore: Λ M ≤ 2 M-1 max i∈[1,M]
q i L 2 (Ω) . Note that the norms of the rectified basis function q i verify q i L 2 (Ω) ≥ 1 from the hypothesis done on the norm of the σ i .

Numerical results

The results that we present here to illustrate the GEIM are based on data acquired in silico using the finite element code Freefem [3] on the domain represented on figure 1.

Figure 1: The domain Ω and its mesh.

We consider over the domain Ω ∈ I R 2 the Laplace problem :

-∆ϕ = f, in Ω (9) f = 1 + (α sin(x) + β cos(γπy))χ 1 (x, y)
complemented with homogeneous Dirichlet boundary conditions. Here α, β and γ are 3 parameters freely chosen in given intervals in I R that modulate the forcing term on the right hand side. We assume that the forcing term only acts on a part of Ω named Ω 1 (Ω 1 = support(χ 1 )) and we denote as Ω 2 the remaining part Ω 2 = Ω \ Ω 1 .

The easy observation is that the solution ϕ, depends on the parameters α, β, γ : we plot here one of the possible solutions Figure 2: One of the solutions, we note that the effect of the forcing is mainly visible on domain Ω 1 on the left hand side.

We also note that the restriction ϕ |Ω2 to Ω 2 is indirectly dependent on these coefficients and thus is a candidate for building a set (when the parameters vary) of small Kolmogorov width. This can be guessed if we look at the numerical simulations obtained for three representative choices for α, β, γ

For the GEIM, we use moments computed from the restriction of the solution ϕ(α, β, γ) over Ω 2 multiplied by localized functions with small compact support over Ω 2 . The reconstructed solutions with the GEIM based on only 5 interpolating functions is 10 14 time better than the reconstructed function with 1 interpolating function illustrating the high order of the reconstruction's convergence.

In the next example, we choose a similar problem but the shape of domain Ω 2 is a further parameter In order to get an idea of the Kolmogorov width of the set {ϕ |Ω2 (α, β, γ, Ω 2 ), we perform two Singular Value Decompositions (one in L 2 , the other in H 1 ) over 256 values (approximated again with Freefem) and plot the decay rate of the eigenvalues ranked in decreasing order: the results are shown on figure 5 We note that after about 9 eigenvalues, the finite element error dominates the decay rate of the true eigenvalues. The GEIM is built up again with captors represented as local weighted averages over Ω 2 . The interpolation error is presented on the next figure (figure 6) and we note that the decay rate, measured both in L 2 and H 1 is again quite fast. In order to compare with the best fit represented by the projection, in L 2 or in H 1 , we use the SVD eigenvectors associated with the first M eigenvalues and compare it with J M , for various values of M . This is represented on figure 7.

The very good comparison allow to expect that the Lebesgue constant is much better than what is announced in lemma 4. A computational estimation of Λ M has been carried out: 3 Coupling of deterministic and assimilation methods

Λ M = max i∈[1,256] I M [u i ] L 2 (Ω) u i L 2 (Ω) 0

The framework

Imagine that we want to supervise a process in real-time for which we have a parameter dependent PDE. Assume that the computation of the solution over the full domain Ω is too expensive but we are in a situation where the domain Ω can be decomposed, as before, into two non overlapping subdomains Ω 1 and Ω 2 and that

• Ω 1 is small subdomain but the set of the restriction of the parameter dependent solutions has a large Kolmogorov width.

• Ω 2 is a big subdomain but the set of the restriction of the parameter dependent solutions has a small Kolmogorov n-width

In addition assume that it is possible to get outputs from sensors based in Ω 2 . The GEIM allows to reconstruct accurately the current solution associated to some parameters over Ω 2 and thus is able to build the boundary condition necessary over the interface between Ω 1 and Ω 2 that with the initially given boundary condition over ∂Ω to be the necessary boundary condition over ∂Ω 1 that complement the original PDE set now over Ω 1 and not Ω as is illustrated in the next figures.

The combined approach -numerical results

We take over the numerical frame of the previous section and go further. We want to apply the GEIM to have a knowledge of the solution ϕ |Ω2 and want to use the trace of the reconstruction on the interface to provide the boundary condition, over ∂Ω 1 to the problem

-∆ϕ = f, in Ω 1 f = 1 + (α sin(x) + β cos(γπy))χ 1 (x, y)
derived from [START_REF] Patera | Reduced Basis Approximation and A Posteriori Error Estimation for Parametrized Partial Differential Equations, Version 1.0[END_REF].

The results are presented in figure 11 where both the H 1 error on ϕ |Ω1 and ϕ |Ω2 are presented as a function of M being the number of interpolation data that are used to reconstruct ϕ |Ω2 . This illustrates that the use of the small Kolmogorov width of the set {ϕ |Ω2 } as the parameters vary (including the shape of Ω 2 ) can help in determining the value of the full ϕ all over Ω.

About noisy data

In practical applications, data are measured with an intrinsic noise due to physical limitations of the sensors. In some sense, the noisy data acquired from the sensors are exact acquisitions from a noisy function that we consider to be a Markovian random field with spacial values locally dependent (on the support of the sensor) and globally independent (from one sensor to the others). An extension of the previous development needs therefore to be done in order to take this fact under consideration.

Let us assume that all the sensors are subject to the same noise, i.e. provide averagesor some moments -computed, not from ϕ, but from a random process ϕ ε ≃ N (ϕ, ε 2 ). The norm of the GEIM operator being equal to Λ M the GEIM-reconstruction forms a random process J

M [ϕ ε ] ≃ N (J M [ϕ], Λ 2
M ε 2 ) due to linearity. Even though the Lebesgue constant seems to be small in practice, we would like to use all the data that are available in order to get a better knowledge of ϕ. For the definition of J M we indeed only use M data selected out of a large set of all data. For this purpose, let us consider that, with some greedy approaches, we have determined P independent series of M different captors {σ follows a Gaussian Markov random field of parameters N (J N (ϕ), ǫ 2 λ 2 P ). A realization of this random process could be chosen for an improved estimate of J M (ϕ). Indeed, the law of the error follows N (0, ǫ 2 λ 2 P ) and its variance can be less than the size of the initial noise on the captors (ǫ) provided that Λ (p) N < √ P , ∀1 ≤ p ≤ P , which, from the numerical experiments, seems to be the case.

Conclusions

We have presented a generalization of the Empirical Interpolation Method, based on ad'hoc interpolating functions and data acquired from sensors of the functions to be represented as those that can arise from data assimilation. We think that the GEIM is already interesting per se as it allows to select in a greedy way the most informative sensors one after the other. It can also propose, in case this is feasible, to build better sensors in order to complement a given family of existing ones and/or detect in which sense some of them are useless because redundant. Finally we also explain how noise on the data can be filtered out.

The coupled use of GEIM with reduced domain simulation is also proposed based on domain decomposition technique leading to a small portion where numerical simulation is performed and a larger one based on data assimilation.

We think that the frame presented here can be used as an alternative to classical Bayesian or frequentistic statistic where the knowledge developed on the side for building mathematical models and their simulations can be fully used for data mining (we refer also to [START_REF] Patera | Regression on Parametric Manifolds: Estimation of Spatial Fields, Functional Outputs, and Parameters from Noisy Data[END_REF] and [START_REF] Rozza | Reduction strategies for PDE-constrained optimization problems in haemodynamics[END_REF] for recent contributions in this direction).

Figure 3 :

 3 Figure 3: Three generic solutions restricted on the sub-domain Ω 2 .

Figure 4 :

 4 Figure 4: Two generic solutions when shape of the sub-domain Ω 2 varies.

2 Figure 5 :Figure 6 :

 256 Figure 5: Two SVD ( in L 2 and in H 1 ) of the set of solutions over Ω 2 .

Figure 7 :

 7 Figure 7: Evolution of the GEIM error versus the best fit error, both in L 2 and in H 1 -norms.

Figure 8 :

 8 Figure 8: Evolution of the Lebesgue constant, i.e. the norm of the GEIM operator, both in L 2 and in H 1 .

Figure 9 :

 9 Figure 9: Schematic representation of the reconstruction over Ω 2 .

Figure 10 :

 10 Figure 10: Schematic representation of the recovery over Ω 1 thanks to the knowledge of the interface condition.

Ω 1 Figure 11 :

 111 Figure 11: Reconstructed analysis -error in H 1 -norm over Ω 1 and Ω 2 .

M.

  }, ∀1 ≤ p ≤ P . For each of these series, the GEIM applied to ϕ is noisy and each application providesJ p M [ϕ ε ] ≃ N (J p M [ϕ], Λ p M 2 ε 2). We shall use these P reconstructions by averaging them and expect to improve the variance of the reconstruction.Let λ -1 = 1 Since the P realizations : {J p M [ϕ ε ]} p are independent, then the ran-

Acknowledgements

This work was supported in part by the joint research program MANON between CEA-Saclay and University Pierre et Marie Curie-Paris 6. We want to thank G. Biot from LSTA and G. Pagès from LPMA for constructive discussions on the subject.