Keywords: connectivité fonctionnelle, recalage non-linièaire, modèle deformables, optimisation, machine learning, IRMf functional connectivity, nonlinear registration, deformable models, optimization, machine learning, fMRI networks, and functional connectivity . . . . . . . . . . . . . . subject functional variability . . . . . . . . .

Titre : Amélioration de connectivité fonctionnelle par utilisation de modèles déformables dans l'estimation de décompositions spatiales des images de cerveau

Dedicated to my parents: Florence & Henry Dohmatob.

And to my sister and brothers: Kah, Kingsly, Valy.

-The most wonderful people I know!-3.

Notation

Throughout the manuscript, we shall use the following standard notation and terminology without further explanation.

Notation

Name / synopsis De nition

[[k]]

Integers from 1 to k (inclusive) {1, 2, . . . , k} v r r norm for vectors in a nite-dimensional Euclidean space

       r i | i | r , if 1 ≤ r < ∞, max i | i |, if r = ∞ H
Arbitrary Hilbert space with norm inner product ., . H and norm v → v H := v, v H Normed space containing its Cauchy limits

H (X )
Entropy of a random variable X ∼ p X x p X (x ) log(p X (x ))

MI (X 1 , X 2 ) Mutual Information between two random variables X 1 ∼ p X 1 and Mixed-norm of matrix X ∈ R n×m [ X 1 r , X 2 r , . . . , X n r ] s I n Identity matrix of size n I i j = δ i j , ∀1 ≤ i, j ≤ n

X 2 ∼ p X 2 H (X 1 ) -H (X 2 |X 1 ) N MI (X 1 , X 2 ) Normalized Mutual Information between two random variables X 1 ∼ p X 1 and X 2 ∼ p X 2 I (X 1 , X 2 )/ H (X 1 )H (X 2 ) B r
1 n
Vector of ones of size n

1 i = 1, ∀1 ≥ i ≥ n X † Moore-
       0, if x ∈ C ∞, otherwise σ C Support function of C σ C (x) := sup z∈C x T z dom( f ) E ective domain of f : H → (-∞, +∞] {x ∈ R n | f (x) < +∞} ∂ f (x ) Subdi erential of f at x ∂ f (x) := {v ∈ H | f (z) ≥ f (x) + v T (z -x) ∀z ∈ H } prox f (x ) Proximal operator of f at x arg min p 1 2 p -x 2 2 + f (p) proj C (x ) Orthogonal projection of x onto C arg min p∈C 1 2 p -x 2 2 = prox i C (x) f * Convex conjugate of f f * (x) := sup y x T z -f (z) L F Lipschitz constant of F : H 1 → H 2 inf {C ≥ 0| F (x) -F (y) H 2 ≤ C x -y H 1 ∀x, y ∈ H 1 }
∇ Discrete spatial gradient operator. This de nes a linear operator from R p to R 3p , where p is the number of voxels in the image At a voxel j, the spatial gradient of an image w is a vector ∇w(j) := [∇ x w(j), ∇ w(j), ∇ z w(j)], ∀w ∈ R p ∆ Discrete spatial image Laplacian operator -∇ T ∇ ∈ R p×p

∇ ρ

The identity-augmented version of the discrete spatial gradient operator

∇ ρ w := [(1 -ρ)∇w, ρw] ∈ R 4p , ∀w ∈ R p ∆ ρ
Laplacian operator corresponding to the identity-augmented spatial gradient operator ∇ ρ . This de nes a linear operator from R p to R 4p

ρ 2 I + (1 -ρ) 2 ∆ ∈ R p×p Lap(w)
Laplacian regularization of a 3D image w 1.3 Organization of the manuscript . . . . . . . . 14

Context

A of the neurosciences is to understand the structure, function, and variability of the human brain, and how these give rise to the complex high-level behavior of human beings. One is typically interested in questions such as:

• Which parts of the brain are in-charge of processing mathematical formulae as opposed to ordinary natural language ?

• Which parts of the brain increase/decrease their activity when the brain is at rest ?

• What are the neuro-biological markers of neurological or psychiatric mental illness ?

• How does the brain structure (sulci, gyri, etc.) and function change during aging ? etc.

• How do the language-responsive regions of one subject compare with that of another ? Can they be registered anatomically ?

• How are the di erent motor or cognitive functions (language, emotion, etc.) distributed over the brain, in terms of regions and networks of regions ?

• How are numbers represented and manipulated in the brain ?

• How does the brain and behavior change under the attack of a disease (e.g schizophrenia or a neuro-degenerative disease)

Note that this list is by no way exhaustive.

In the last three decades, mapping brain functional connectivity from functional Magnetic Resonance Imaging (MRI) data has become a very active eld of research. However, analysis tools are limited and many important tasks, such as the empirical de nition of brain networks, remain Dendrites Cell body Axon Figure 1.1: Views of the brain at di erent levels of detail. The brain is composed of (spatially connected) regions and such regions are in turn composed of populations of neurons. Left: Simpli ed view of a neuron. A neuron (there are many types) has a cell body called the soma, many regions for receiving information from other neural cells called dendrites, and often an axon (nerve ber) for transmitting information to other cells (an axon can be longer than 1 meter in humans). The information in the axon is transmitted through an electrical signal called action potential, which is based on the electrical properties of the neuronal membrane. Adapted from http:// commons.wikimedia.org/. Right: Each region is associated with a particular function such as sensory areas (e.g. visual cortex, auditory cortex) that receive and process information from sensory organs, motors areas (e.g. primary motor cortex, premotor cortex) that control the movements of the subject, and associative areas (e.g. Broca's area, Wernicke's area) that process the high-level information related to language production and understanding or the Intra Parietal Sulcus -IPS-that processes spatial information. Adapted from http://agaudi.files. wordpress.com/.

di cult due to the lack of a good framework for the statistical modeling of the data used to de ne these networks.

Objectives. The goal of this PhD thesis is to develop new statistical methods for studying inter-subject variability (eg. amplitude of activation, size of activation clusters, topography of activation maps, etc.), the prime goal being to improve the analysis of functional connectivity in the human brain at the population level. It turns out that these concerns naturally lead to problems related to data-driven extraction of functional atlases, multivariate models for brain decoding and segmentation, and inter-subject registration of functional MRI images.

Sketch of contributions

During the preparation of this PhD project, I have authored and co-authored a number of papers in conferences and journals (including NIPS, ICASSP, MICCAI, Frontiers in Neuroscience, etc.). A complete least of my publications can be found on my Google scholar page https://scholar. google.fr/citations?user=FDWgJY8AAAAJ&hl=fr. In gures,

• Total citations ≥ 194.

• Total papers (including co-authored papers) ≥ 15.

• h index ≥ 4.

• 110 index ≥ 3.

Below, I have roughly classi ed my main contributions under their respective sub-elds of relevance. Viz,

• Sparsity and spatial regularization: [Dohmatob et al., 2014], [Dohmatob et al., 2015b], [Abraham et al., 2014], [Eickenberg et al., 2015], [Pellé et al., 2016] • Registration of brain images: [Dohmatob et al., 2016a] • Optimization: [Dohmatob et al., 2015a], [Varoquaux et al., 2015], [START_REF] Dohmatob | A simple algorithm for computing Nash-equilibria in incomplete information games[END_REF] • Modeling inter-subject functional variability: [Dohmatob et al., 2016b] • Neuroscience: [START_REF] Rahim | Integrating Multimodal Priors in Predictive Models for the Functional Characterization of Alzheimer's Disease[END_REF], [Thirion et al., 2014] There are also a number of preprints currently being prepared for journal publication:

• Sparsity and spatial regularization: "Structured penalties for brain decomposition and decoding: a uni ed view"

• "Inter-subject registration of functional images: do we need anatomical images ?"

• "Enhanced prediction of task-based activation maps from resting-state data"

Organization of the manuscript

In this report, I shall present a selection 1 of the work I have done during 1 For example, I shall not talk on excursional work I did on algorithmic non-cooperative game theory [START_REF] Dohmatob | A simple algorithm for computing Nash-equilibria in incomplete information games[END_REF].

the preparation of my PhD project. This selection will be centered around

• Part I: General preliminaries on neurosciences and neuro-imaging methodology =⇒ chapter 2.

• Part II: Structured penalties for brain decoding =⇒ chapters 3, 4, 5, 6, 7.

• Part III: Functional inter-subject variability =⇒ chapters 9, 10, 8, 11

• Conclusion: Summary and concluding remarks =⇒ chapter 12.

My precise contributions in these domains will be comprehensively outlined as we proceed.

Bibliography

Alexandre Abraham, Elvis Dohmatob, Bertrand Thirion, Dimitris Samaras, and Gael Varoquaux. Region segmentation for sparse decompositions: better brain parcellations from rest fMRI. In Sparsity Techniques in Medical Imaging, 2014.

Elvis Dohmatob. A simple algorithm for computing Nash-equilibria in incomplete information games. In OPT2016 -NIPS workshop on optimization for machine learning, 2016.

Elvis Dohmatob, Alexandre Gramfort, Bertrand Thirion, and Gaėl Varoquaux. Benchmarking solvers for TV-l1 least-squares and logistic regression in brain imaging. In PRNI. IEEE, 2014.

Elvis Dohmatob, Michael Eickenberg, Bertrand Thirion, and Gael Varoquaux. Local Q-Linear Convergence and Finite-time Active Set Identication of ADMM on a Class of Penalized Regression Problems. In ICASSP 2016, 2015a.

Elvis Dohmatob, Michael Eickenberg, Bertrand Thirion, and Gaël Varoquaux. Speeding-up model-selection in GraphNet via early-stopping and univariate feature-screening. In Pattern Recognition in NeuroImaging (PRNI), 2015 International Workshop on. IEEE, 2015b.

Elvis Dohmatob, , Gaël Varoquaux, and Bertrand Thirion. Inter-subject highres EPI-to-EPI direct nonlinear registration outperforms classical T1based method. In Annual meeting of the Organization for Human Brain Mapping -2016, 2016a.

Elvis Dohmatob, Arthur Mensch, Gaël Varoquaux, and Thirion Bertrand. Learning brain regions via large-scale online structured sparse dictionary-learning. In NIPS, 2016b. N as a distinguished data-acquisition and set of analysis techniques for probing and observing brain activity. Acquisition of the data goes hand-in-hand with statistical analysis methods for analyzing the data, in view of making speci c quanti able claims. These techniques operate at a scale much coarser than that of the neuron: one is interested physiological e ects which are ultimately aggregates of activity over large population of neurons (see Fig. 1.1).

In this introductory chapter, I review the relevant theory su cient to situate my own work in a larger scienti c context. Section 2.1 will focus on imaging the human brain and preprocessing of the collected data and also classical methodologies for analyzing the data. Section 2.2 will present another celebrated way of probing brain function, namely resting-state fMRI -or the study of background spontaneous brain activity at rest.

Functional magnetic resonance imaging

Human neuroimaging consists in acquiring ex-vivo (non-invasively) image data from normal and diseased human populations. Several types of functional imaging techniques have been developed. Electro-encephalography (EEG) and magneto-encephalography (MEG) measure the super cial corti-cal neural activity of the brain with a high temporal resolution. Functional magnetic resonance imaging (fMRI) [Ogawa et al., 1990a,b] uses strong magnetic elds to measure changes in oxygen ow in the brain that correlates with synaptic activity in the brain. This technique yields information on brain structure, variability, and function at high spatial resolution. Finally, invasive techniques have been developed such as positron emission tomography (PET) that relies on a radioactive tracer to track glucose consumption.

Figure 2.1: Imaging modalities for the brain. Left: The di erent imaging modalities for brain mapping. MRI and functional MRI have the unique property to yield high-resolution information while being minimally invasive. Unlike other modalities, MRI allows whole brain imaging. Right: Typical example of T1 / anatomical MRI (top), preprocessed Di usion-Weighted (DW) MRI middle and fMRI bottom images, presented in axial views. These images are from the Neurospin 3T scanner. For the DW-MRI image, the main direction of water di usion is colorcoded: green for antero-posterior di usion, red for lateral di usion, blue for vertical di usion. The functional image has been analyzed to yield the regions activated in an auditory task. Adapted with permission from [START_REF] Thirion | Structural and probabilistic methods for group analysis in functional neuroimaging[END_REF].

The BOLD signal

When a brain area is solicited, the brain res chemical signals to report the consumption of oxygen and sugar. Nearby blood capillaries dilate to increase the quantity of owing blood and provide these resources. This phenomenon is called the haemodynamic response. As a result, we expect a higher concentration of oxygenated hemoglobin in a given brain area soon after its activation. fMRI imaging can be used to measured this e ect, called the BOLD (Blood Oxygen-Level-Dependent) signal [Ogawa et al., 1990a,b], at a spatial resolution of 1.5 to 3mm, and a temporal resolution of 1-3s, typically. This yields a spatially resolved image of brain functional networks that can be modulated either by speci c cognitive tasks or appear as networks of correlated activity. This method is subject to several physical and physiological noises. First, some artifacts may be induced by radio transmitters or other equipment. Then, spurious activations are naturally introduced by arteries present in the brain, heart beats and breathing movements. Finally, the brain can be shifted if the subject makes large movements in the scanner.

Preprocessing and analysis of fMRI data

Raw fMRI images are not intepretable with bare eyes. In particular because we are interested in small signal co-variation between voxels 1 and 1 Voxel stands for volume element.

It refers to a point in a 3D image, just as pixel refers to a point in a 2D image.

not by the values themselves. The human eye, however, is good at perceiving global artifacts in the data such as movements, ghost or scanner coils.

Quality assessment of preprocessed fMRI data is done by eye and by relying on dedicated medical imaging software. In order to prepare the data for further statistical analysis, some preprocessing steps are required. Below, we outline the main ones. Viz, Data acquisition. The resolution of fMRI is usually between 1mm 3 and (3mm) 3 . In a single 3D scan, the brain represents 10 4 to 10 6 voxels. A run contains usually from 100 to 1000 scans. Functional MRI scans are acquired by slices, usually in the axial direction. The time required to acquire one slice is called echo time (TE) and is in the order of tens of milliseconds. The time required to acquire a whole 3D volume is called repetition time (TR) and is in the order of seconds. Typical values for a 3D volume of 60 slices are TE=33ms and TR=2s for a 3T (Tesla) scanner.

Motion-correction and coregistration to the anatomy. Head movement has a big impact on fMRI. A movement with an amplitude higher than the voxel resolution (i.e. 2 to 3mm) can shift the signal of the entire brain. Moreover, the worst impact of motion is in ow e ects, i.e. artefactual signals. In the scanner, the head of the subject is xed using cushion pads to avoid movements and the subject is asked to stay as still as possible. Yet, it is impossible to completely avoid head movement. In order to mitigate the effect of movement, the 3D scans are realigned on a reference scan -usually the one in the middle of the sequence-using rigid body transformation (translation and rotation, without change of scale). This is usually followed by an a ne registration of the motion-corrected images to the anatomical (T1) image of the subject, in view of subsequent inter-subject preprocessing and analysis, like registration onto a group template (more on this later).

Slice-timing correction. As stated before, brain slices are not acquired at the same time. This introduces a shift in the haemodynamic response associated to each of them. The problem can be solved by interpolating the signal of each slice so that all of them can be considered as acquired at the same time. [START_REF] Sladky | Slice-timing e ects and their correction in functional MRI[END_REF] showed that depending on repetition time and paradigm design, slice-timing e ects can signi cantly impair fMRI results and slice-timing correction methods can successfully compensate The bottom plot show estimated motion parameters for the subjects-head motion.

Here all movements are well below 0.5mm, which is generally considered as ne. The preprocessing and plots displayed here were done using Pypreprocess https://github.com/ neurospin/pypreprocess, an open-source Python wrapper built on standard toolkits like SPM [Friston et al., 1994].

for these e ects and therefore increase the robustness of subsequent data analysis.

Registration of brain data into a common reference space. Each brain is of di erent size and shape. In order to compare brain activations across several individuals, we need to normalize them by registration to a common template [Friston et al., 1995[START_REF] Ashburner | Uni ed segmentation[END_REF], Ashburner, 2007, Klein et al., 2009]. This template can be a reference template used in the community (MNI for example). It is also possible to compute a template directly from the data. Once a template is chosen, for each subject, we perform two successive registrations. First, the anatomical scan acquired in the subject is registered to the MNI template. Then, the fMRI data are registered to the anatomical scan. After that, the two transformation matrices are combined in order to normalize the fMRI data to the template 2 . Estimation of the deformations necessary to warp a subject's 2 In chapter 8, we study the possibility of by-passing the anatomical image, when normalizing functional data.

brain anatomy onto a template is usually done alongside the classi cation of individual voxels into di erent classes: white matter (wm), grey matter (gm), and cerebro-spinal uid (csf) producing so-called tissue probability maps (TPMs) [START_REF] Ashburner | Uni ed segmentation[END_REF].

Statistical analysis of brain data

Forward inference made on fMRI data (e.g. prediction of brain activation from the stimuli) can be conceptualized as the encoding of perceptual, motor or cognitive parameters into brain signals. The inverse model, that predicts behavioral data from brain activation is called decoding, and will be the subject matter of chapter 3. Two main paradigms allow to experimentally study brain signals: either we study them in controlled condition on a particular task -this is the task paradigm-or we study the spontaneous activity of the brain in order to uncover its organization: this is the resting-state paradigm.

The task paradigm and the general linear model

Using an experimental design, it is possible to relate the BOLD signal with speci c tasks performed by the subject. For example, a sound can be played in the left or the right ear of the subject. By comparing brain activation between resting state and when the sound is playing, we can isolate the auditory cortex of the brain. Statistically, we do that by crafting a design matrix corresponding to the experiment: one column of the matrix represents an ideal response to one of the presented conditions [Friston et al., 1994]. Columns corresponding to known artifacts of the BOLD signal, such as heart beats or movements, can be added in the design matrix in order to regress out the part of the signal related to them. We then use a general linear model (GLM) to recover the brain maps corresponding to each of the columns in the design matrix X ∈ R T ×k , where k is the number of conditions and T is the number of time points (times of repetition -TR). It is then supposed that for each voxel , the measured BOLD signal y ∈ R T is a linear combination of the columns of X, i.e of the experimental conditions, that is

y = Xβ + ϵ , ( 2.1)
where β ∈ R k are regression coe cients and ϵ = (ϵ ,1 , . . . , ϵ ,T ) ∈ R T is a non-iid vector of normally distributed noise. Such a problem is well-posed and weighted least-squares (WLS) are used to obtain a solution 3 , to obtain 3 There are usually more time points than experimental conditions, and so the design matrix is full-rank.

β = X † y . Stacking these coe cients across all voxels per-brain correspond to k so-called β-maps [Friston et al., 1994]. For a given combination of experimental conditions 4 c = (c 1 , c 2 , . . . , c k ) ∈ R k , one can compute a 4 For example, for k = 3 conditions, one may be interested in take c = (1, -1, 0), meaning we wish the nd the e ect of the rst condition relative to the second, or c = (1, -1/2, -1/2) corresponding to the e ect of the rst condition w.r.t the average e ect. statistic

t ,c := c T β ar (ϵ )c T (X T X) -1 c . (2.2)
Under the null hypothesis that the e ect were are interested in is zero, i.e

H 0 : c T β = 0, (2.3)
the above statistic is student-t distributed with Tk degrees of freedom, and one can analytically obtain p-values and con dence intervals for inference. Projecting these values unto the brain (one value per voxel) yields a so-called activation map. Such maps are the main output of any forward analysis in task-based fMRI studies. Subsequent statistical inference su ers from heavy multiple comparison issues in these so-called mass-univariate methods. The problem is further confounded by the fact that there are correlations between neighboring voxels, leading to situation where the Bonferoni and similar correction procedures, usually used to deal with these issues, may be too conservative and destroy the the sought-for e ects. An alternative is to use multi-variate methods which directly model the spatial interactions between the voxels. Such methods will be the subject of chapter 3.

Resting-state fMRI, brain networks, and functional connectivity

Resting-state fMRI -or rsfMRI for short-uses the same acquisition method as task fMRI. However, instead of giving a particular task to the subject, they are asked to let their mind wander without sleeping. By studying this background activity of the brain, it is possible to uncover its underlying organization [START_REF] Raichle | Two views of brain function[END_REF]. Unlike the techniques described previously where the aim was to localize regions of activation for a given set of conditions, in functional connectivitly analysis were are interested in infering connections between such regions. Depending on the protocol, the subject can be asked to keep eyes closed or to contemplate a xation cross. The xation cross prevents random eye movements and helps the subject not to sleep. In rsfMRI, we do not study the signal of each voxel itself but the interactions between the brain voxels. In particular, we study the functional connectivity of the brain, i.e. the similarity of activation patterns between brain regions that share a common functional role. Since there is no design matrix in rest fMRI, one must be careful to properly regress out physiological noises or spurious correlations may appear between brain regions, in particular longitudinally [START_REF] Power | Spurious but systematic correlations in functional connectivity MRI networks arise from subject motion[END_REF][START_REF] Van Dijk | The in uence of head motion on intrinsic functional connectivity MRI[END_REF]. A rst approach of functional connectivity is the voxel-to-voxel approach in which the similarity is measured between each pair of voxels. This method is not only computationally expensive, given the number of voxels in the brain, but it is also unfounded from the statistical standpoint: it requires the estimation of millions of parameters (one for each voxel pair), much more than the number of observations supports. As a consequence, some form of dimensionality reduction -a feature selection or extraction-is necessary to study connectivity.

Inter-subject functional variability

As noted in [Thirion et al., 2007, Thyreau et al., 2012[START_REF] Xu | Modeling inter-subject variability in fMRI activation location: A Bayesian hierarchical spatial model[END_REF], the inter-subject variability in GLM results (see Fig. 2.6) is not due to misregis- tration, but intrinsic subject di erences with a more physiological nature: the size of e ects and the anatomical localization are subject-speci c. Also, [Tavor et al., 2016] used dual regression [START_REF] Filippini | Distinct patterns of brain activity in young carriers of the APOE-ε4 allele[END_REF] to provide quantitative evidence that inter-subject di erences in task-based brain activations are largely physiological -in contrast to being driven by subjects' brain morphological di erences.

Chapters 9 and 11 will present generative models for understanding intersubject variability at the functional level.

Introduction to brain decoding

As already discussed in chapter 2, functional brain imaging provides a distinctive opportunity to study brain functional architecture, while being minimally invasive, and is thus well-suited for the challenging study of the spatial layout of neural coding. Di erent modalities exist, each one having speci c spatial and temporal resolutions; among them Functional Magnetic Resonance Imaging (fMRI) [Ogawa et al., 1990a,b] has emerged as a fundamental modality for brain imaging, striking a good balance between spatial and temporal resolution. fMRI images are pre-processed, and modeled through a general linear model (GLM), that takes into account the di erent experimental conditions and the dynamics of the haemodynamic response in the design matrix. The resulting model parameters, a.k.a. activation maps, represent the in uence of the di erent experimental conditions on local fMRI signals.

The standard approach. The standard approach used for analyzing these activation maps is called classical inference, and relies on a mass-univariate statistical tests (one for each voxel), yielding the so-called statistical parametric maps (SPMs) [Friston et al., 1994]. Such maps are of particular interest in cognitive neuroscience, as they open the door to localizing the voxels that are signi cantly active for any combination of experimental conditions, and thus are probably implied in the underlying neural code of the cognitive processes. However, this classical inference su ers from multiple comparisons issues. Also, it does not take into account the multivariate structure of the fMRI data. Inverse inference (or "brain reading"). An alternative approach called inverse inference (or "brain-reading") [Dehaene et al., 1998, Cox andSavoy, 2003], has been proposed in order to cope with the limitations of the aforementioned classical inference. Inverse inference relies on a pattern recognition, and aims at decoding the neural code by using machine-learning methods. Based on a set of brain activation maps, inverse inference builds a predictive model that can be used for predicting a behavioral target (age, sex, IQ, etc.) for a new set of images. The prediction accuracy of the model is used as a measure of the quantity of information about the cognitive task shared by the voxels. By construction, this approach is multivariate, and can provide more sensitive analysis than standard statistical parametric mapping procedure [Kamitani 05, Haynes 06] Several methods have been tested for classi cation or regression of activation images (Linear Discriminant Analysis -LDA, Support Vector Machines -SVM, Lasso, Elastic net regression, and many others), but, in this problem, the major bottleneck remains the localization of predictive regions within the brain volume. Additionally, we have to deal with the curse of dimensionality, as the number of features (voxels, regions) is much larger (∼ 10 5 -10 6 ) than the numbers of sample (images) (∼ 10 2 -10 3 ), the latter being limited by the cost of acquisition.

Thus the prediction method may over t the training set and thus not generalize well to new samples.

Sparsity and structure-inducing priors: towards intepretable multi-variate models

To cope with the high dimensionality of the data, the learning method has to be regularized. However, the spatial structure of the image is not taken into account in standard regularization methods, so that the extracted features are often hard to interpret. Sparsity and spatial smoothness inducing priors can be used to perform jointly the prediction of a target variable and region segmentation in multivariate analysis settings. Sparsity can be enforced by penalizing the (sum of) absolute values of the regression coe cients, leading to the so-called Lasso model. Smoothness can achieved in penalizing the spatial gradient of the regression coe cients, to enforce smooth regions ("blobs"). The Total-variation (TV) [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] penalty has proven to be particularly powerful for realizing such e ects. Laplacian regularization is an easier means to this end (because in leads to a di erentiable problem), but have sub-optimal rates for noisy signal recovery [Sadhanala et al., 2016], and the visual e ect is less appealing.

In the context of neuro-imaging, sparsity and smoothness have been compiled to yield regression coe cients which are faithful to known neurobiological organization of the brain, while alleviating the risk of over-tting due to inherently small-sample settings. Speci cally, it has been shown that one can employ priors like TV-1 [Baldassarre et al., 2012, Gramfort et al., 2013], TV-ElasticNet [START_REF] Dubois | Predictive support recovery with TV-Elastic Net penalty and logistic regression: an application to structural MRI[END_REF], and GraphNet [Grosenick et al., 2013] (aka Smooth-Lasso [Hebiri and van de Geer, 2011]) to regularize regression and classi cation problems in brain imaging. TV has also been employed to enhance the estimation of the voxel-wise Hurst exponent 1 , as

1 H := lim N →∞ E(r N /σ N ) log N
, where r N is the empirical range (i.e max value minus min value) of the rst N values in a time-series, and σ N is their standard deviation. For example in 1D, white noise has H = -1/2. a measure of temporal self-similarity in brain dynamics [Pellé et al., 2016].

Notation. We denote by y ∈ R n the targets to be predicted (age, sex, IQ, etc.); the design matrix X ∈ R n×p are the masked (see Fig. 3.2) brain images related to the presentation of di erent stimuli, or other brain acquisition (e.g gray-matter concentration maps from anatomy, etc.). The integer p is the number of voxels, and n the number of samples (images). In brain imaging, n p; typically, p ∼ 10 3 -10 6 (in full-brain analysis), while n ∼ 10 -10 3 (n being limited by the cost of acquisition, etc.). ∇ x will denote the discrete spatial gradient operator along the x-axis, ∇ along the -axis, etc. Thus, at a voxel j, the spatial gradient of an image w is a vector ∇w(j) := [∇ x w(j), ∇ w(j), ∇ z w(j)], ∀w ∈ R p . This de nes a linear operator ∇ ∈ R 3p×p (the discrete 3D spatial gradient operator) from R p to R 3p . For a mixing constant ρ ∈ [0, 1], ∇ ρ ∈ R 4p×p will denote the identity-augmented version of ∇, de ned by ∇ ρ w := [(1 -ρ)∇w, ρw] ∈ R 4p . Figure 3.2: Masking of volumic brain data (4D = 3D space + 1 time or samples) to produce a design matrix required in standard machine learning (clustering, classi cation, regression, etc.). Each 3D volume considered is a sample point. The values of the voxels in this volume that lie in the mask are collected into a feature vector. All these vectors are vertically stacked to produce an nby-p design matrix X, where p is the number of voxels in the mask. The mask can be just the region of the 3D cube occupied by the brain, or a subset of such. In the latter case, this typically corresponds to Region-of-Interests (ROIs) deemed to be interesting for an experiment. The former case is referred to as "full brain", and the mask typically contains up to p = millions of voxels. See [Abraham et al., 2014] for more details.

Prerequisites. Given that this chapter and many others will be quite heavy on proximal calculus and sparse modelling, we would like to suggest the following references as a good starting point for the non-expert reader on these subjects:

• Proximal calculus [START_REF] Patrick | Signal recovery by proximal forward-backward splitting[END_REF], Beck and Teboulle, 2009, Combettes and Pesquet, 2011].

• Sparse modelling [START_REF] Mairal | Sparse Modeling for Image and Vision Processing[END_REF], Bach et al., 2012].

That notwithstanding, we shall endeavour to develop the material bottom-up assuming as much as possible only a bare minimum prerequisite knowledge on very specialized topics.

SpaceNet: sparse structured models for brain data

We now describe the family of structured models which have been proposed for enhanced multivariate analysis in neuro-imaging, namely: Total-Variation (TV) [Michel et al., 2011], TV-1 [Baldassarre et al., 2012, Gramfort et al., 2013], GraphNet [Grosenick et al., 2013, Hebiri andvan de Geer, 2011], TV-ElasticNet [START_REF] Dubois | Predictive support recovery with TV-Elastic Net penalty and logistic regression: an application to structural MRI[END_REF], Sparse Variation [Eickenberg et al., 2015], and social-sparsity [Kowalski et al., 2013, Varoquaux et al., 2016]. These can all be synthesized into a general framework, referred to as SpaceNet [Dohmatob et al., 2015b], as follows (3.1)

The coe cients w de ne a spatial map in over the brain (one value per voxel). The term (y, Xw) is the loss / data-t term. Popular choices include:

(y, Xw) = 1 n n i=1                      1 2 (X T i w -i ) 2 , for least-squares regression log(1 + exp(-i X T i w)), for logistic regression, ( 1
-i X T i w) + ,
for hinge loss (used in SVMs) . . .

In the above general model, P (w) is the penalty term, which simultaneously imposes both sparsity and structure (blobs). The di erent spatial regularization methods that have appeared in neuro-imaging literature can be cast into this correspond to di erent choices of the convex penalty P acting on the extended gradient of the coe cients w. Viz,

P (w) =                            ρ w 1 + 1 2 (1 -ρ) ∇w 2 Fro = j ∈[[p]] ρ|w j | + 1 2 (1 -ρ) (∇w) j 2 2 ,
for GraphNet ,

∇ ρ w 1+2,1 = ρ w 1 + ∇w 2,1 = j ∈[[p]] ρ|w j | + (1 -ρ) (∇w) j 2 , for isotropic TV-1 , ∇ ρ w 1,1 = ρ w 1 + (1 -ρ) ∇w 1,1 = j ∈[[p]] ρ|w j | + (1 -ρ) (∇w) j 1 , for anisotropic TV-1 , ∇ ρ w 2,1 = j ∈[[p]] (∇ ρ w) j 2 ,
for Sparse Variation , . . . (3.2) where • α > 0 is a regularization parameter controls the total amount of regularization;

• ρ (0 < ρ ≤ 1) is a mixing constant between the sparsity-inducing 1 part and the cluster-promoting part of the penalty term. The particular case ρ = 1 corresponds to the usual Lasso. Vanilla TV [Michel et al., 2011] corresponds to TV-1 with ρ = 0.

• The matrix ∇ ρ is the extended discrete gradient operator de ned in Table 1.

Bayesian interpretation of SpaceNet models. The penalties P (w) in (3.1) admit a Bayesian interpretation as a prior on the distribution of the coe cients w p α ,ρ (w) ∝ exp(-P (w)).

(3.3)

For example, the GraphNet [Grosenick et al., 2013, Hebiri andvan de Geer, 2011], the penalties P (w) penalty corresponds to SpaceNet models (3.1) result in brain maps which are both sparse (i.e regression coe cients w are zero everywhere, except at predictive voxels) and structured (blobby). See Fig. 3.4. The superiority of such methods over methods without structured priors like the Lasso, ANOVA, Ridge, SVM, etc. for yielding more intepretable maps and improved prediction scores is now well established. See for example [Baldassarre et al., 2012, Gramfort et al., 2013]. These priors are fast becoming popular for brain decoding and segmentation. Indeed, they leverage a feature-selection function (since they limit the number of active voxels), and also a structuring function (since they penalize local di erences in the values of the brain map). For example, see Fig. 3.6. Also, such priors produce state-of-the-art methods for automatic extraction of functional brain atlases [START_REF] Abraham | Extracting brain regions from rest fMRI with Total-Variation constrained dictionary learning[END_REF].

Submodular interpretation of TV. We note that anisotropic TV penalty (3.2) on an arbitrary (undirected) graph G = (V , E) is the Lovasz extension of the cut-function F : 2 V → N, S → "number of edges between S and V \ S", de ned by

F L (x ) := E λ∼U ([0,1]) [F ({ ∈ V |x ≥ λ})], for all x ∈ [0, 1] #V .
Thus anisotropic TV minimization can be seen as a graph-cut problem for which e cient algorithms exist [Bach, 2013, Landrieu and[START_REF] Landrieu | Cut Pursuit: Fast Algorithms to Learn Piecewise Constant Functions[END_REF].

Methods

The SpaceNet model leads to di cult non-smooth mathematical optimization problems making their implementation and practical usability challenging. [Dohmatob et al., 2014] benchmarked a rich variety of cutting-edge solvers for such problems, and gave crucial recommendations on how to e ectively implement these algorithms in practice. In these benchmarks, the FISTA algorithm emerged as the go-to algorithm for the TV-L1 problem [Dohmatob et al., 2015a]. These hints have been carefully used in implementing SpaceNet. Also as a preprocessing step, we use univariate featurescreening (ANOVA) to eliminate voxels which are irrelevant to the learning problem, thus reducing the size of the problem. As a result the implementation of SpaceNet is fast, robust, and automatically sets its hyper-parameters (internal cross-validation). All these technical details will be properly presented in the next few chapters.

Cross-validation

Cross-validation (e.g see [START_REF] Stone | Cross-Validatory Choice and Assessment of Statistical Predictions[END_REF]) is a technique used to protect against over tting in a predictive model, particularly in a case where the amount of data may be limited. In cross-validation, you make a xed number of folds (or partitions) of the data, run the analysis on each fold, and then average the overall error estimate. This gives an (asymptotically) unbiased estimate of the true generalization error of the model. Two major types of cross-validation are K-Fold and Leave-One-Out (LOO).

K-Fold cross-validation. One iteration of the K-fold cross-validation is performed in the following way: First, a random permutation of the sample set is generated and partitioned into K subsets ("folds") of about equal size. Of the K subsets, a single subset is retained as the validation data for testing the model (this subset is called the "testset"), and the remaining K -1 subsets together are used as training data ("trainset"). Then a model is trained on the trainset and its accuracy is evaluated on the testset. Model training and evaluation is repeated K times, with each of the K subsets used exactly once as the testset. The case of a 5-fold cross-validation with 30 samples is illustrated in Fig. 8.4.

Leave-One-Out cross-validation. As the name suggests, leave-one-out cross-validation involves using a single sample from the original sample set as the validation data, and the remaining samples as the training data. This is repeated such that each sample in the sample set is used exactly once as the validation data. This is the same as K-fold cross-validation where K is equal to the number of samples in the sample set. In LOO, there is no need in generating random permutations and in repeating it, because the training and validation datasets for each of the folds are always the same, and therefore the result of the accuracy estimation is determined.

Model-selection via cross-validation. One can instrument cross-validation to tune the hyper-parameters of a model like SpaceNet (7.1), by selected the con guration of model parameters with least cross-validation error. The number of models tted is proportional to the size of the parameter grid -i.e exponential in the number of parameters to tune-and therefor can become prohibitive in case there are many free hyper-parameters in the model. Also, since some of the data has to be set aside for validation, cross-validation in very small sample settings (e.g e few tenths, as is the in some neuroimaging experiments) may be troublesome as the error estimates then have very high variance. A reasonable alternative in such situations are SURE (short for Stein's Unbiased Risk Estimator [START_REF] Charles | Estimation of the mean of a multivariate normal distribution[END_REF])-based methods, which are applicable whenever a procedure for obtaining (an unbiased estimate of) the number of degrees of freedom of the model is available. This is the case with the models like the ElasticNet and GraphNet [Hebiri and van de Geer, 2011]. Recently, [Deledalle et al., 2014] has proposed a SURE-like technique for structured models with many hyper-parameters. (b) L × M grid over which to search for optimal conguration in a model with two hyper-parameters α and ρ. For a model like SpaceNet (7.1), the grid is is constrained to verify 0 ≤ α L < . . . < α 1 = α max and 0 ≤ ρ M < . . . < ρ 1 ≤ 1, with L = 10 and M = 3 typically, given a total of LM = 30 models to compare. In chapter 5, we show how early-stopping and other heuristics can be used to make the total cost much more e ective than just tting LM models in a CV loop.

(α 1 , ρ 1 ) (α 1 , ρ 2 ) . . . (α 1 , ρ M ) (α 2 , ρ 1 ) (α 2 , ρ 2 ) . . . (α 2 , ρ M ) (α 3 , ρ 1 ) (α 3 , ρ 2 ) . . . (α 3 , ρ M ) . . . . . . . . . . . . . . . (α L , ρ 1 ) (α L , ρ 2 ) . . . (α L , ρ M ) (b)

How SpaceNet compares against classical unstructured models

Classification. We compared SpaceNet (TV-L1 and GraphNet / Smooth-Lasso priors) with an SVM (Support Vector Machine) on the visual-recognition dataset [Haxby et al., 2001]. This dataset consists of 6 subjects with 12 runs per subject. In each run, the subjects passively viewed images of eight object categories, grouped in 24-second blocks separated by intermittent rest periods. This experiment is a classi cation task: predicting the object category. The design matrix is made of time-series from the full-brain mask of p = 23, 707 voxels over 216 TRs (Repetition Times), of a single subject (subj1). 126 TRs were used for training all the models, whilst testing was done on 90 left-out TRs. The results are depicted in Figures 3.6 and 3.7.

Regression. In [Gramfort et al., 2013], the authors compared several models on a dataset in which subjects were presented with mixed (gain/loss) gambles, and decided whether they would accept each gamble [Jimura and Poldrack, 2012]. No outcomes of these gambles were presented during scanning, but after the scan three gambles were selected at random and played for real money. The prediction task here is to predict the magnitude of the gain and thus a regression on a continuous variable. The full dataset of 16 subjects with 48 3D scans each, making up for a total of n = 768 samples with approximately p = 3.3 × 10 4 voxels. The prediction here is inter-subject: the estimator learns on some subjects and predicts on left out subjects. The results are shown in Fig. 3.8.

Figure 3.6: The gure shows results of comparing the SpaceNet models TV-1 and Graph-Net against an SVM (Support Vector Machine) classi er on the visual-recognition dataset [Haxby et al., 2001] As can be seen from the gure, SpaceNet priors (TV-1 , GraphNet/Smooth-Lasso, etc.) yield stable and more intepretable maps by enforcing smoothness on the coe cients while segmenting predictive regions (blobs) from noisy background. [Haxby et al., 2001]. We see that the highly structured maps produced by SpaceNet models (3.1) (e.g see Fig. 3.6) are not at the expense of model accuracy.

Conclusion

We have presented SpaceNet, a family of priors for brain decoding that enforce both sparsity and structure, leading to better prediction scores and intepretable brain maps. We believe that such priors will become commonplace in future. In the next few chapters, we open the "black-box" and develop from ground-up, the details of such models, including their practical implementation on a computer.

Bibliography

Alexandre Abraham, Elvis Dohmatob, Bertrand Thirion, Dimitris Samaras, and Gael Varoquaux. Extracting brain regions from rest fMRI with Total-Variation constrained dictionary learning. In MICCAI. 2013. T S N introduced in equations (3.1) of chapter 3 lead to superior estimators compared to classical estimators (Ridge regression, SVM, etc.) without spatial penalization, they are considerably harder to optimize than these classical models. Indeed, the corresponding optimization problems is non-separable in the model coe cients, and except for the case of GraphNet [Hebiri andvan de Geer, 2011, Grosenick et al., 2013] and social-sparsity [Kowalski et al., 2013, Varoquaux et al., 2016], the penalty term P (w) is neither smooth nor proximable 1 . For the penalty to 1 A function f is said to be proximable if its operator prox γ f is easy to compute. This is the case for p -norms (with p ≥ 1, to ensure convexity) and indicator functions of simple closed convex sets like balls, simplexes, half-spaces, etc.

fully exercise its structuring e ect on the maps, this optimization problem must be solved to a good tolerance resulting in a computational challenge. Lack of good solver and explicit control of tolerance can lead to brain maps and conclusions that re ect properties of the solver more than of model coe cients, as illustrated in Fig. 4.1.

Solving TV-L1 regularized problems

The optimization problem (3.1) is very challenging: it is non-smooth (except in the case of Laplacian regularization), non-separable and heavily illconditioned. For the penalty to fully exercise its structuring e ect on the maps, this optimization problem must be solved to a good tolerance resulting in a computational challenge. In [Dohmatob et al., 2014], we did an extensive study of all solvers applicable to the problem in TV-1 special case (which happens to be the most di cult scenario). Our results outlined the best strategy: a double FISTA loop, where the inner loop computes the proximal operator of the penalty term, with approximate precision on the duality-gap. This was further re ned and implemented in [Varoquaux et al., 2015].

Figure 4.1: TV-1 maps for the facehouse discrimination task on the visual recognition dataset. Note that the stopping criterion is dened as a threshold on the energy decrease per one iteration of the algorithm, and thus differs from the tolerance displayed in gure 4.1. This gure shows the importance of convergence for problem (3.1), and motivates the need for fast solvers for SpaceNet priors, especially the non-smooth ones like TV-1 and Sparse Variation. See [Dohmatob et al., 2014] for details.

The algorithms

ISTA/FISTA. ISTA [Daubechies et al., 2004], and its accelerated variant FISTA [Beck and Teboulle, 2009a], are proximal gradient approaches: the go-to methods for non-smooth optimization. In their seminal introduction of TV for fMRI, [Michel et al., 2011] relied on ISTA. The challenge of these methods for TV is that the proximal operator of TV cannot be computed exactly; we approximate it in an inner FISTA loop [Beck andTeboulle, 2009b, Michel et al., 2011]. Here, for all FISTA implementations we use the faster monotonous FISTA variant [Beck and Teboulle, 2009b]. We control the optimality of the TV proximal via its dual gap [Michel et al., 2011] and use a line-search strategy in the monotonous FISTA to decrease the tolerance as the algorithm progresses, ensuring convergence of the TV-1 regression with good accuracy. See [Dohmatob et al., 2014, Varoquaux et al., 2015].

ISTA/FISTA with backtracking. A key ingredient in FISTA's convergence is the Lipschitz constant L ∇ , of the derivative of smooth part of the objective function . The tighter the upper bound used for this constant, the faster the resulting FISTA algorithm. In FISTA, the main use of L ∇ is the fact that: for any stepsize 0 < t ≤ 1/L ∇ and for any point z,

(p t (z)) ≤ (z) + r T t ∇ (z) + 1 2t r t 2 2
, where

p t (z) := prox α t P (z -t∇ (z)) and r t := p t (z) -z (4.1)
In least-squares regression, L ∇ is precisely the largest singular value of the design matrix X. For logistic regression however, the tightest known upper bound for L ∇ is X X T , which performs very poorly locally (i.e, stepsizes ∼ 1/L ∇ are sub-optimal locally). A way to circumvent this di culty is backtracking line search [Beck and Teboulle, 2009a], where one tunes the stepsize t to satisfy inequality (4.1) locally at point z.

ADMM: Alternating Direction Method of Multipliers. ADMM is a Bregman Operator Splitting primal-dual method for solving convex-optimization problems by splitting the objective function in two convex terms which are functions of linearly-related auxiliary variables [Boyd et al., 2010]. ADMM is particularly appealing for problems such as TV regression: using the variable split z ← ∇w, the regularization is a simple 1 / 2 norm on z for which the proximal is exact and computationally cheap. However, in our settings, limitations of ADMM are:

• The w-update involves the inversion of a large p-by-p ill-conditioned linear operator (precisely a weighted sum of X T X, the laplacian ∆, and the identity operator).

• The dual stepsize parameter ν in the penalization of the split residual 1 2 ν z -∇w 2 2 is hard to set (this is still an open problem), and though under mild conditions ADMM converges for any value of ν, the convergence rate depends on ν . In chapter 7, we study the rate of convergence of ADMM on the kinds of penalized least squares regression problem considered in this manuscript, and derive some theoretical results.

Primal-Dual algorithm of Chambolle and Pock [Chambolle and Pock, 2011]. This scheme is another method based on operator splitting. Used for fMRI TV regression by [Gramfort et al., 2013], it does not require setting a hyperparameter. However it is a rst-order single-step method and is thus more impacted by the conditioning of the problem. Note that here we explore this primal-dual method only in the squared loss setting, in which the algorithm can be accelerated by precomputing the SVD of X [Gramfort et al., 2013] .

HANSO [START_REF] Lewis | Nonsmooth Optimization via BFGS[END_REF]. a modi ed LBFGS scheme based on gradient sampling methods [Burke et al., 2005] and inexact line-search. For non-smooth problems as in our case, the algorithm relies on random initialization, to avoid singularities with high probability. Here, we used the original authors' implementation.

Uniform approximation by smooth convex surrogates. The 1 norm (resp. TV semi-norm) is di erentiable everywhere with gradient

w j /|w j | j ∈[[p]] (resp. -div(((∇w) j / (∇w) j 2 ) j ∈[[p]]
))), except when some voxels are inactive with w j = 0 (resp. (∇w) j = 0), corresponding to black spots (resp. edges). A convenient approach (see for example [Bobin et al., 2011, Nesterov, 2005a,b, Beck and Teboulle, 2012]) for dealing with such singularities is to uniformly approximate the o ending function with smooth surrogates that preserve its convexity. Given a smoothing parameter µ > 0, we de ne smoothed versions of 1 and TV:

w 1,µ := j w 2 j + µ 2 , w TV,µ := j (∇w) j 2 2 + µ 2 (4.2)
These surrogate upper-bounds are convex and everywhere-di erentiable with gradients that are Lipschitz-continuous with constants 1/µ and ∇ 2 (1/µ) = 12/µ respectively. They lead to smoothed versions of problem (3.1):

ŵµ := arg min w {E µ (w) := (w) + α P TV-L1,µ (w)}, (4.3)
where

P TV-L1,µ (w) := ρ w 1,µ + (1 -ρ) w TV,µ .
To solve (3.1), we consider problems of the form (4.3) with µ → 0 + : we start with a coarse µ (= 10 -2 , e.g) and cheaply solve the µ-smoothed problem (4.3) to a precision ∼ µ using a fast iterative oracle like the LBFGS [START_REF] Zhu | L-BFGS-B-FORTRAN SUBROUTINES FOR LARGE-SCALE BOUND CON-STRAINED OPTIMIZATION[END_REF]; we obtain a better estimate for the solution; then we decrease µ by a xed factor, and restart the solver on problem (4.3) with this solution; and so on, in a continuation process [Bobin et al., 2011] detailed in Alg. 1. This algorithm is not faster than O(1/ϵ ): indeed a rst-order algorithm for the sub-problem (4.3) has optimal worst-case iteration complexity O( L µ /ϵ ) [START_REF] Nesterov | A method of solving a convex programming problem with convergence rate O(1/k 2 )[END_REF], and L µ ∼ 1/µ ∼ 1/ϵ. We believe that this bound is tight but a detailed analysis is beyond the scope of this manuscript.

Algorithm 1: LBFGS algorithm with continuation Require: ϵ > 0 the desired precision, β (0 < β < 1) the rate of decay of the smoothing parameter µ, and γ > 0 be a constant. Finally, let LBFGS: (E µ , w (0) , ϵ ) → w be an oracle which when warm-started with an initial guess w (0) , returns an ϵ-optimal solution (i.e

E µ (w) -E * µ < ϵ) for problem (4.3). 1: Initialize 0 < µ (0) (= 10 -2 , e.g), w (0) ∈ R p , and k = 0. 2: while γ µ (k ) ≥ ϵ do 3: w (k +1) ← LBFGS(E µ (k ) , w (k ) , γ µ (k ) ) 4: µ (k +1) ← β µ (k ) 5: k ← k + 1 6: end while

Experiments on fMRI datasets

We now detail experiments done on publicly available data. All experiments were run full-brain without spatial smoothing.

Visual recognition. Our rst benchmark dataset is a popular block-design fMRI dataset from a study on face and object representation in human ventral temporal cortex [Haxby et al., 2001]. It consists of 6 subjects with 12 runs per subject. In each run, the subjects passively viewed images of eight object categories, grouped in 24-second blocks separated by intermittent rest periods. This experiment is a classi cation task: predicting the object category. We use a two-class prediction target: y encodes faces versus houses. The design matrix X is made of time-series from the full-brain mask of p = 23 707 voxels over n = 216 TRs, of a single subject (subj1).

Mixed Gambles. Our second benchmark dataset is a study in which subjects were presented with mixed (gain/loss) gambles, and decided whether they would accept each gamble [START_REF] Tom | The Neural Basis of Loss Aversion in Decision-Making Under Risk[END_REF]. No outcomes of these gambles were presented during scanning, but after the scan three gambles were elected at random and played for real money. The prediction task here is to predict the magnitude of the gain and thus a regression on a continuous variable [Jimura and Poldrack, 2012]. The data The are pulled from 16 subjects with 48 3D scans each, making up for a total of n = 768 samples with approximately p = 3.3 × 10 4 voxels.

We study the convergence of the algorithms for parameters close to the optimal parameters set by 10-fold cross-validation to maximize prediction accuracy.

Results: convergence times

Here, we present benchmark results for our experiments. Figure 4.2 gives results for the logistic regression run on the visual recognition dataset: convergence plots of energy as a function of time show that all methods are asymptotically decreasing. The left part of Fig. 4.2 shows the time required to give a convergence threshold, de ned as a given excess energy compared to the lowest energy achieved by all methods, for di erent choices of regularization parameters. Similarly, the right part of Fig. 4.2 shows convergence times for squared loss on both datasets. For these gures, each solver was run for a maximum of 1 hour per problem. Solvers that do not appear on a plot did not converge for the corresponding tolerance and time budget.

For logistic loss, the most serious contender is algorithm 1, LBFGS applied on a smooth surrogate, followed by ADMM, however ADMM performance varies markedly depending on the choice of ν (more on this in chapter 7). For the squared loss FISTA and algorithm 1 are the best performers, with FISTA achieving a clear lead for the larger mixed-gambles dataset. Note that in the case of strong regularization the problem is better conditioned, and rst-order methods such as the primal-dual approach can perform well.

Bibliography

A. Beck [Dohmatob et al., 2015], we developed some heuristics for speeding up the overall optimization process: (a) Earlystopping, whereby one halts the optimization process when the test score (performance on left-out data) for the internal cross-validation for modelselection stops improving, and (b) univariate feature-screening, whereby irrelevant (non-predictive) voxels are detected and eliminated before the optimization problem is entered, thus reducing the size of the problem. Empirical results with GraphNet on real MRI (Magnetic Resonance Imaging) datasets indicated that these heuristics are a winning strategy, as they add speed without sacri cing the quality of the predictions / classi cations.

One notes that in the case of GraphNet, the penalty term of problem (3.1), the ∇w 2 2 sub-term is smooth (i.e di erentiable) with Lipschitz gradient, whilst the 1 term though nonsmooth, is proximable by means of the soft-thesholding operator [Daubechies et al., 2004]. Thus problem (3.1) is amenable to the FISTA (Fast Iterative Shrinkage-Thresholding Algorithm) [Beck and Teboulle, 2009], with a provable O(1/ √ ϵ ) convergence rate. Our implementation of FISTA uses technical recommendations (line-searching, parametrization, etc.) which were provided in [Dohmatob et al., 2014], in the context of TV-L1 [Baldassarre et al., 2012, Gramfort et al., 2013]. The model parameters α and ρ in (3.1) are set by internal cross-validation. [Hebiri andvan de Geer, 2011, Grosenick et al., 2013] problem (3.1) on di erent datasets. This gure shows spatial maps of X T j y, thresholded so that only voxels j with (from left to rightmost column)

|X T j y| ≥ p 10% (|X T y|), |X T j y| ≥ p 20% (|X T y|), |X T j y| ≥ p 50% (|X T y|), and |X T j y| ≥ p 100% (|X T y|) (full- brain
) respectively, survive. The green contours enclose the elite voxels which are selected by the screening procedure at the respective threshold levels. (a): Mixed Gambles dataset [Jimura and Poldrack, 2012].

Weights maps obtained for the GraphNet model (3.1) with these di erent screening-percentiles are shown in Figure 5.4. (c): OASIS dataset [START_REF] Daniel S Marcus | Open Access Series of Imaging Studies (OASIS): cross-sectional MRI data in young, middle aged, nondemented, and demented older adults[END_REF] with VBM. See Figure 5.2 for weights maps and age predictions obtained using these di erent screeningpercentiles.

Methods

Univariate feature-screening

In machine-learning, feature-screening aims at detecting and eliminating irrelevant (non-predictive) features thus reducing the size of the underlying optimization problem (here problem (3.1)). The general idea is to compute for each value of the regularization parameter, a relevance measure for each feature, which is then compared with a threshold (produced by the screening procedure itself). Features which fall short of this threshold are detected as irrelevant and eliminated. For the Lasso and similar models (including Group Lasso), exact screening techniques (i.e, techniques which don't mistakenly discard active predictive features) include those developed in [START_REF] Ghaoui | Safe Feature Elimination in Sparse Supervised Learning[END_REF][START_REF] Lee | Exact post model selection inference for marginal screening[END_REF][START_REF] Liu | Safe Screening With Variational Inequalities and Its Applicaiton to LASSO[END_REF], Wang et al., 2015]. Inexact screening techniques (e.g [START_REF] Tibshirani | Strong rules for discarding predictors in lasso-type problems[END_REF]) have also been proposed in the literature.

Our proposed heuristic screening technique is inspired by the Marginal screening technique developed in Algorithm 1 of [START_REF] Lee | Exact post model selection inference for marginal screening[END_REF], and operates as follows. The data (X, y) are standardized so that y has unit variance and zero mean, likewise each row of the design matrix X. To ensure obtention of a smooth mask, a Gaussian-smoothed version of X is used in the screening procedure (but not in the actual model t). For each voxel j (voxels are the features here) the absolute dot-product |X T j y| of y with the jth column of X is computed. For a given screening-percentile sp ∈ [0, 100] , the spth percentile value of the vector |X T y| := (|X T 1 y|, …, |X T p y|), denoted p sp (|X T y|), is computed. The case sp = 100 corresponds to fullbrain analysis with no screening. sp = 25 means we keep the quarter of the brain made of voxels with the highest |X T j y| values, and so on. A brain-mask is then formed containing only those voxels j for which |X T j y| ≥ p sp (|X T y|). Next, this brain-mask is morphologically eroded and then dilated, to obtain a more structured mask. Figure 5.1 shows results of applying this screening heuristic to various datasets, prior to model tting.

Early-stopping

Optimization is a means to an end and not an end on its own. The only incentive for optimizing a model is to improve its generalization power: performance on unseen data. If this performs stops improving during training (statistical convergence), as measured on a left-out subset of data , then we may as well abrupt the optimization algorithm. We implement this principle heuristically as follows. In each train sub-sample of the internal crossvalidation loop for setting the parameters of the GraphNet model (3.1), a pass is done on the 2-dimensional parameter grid (see Fig. 8.4) and each parameter pair (α, ρ) is scored according to its prediction / classi cation performance. For a xed parameter pair (α, ρ), an instance of problem (3.1) is solved iteratively using FISTA Beck and Teboulle [2009]. At each iteration, the prediction / classi cation performance of the current (not yet optimal) solution ŵk in (3.1) is computed. If in a time-window of 5 iterations this score has not increased above an a priori xed threshold, called the early-stopping tolerance (es tol), then the optimization process is halted for the current model parameter pair (α, ρ) under inspection. This heuristic is motivated by the intuition that, for a particular problem, sub-optimal solutions ŵk can give the same score as an optimal solution ŵ (i.e "statistical convergence" happens before numerical convergence). By default we set this early-stopping tolerance to -10 -4 for classi cation and -10 -2 for regression problems. A value of +∞ (in fact, any value above 10, say) corresponds to no early-stopping at all (i.e, solve problem (3.1) until numerical convergence).

Experiments

We experimented our early-stopping and (separately) feature-screening heuristics on di erent MRI datasets. All experiments were run using a single core of a laptop.

Regression. The OASIS dataset [START_REF] Daniel S Marcus | Open Access Series of Imaging Studies (OASIS): cross-sectional MRI data in young, middle aged, nondemented, and demented older adults[END_REF] consists of a crosssectional collection of 416 subjects aged 18 to 96. For each subject, 3 or 4 individual T1-weighted MRI scans obtained in single scan sessions are included. A natural regression problem for this dataset is to predict the age of a subject from their anatomical data. To this end, we segmented the gray-matter from the anatomical data of each subject (obtained from the T1 images), and used the gray-matter maps as features for predicting age. We split the 416 subjects into two equally-sized and age-balanced groups: a train set and a validation set. The GraphNet model [Hebiri andvan de Geer, 2011, Grosenick et al., 2013] was tted on the train set, with parameters (α and ρ in (3.1)) set internally via 8-fold cross-validation. The results for this experiment are shown in Figure 5.2.

Classification. The visual recognition dataset [Haxby et al., 2001] is a popular block-design fMRI dataset from a study on face and object representation in human ventral temporal cortex. It consists of 6 subjects with 12 sessions / runs per subject. In each run, the subjects passively viewed images of eight object categories, grouped in 24-second blocks separated by intermittent rest periods. This experiment is a classi cation task: predicting the object category y. We use a One-versus-Rest (OvR) strategy. The design matrix X is made of time-series from the full-brain mask of p = 23 707 voxels over n = 216 TRs, of a single subject (subj1). We divided the 12 runs into 6 runs for training and 6 other runs for validation. Leave-one-session-out 1 1 One session is held out and the other S -1 sessions are used to train a model which is validated on the left-out session. This is repeated for all the sessions, yield an estimate -with error bars-on the generalization error of the model.

cross-validation was used for selecting the model parameters (α, ρ). The results are depicted in Figure 5.4.

Results

We now summarize and comment the results of the experiments (refer to section 4.1.2). Figure 5.2 shows the e ects of early-stopping heuristic and feature-screening heuristic on age prediction scores on the OASIS dataset [START_REF] Daniel S Marcus | Open Access Series of Imaging Studies (OASIS): cross-sectional MRI data in young, middle aged, nondemented, and demented older adults[END_REF] (416 subjects). We see that in the internal cross-validation, stopping the optimization procedure for xed (α, ρ) pair of regularization parameters, when test score increases by -10 -2 or more is a good heuristic, and does just as good as running the optimization until numerical convergence. Also (and independently), one gets similar prediction scores using as little as a fth of the brain volume (sp = 20), compared to using the fullbrain (sp = 100). Figure 5.4 reports similar results for classi cation on the visual recognition dataset [Haxby et al., 2001]. Overall, we see from Figures 5.4 and 5.2 that we can achieve up to 10-fold speedup with the proposed heuristics, with very little loss in accuracy. Also, we see that contiguous groups of bars are roughly at at the top, with a sligh increase from lower to high screening-percentile values. The case "chair vs scramped" is an exception, where a slightly reverse tendency if observed. A possible explanation is that 20th percentile feature-screening already selects the right voxels (quasi-exact support recovery), and so including more voxels in the model can only hurt its performance.

Conclusion

We have presented heuristics that provide speedups for optimizing Graph-Net [Hebiri andvan de Geer, 2011, Grosenick et al., 2013] in the di cult context of brain data. These heuristics are a winning strategy as they add speed without sacri cing the quality of the predictions / classi cations. In practice, we do a 20univariate feature-screening by default, which ensures a 5-fold speedup over full-brain analysis, and independently of an approximately 2-fold speedup obtained by the early-stopping heuristic, leading to an overall 10-fold speedup. Our results have been veri ed empirically on di erent MRI datasets. Our heuristics should be applicable to other hard-tooptimize models like TV-L1 [Baldassarre et al., 2012, Gramfort et al., 2013]. The result of these numerous contributions on optimizing the SpaceNet model (3.1) have been implemented as part of the Nilearn package [Abraham et al., 2014].

On the equivalence of TV-L1 and iteratively-reweighted GraphNet We present a result that shows that TV-L1 regularized regression problems (3.1) can also be solved through an iteratively reweighted GraphNet problem. Among other things, this provides a long-awaited statistical interpretation of the TV-L1 penalized models (3.1). The method dubbed iGraph-Net, solves TV-L1 penalized model by considering modi ed GraphNet subproblems corresponding to the minimization of the energy E γ GraphNet (w) de ned in (6.5). These sub-problems are very well-conditioned and are quadratically easier to solve than TV-L1 itself. The limit of this sub-problems is solve the exact TV-L1 penalized problem.

This work follows the spirit of [START_REF] Emmanuel | Enhancing Sparsity by Reweighted l1 Minimization[END_REF] which proposed an enhanced Lasso problem built iteratively from surrogate Ridge regression problems with inhomogeneous feature penalty parameters. However, unlike [START_REF] Emmanuel | Enhancing Sparsity by Reweighted l1 Minimization[END_REF], we leave the Lasso part of the TV-L1 penalty (3.2) untouched and instead derive a surrogate on the TV part, which turns out to be a GraphNet problem with inhomogeneous penalty parameters. See Figure 6.1. Pending gures comparing (maps, scores, and runtime) Graph-Net, iGraphNet, and the baseline TV-L1 implementation via double-FISTA implementations [Dohmatob et al., 2014, Varoquaux et al., 2015].

Derivation

Invoking the following well-known elementary result 1 1 To prove it, one simply uses the fact that

w 2 u 2 + 1 -2wu = (wu - 1) 2 ≥ 0, with equality i wu = 1. ∀u, w > 0, u ≤ wu 2 + w -1 2 , with equality i w = u -1 , (6.1)
we can rewrite the TV semi-norm as follows,

w TV := j ∈[[p]], (∇w) j 2 >0 (∇w) j 2 ≤ 1 2 j ∈[[p]], (∇w) j 2 >0 γ j (∇w) j 2 2 + γ -1 j , ∀γ ∈ R p ++ , (6.2)
with equality i

γ j = (∇w) j -1 2 , ∀j ∈ [[p]] s.t (∇w) j 2 > 0. (6.3)
Thus,

w TV = min γ ∈R p ++ 1 2 j ∈[[p]], (∇w) j 2 >0 γ j (∇w) j 2 2 + γ -1 j , (6.4)
with the optimal scaling vector γ ∈ R p ++ given by (6.3). Whence, the minimizers of the TV-L1 energy E TV-L1 (w) := (y, Xw) + α P TV-L1 (w) coincide with the minimizers of the rescaled GraphNet energy

E γ GraphNet (w) = (y, Xw) + α ρ w 1 + 1 2 α (1 -ρ) j ∈[[p]], (∇w ) j 2 >0 γ j (∇w) j 2 2 +γ -1 j ,
(6.5) where the minimization is done both over regression coe cients w and the scaling parameters γ 1 , . . . , γ p > 0.

Algorithm 2: iGraphNet: iteratively-reweighted GraphNet solver for the TV-L1 model Require: Values for the model-tuning parameters λ > 0, and 0 ≤ ρ ≤ 1; initial brain-map w (0) ∈ R p (e.g, the zero vector); tolerance threshold ϵ > 0 (say 10 -5 ); maximum number of outer iterations K.

Ensure: An optimal vector ŵTV of regressor coe cients (an approximation of) for the TV-L1 model.

1: Initialize: k ← 0; µ ← 10 -4 2: while w (k +1) -w (k ) ∞ ≥ ϵ do 3:
Recompute scaling:

γ (k ) j ← ( (∇w (k ) ) j 2 2 + µ 2 ) -1 2
, for every voxel j 4:

Recompute coe cients: w (k +1) ← arg min w∈R p E γ (k ) GraphNet (w), with energy tolerance ∼ µ. The solver for this sub-problem is warmstarted with w = w (k ) .

5:

Goto next iteration: k ← k + 1 6: end while As a function of the regressor coe cients w, the energy in (6.5) corresponds to a modi ed GraphNet model in which per-voxel penalty parameters α (1 -ρ)γ j given by (6.3) replace the constant α (1 -ρ) factor in the pure GraphNet model (3.1), or equivalently the ∇ is pre-whitened by the diagonal matrix

Γ := diag( √ γ 1 , . . . , √ γ p ).
This energy is optimized by an alternating scheme cyclically switching between optimizing w.r.t the regressor coe cients w and then w.r.t then rescaling parameters γ 1 , . . . , γ p (in closed form, via formula (6.3)). The algorithm so-obtained (detailed in section 6.2) alternates between minimization over the scaling parameters γ and minimization over the coe cients w.

The algorithm: iGraphNet

We now present iGraphNet, an iteratively-reweighted scheme for solving the TV-L1 model, based on modi ed GraphNet (3.1) sub-problems corresponding to the minimization of the energy E γ GraphNet (w) de ned in (6.5).

These sub-problems are very well-conditioned and are quadratically easier to solve than TV-L1 itself, and can be solved by a fast rst-order method like FISTA or LARS. The algorithm is presented in Alg. 2.

Overall, for a tolerance ϵ > 0, Alg. 2 converges in O(1/ϵ ) basic iterations (i.e counting all the iterations run in a rst-order method for solving the GraphNet sub-problem), though its observed runtime is in the order of about K times the time taken by a run of a solver for the GraphNet subproblem. Practical details (like handling a brain mask, automatic model parameter selection via cross-validation and bagging, early-stopping, etc.) that go in the implementation of the optimization algorithms like the one just presented can be found in [Dohmatob et al., 2014].

Generalization to other complex non-smooth models. Similarly, one can show that Sparse-Variation [Eickenberg et al., 2015] can be solved via an IRLS (iteratively-reweighted Least Squares) scheme, where the weights are computed via (6.3), with the ∇ operator replaced with an identity-augmented version. Indeed, thanks to the inequality (6.1), it turns out that most complex rich non-smooth p -norm-based models are just iteratively-reweighted versions of much simpler counterparts like Ordinary Least Squares, Lasso, ElasticNet, GraphNet, etc.

Experimental results

Preliminary experimental results are shown in Fig. 6.1. We run our iGraph-Net procedure (Alg. 2) on data for the Face vs House condition of the visual recognition dataset [Haxby et al., 2001]. Model coe cients and accuracies on held-out data are shown. We monitor the evolution of the model as a function of the number of iGraphNet iterations k = 0, 1, 2, . . .. We see that as more and more iterations of iGraphNet are run, the coe cients become more and more spatially denoised and localized (and therefore more intepretable), without deterioration of prediction accuracy. A result on the rate of convergence of the ADMM algorithm

Introduction

T ADMM [START_REF] Glowinski | Sur l'approximation, par élémentsnis d'ordre un, et la résolution, par pénalisation-dualité d'une classe de problemes de Dirichlet non linéaires[END_REF][START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via nite element approximation[END_REF][START_REF] Eckstein | On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF]] is an operator-splitting optimization method which is easy to implement and well-adapted for large-scale optimization problems [START_REF] Boyd | Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers[END_REF]. ADMM can provide a distinctive advantage over proximal gradient methods such as FISTA [Beck and Teboulle, 2009] when there is no closed-form expression for the proximal operator. Indeed, ADMM can avoid this di culty by introducing a "split" variable, for which the proximal operator results in updates computable in closedform. This is typically the case in analysis sparsity regularization, that impose sparsity on a transformation of the optimization variable. However, the theory of the convergence rate of ADMM is not complete [START_REF] Boyd | Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers[END_REF].

In our ICASSP 2016 paper [Dohmatob et al., 2015], we studied the convergence of the ADMM (Alternating Direction Method of Multipliers) algorithm on a broad range of penalized regression problems including the Lasso, Group-Lasso and Graph-Lasso,(isotropic) TV-L1, Sparse Variation, and others, that can be written in the form minimize

(w,z) ∈R p ×R q 1 2
Xwy 2 + λΩ(z) subject to Kwz = 0, (7.1) where X ∈ R n×p is the design matrix; y ∈ R n is a vector of measurements or classi cation targets; K ∈ R q×p is linear operator; λ > 0 is the regularization parameter; and Ω : R p → (-∞, +∞] is the penalty, which is assumed to be a closed proper convex function. This is an instance of the SpaceNet model (3.1) presented in chapter 3. In signal processing literature, (7.1) is an example of what is referred to as a synthesis problem: the penalty Ω is imposed not directly on the image, but on a the output of a dictionary, z = Kw. K is referred to the analysis operator. The case K = I corresponds to the synthesis setting.

The ADMM algorithms

Consider the ADMM algorithm [START_REF] Glowinski | Sur l'approximation, par élémentsnis d'ordre un, et la résolution, par pénalisation-dualité d'une classe de problemes de Dirichlet non linéaires[END_REF][START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via nite element approximation[END_REF][START_REF] Eckstein | On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF][START_REF] Boyd | Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers[END_REF] applied to problem (7.1). Let µ ∈ R q be the dual variable and ν > 0 be the penalty parameter on the splitting residual. The augmented Lagrangian is:

L ν (w, z, µ) = 1 2 Xw -y 2 + λΩ(z) + µ T (Kw -z) + 1 2 ν Kw -z 2 .
Further, introducing the scaled dual variable u := ν -1 µ, which we will use instead of µ from here on, the ADMM iterates for problem (7.1) are given by the following equations:

w (n+1) ← arg min w L ν (w, z (n) , u (n) ) = (ν K T K + X T X) -1 (ν K T (z (n) -u (n) ) + X T y) z (n+1) ← arg min z L ν (w (n+1) , z, u (n) ) = prox (α /ν )Ω (Kw (n+1) + u (n) ) u (n+1) ← u (n) + Kw (n+1) -z (n+1) . (7.2)
Assumptions. We will assume that the matrix sum ν

K T K + X T X is in- vertible.
This assumption is equivalent to ker K T K ∩ ker X T X = {0} (see e.g [Piziak et al., 1999, Theorem 1]), which is reasonable in the context of regularization. Indeed, the idea behind this assumption is that, in highdimensional problems (n p), X typically has a large kernel, and so one would naturally choose K to act on it.

Examples

Problem (7.1) covers a broad spectrum of problems encountered in pattern recognition and image processing. Here are a few:

Classical examples. We have Ω = 1 2 . 2 for Ridge regression; Ω = . 1 : z → j ∈[p] |z j | for Lasso and Fused-Lasso [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF]. For all but the last of these examples, we have K = I. For Group-Lasso, we have

K = I, Ω = the mixed-norm 2,1 = . 2,1 : z → j ∈[[d ]] z j :j+c-1
, where there are d ≥ 1 blocks z j :j+c-1 := (z j , z j+1 , …, z j+c-1 ) each of size c ≥ 1.

Isotropic TV-L1 and Sparse Variation. The di erent extensions of the TV penalty presented in chapter 3 can be posed in the form of the problem above. For example, Sparse Variation [Eickenberg et al., 2015] corresponds to taking K = [ρI, (1 -ρ)∇] T ∈ R 4p×p , where ∇ is the discrete (refer to chapter 3) spatial gradient operator and ρ ∈ [0, 1] is a mixing parameter. For TV-L1 [Baldassarre et al., 2012, Gramfort et al., 2013], the penalty is given by Ω(z) = j ∈[[p]] |z j,1 | + j ∈[p] z j,2:4 (i.e an 1 norm on the rst p coordinates of z and an 2,1 mixed-norm on the last 3p coordinates). In particular, the case ρ = 1 corresponds to the usual 1 norm, while ρ = 0 corresponds to the isotropic TV semi-norm.

In Sparse Variation [Eickenberg et al., 2015], the penalty is modi ed to simply be an 2,1 mixed-norm on d = p blocks of size c = 4 each, i.e Ω(z) = j ∈[p] z j,1:4 . TV-L1 and Sparse Variation combine sparsity (due to the the 1 -norm) and structure (due to the isotropic TV term) to extract local concentrations of spatially correlated features from the data.

Our contributions 7.2.1 Preliminaries

In the spirit of [START_REF] Ghadimi | Optimal parameter selection for the alternating direction method of multipliers (ADMM): Quadratic problems[END_REF], let us start with a simple lemma (proof omitted) which rewrites the ADMM iterates (7.2) as a Picard xed-point process in terms of the (z, u) pair of variables.

Lemma 1. De ne the following objects:

G ν := K(K T K + ν -1 X T X) -1 K T , A ν := [G ν I -G ν ], b ν := ν -1 K(K T K + ν -1 X T X) -1 X T y, Ãν := A ν (.) + b ν , Λ ν := prox (α /ν )φ • Ãν , (I -prox (α /ν )φ ) • Ãν .
Then the z and u updates in the ADMM iterates (7.2) can be jointly written as a Picard xed-point iteration for the operator Λ ν , i.e (z (n+1) , u (n+1) ) ← Λ ν (z (n) , u (n) ).

(7.3)

In the special case where prox (α /ν )φ is a linear transformation -as in Ridge regression or the nonnegative Lasso, for example-the operator Λ ν is linear so that the xed-point iteration (7.3) is a linear dynamical system. Moreover, in such cases one can derive closed-form formulae for the spectral radius r (Λ ν ) of Λ ν as function of ν , and thus recover the results of [START_REF] Ghadimi | Optimal parameter selection for the alternating direction method of multipliers (ADMM): Quadratic problems[END_REF] and [START_REF] Boley | Local Linear Convergence of the Alternating Direction Method of Multipliers on Quadratic or Linear Programs[END_REF]. In the latter simple situations, a strategy for speeding up the ADMM algorithm is then to choose the parameter ν so that the spectral radius of the linear part of the then a ne transformation Λ ν is minimized. The following Corollary is immediate, whose proof is obtainable via the Spectral Mapping Theorem.

Corollary 1. Let G ν , A ν , Ãν , and Λ ν be de ned as in Lemma 1. Then the following hold:

(a) max( G ν , I -G ν ) ≤ 1, ν min * (A ν ) ≥ 1/
√ 2, and A ν ≤ 1 with equality in the last inequality i at least one of G ν and I -G ν is singular.

(b) Λ ν is A ν -Lipschitz. That is, ∀(x 1 , x 2 ) ∈ R q+q × R q+q , Λ ν (x 1 ) -Λ ν (x 2 ) ≤ A ν x 1 -x 2 .
(7.4)

In particular, if A ν < 1, then Λ ν is a contraction and the ADMM iterates (7.2) converge globally Q-linearly to a solution of (7.1). Moreover, this solution is unique.

According to Corollary 1, Λ ν is an A ν -contraction in case A ν < 1, and so we have global Q-linear convergence of the ADMM iterates ( 7.2) at the rate A ν . This particular case is analogous to the results obtained in [START_REF] Nishihara | A General Analysis of the Convergence of ADMM[END_REF] when the loss function or the penalty is strongly

convex. But what if G ν = I -G ν = A ν = 1 ?
Can we still have Q-linear convergence, -at least locally ? These questions are answered in the sequel.

Behavior of ADMM around fixed-points

10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 10 Figure 7.1: Rate of convergence r (Λ ν (x * )) as a function of ν for a Lasso problem with column-rank de cient design matrix X. Taking ν too small leads to badly conditioned problem (as I + (1/ν )X T X is then almost singular), and thus a slow rate of convergence (near 1). On the other hand, the gure suggests that taking ν "too large" is also detrimental. Most remarkable, one notices that the basin of "good" ν values is rather tight, and so care must be taken in choosing the ν parameter.

Henceforth, we consider problem (7.1) in situations where the penalty φ is an 2,1 mixed-norm. Note that the 1 -norm is a special case of the 2,1 mixednorm with c = 1 feature per block, and corresponds to the anisotropic case.

The results presented in Theorem ( 1) carry over e ortlessly to the case where the φ is the concatenation of 2,1 norms, for example as in the the TV-L1 semi-norm. The following theorem -inspired by a careful synthesis of the arguments in [START_REF] Holmes | Smoothness of certain metric projections of Hilbert space[END_REF] and [Bayram and Selesnick, 2010]-is our main result.

Our main results are summarized in Theorem 1 of the aforementioned paper, which we now state.

Theorem 1. Consider the ADMM algorithm (7.2) on problem (7.1), where Ω is an 2,1 mixed-norm on d ≥ 1 blocks each of size c ≥ 1, for a total of q = d × c features. Let the operators A, Ã, and Λ be de ned as de ned above, with the ν subscript dropped for ease of notation. Let For x = (z, u) ∈ R q+q , let Λ 1 (x) ∈ R q denote the rst q coordinates of Λ(x), i.e its z-part. De ne

• supp(z) := {j ∈ [[d]] | z j :j+c-1 0}; • A 1 (z) := {z ∈ R q | supp(z ) = supp(z)}, and A(x) := A 1 (z) ⊕ R q ; • x := ( xj ) j ∈[[d ]] := Ãx, κ := α /ν , ϵ (x) := min j ∈[[d ]] | xj -κ | ≥ 0.
Then the following hold:

(a) A ractivity of supports. For all x ∈ R q+q , we have

Λ( B2q (x, ϵ (x)/ A )) ⊆ B2q (Λ(x), ϵ (x)) ∩ A (Λ(x)).
In particular, if x * is a xed-point of the operator Γ, then (7.5) where F x := [D x I -D x ] T and D x ∈ R q×q is a block-diagonal matrix with block D x,j ∈ R c×c given by

Λ( B2q (x * , ϵ (x * )/ A )) ⊆ B2q (x * , ϵ (x * )) ∩ A (x * ). (b) Fréchet-di erentiability. If x ∈ R q+q with ϵ (x) > 0, then Λ is Fréchet- di erentiable at x with derivative Λ (x) = F x A ∈ R 2q×2q ,
D x,j =        I -κ xj P xj ⊥ , if j ∈ supp(Λ 1 (x)), 0,
otherwise.

(7.6)

In particular, when c = 1, each D x , j reduces to a bit ∈ {0, 1} which indicates whether the jth feature is active, and D x reduces to a diagonal projector matrix with only 0s and 1s.

(c) Let x * = (z * , u * ) ∈ R q+q be any xed-point of Γ.

(1) Finite-time identi cation of active set. If the closed ball B2q (x * , ϵ (x * )/ A ) contains any point of the sequence of iterates x (n) , then the active set A (x * ) is identi ed after nitely many iterations, i.e

∃N x * ≥ 0 s.t x (n) ∈ A(x * )∀n ≥ N x * . (7.7)
In particular, (7.7) holds if x (n) converges to x * .

(2) Local Q-linear convergence. If ϵ (x * ) > 0 and r (Λ (x * )) < 1, then the iterates x (n) converge locally Q-linearly to x * at the rate r (Λ (x * )).

(3) Optimal rates in the anisotropic case. If c = 1 (as in anisotropic TV deconvolution) and ν is large, then the optimal rate of convergence rate is the cosine of the Friedrichs angle between Im K and Im D x * A 1 (z * ).

If in addition K = I (as in synthesis inverse problems like the Lasso, sparse Spike-deconvolution, etc.), then the whole algorithm converges in a nite number of iterations.

Proof. See our ICASSP paper [Dohmatob et al., 2015].

Relation to prior work

Recently, there have been a number of results on the local linear convergence of ADMM on particular classes of problems. Below, we outline the corresponding major works.

Ridge, QP, and nonnegative Lasso

On problems like Ridge regression, quadratic programming (QP), and nonnegative Lasso, [START_REF] Ghadimi | Optimal parameter selection for the alternating direction method of multipliers (ADMM): Quadratic problems[END_REF] demonstrated local linear convergence of ADMM under certain rank conditions which are equivalent to requiring that the p.s.d matrix G ν (de ned in (7.3)) be invertible. The same paper prescribed explicit formulae for optimally selecting the tuning parameter ν for ADMM on these problems. We note that these results can be recovered from our Lemma 1 and Corollary 1 as they correspond to the case where prox (α /ν )φ is a linear operator. Using similar spectral arguments, [START_REF] Boley | Local Linear Convergence of the Alternating Direction Method of Multipliers on Quadratic or Linear Programs[END_REF] demonstrated similar local convergence results for quadratic and linear QP problems.

works, at least on models like isotropic TV-L1, Sparse Variation, etc., since these models are not PSS. Indeed, the convergence rates in [START_REF] Liang | Activity Identication and Local Linear Convergence of Douglas-Rachford/ADMM under Partial Smoothness[END_REF][START_REF] Liang | Activity Identi cation and Local Linear Convergence of Inertial Forward-Backward Splitting[END_REF] do not extend from anisotropic to isotropic TV, for example. Success in the former case is due to the fact that the anisotropic TV seminorm is polyhedral and therefore is of class PSS at each point. By contrast, our framework can handle isotropic TV and similar "entangled" penalty types like isotropic TV-L1, Sparse Variation, etc., but su ers complementary limitations; for example, we were unable to generalize it beyond the squared-loss setting and we can only handle penalties which are a composition of a 2,1 mixed-norm (or a concatenation of such) and a linear operator. The recent work [START_REF] Vaiter | The degrees of freedom of partly smooth regularizers[END_REF] on counting the degrees of freedom of general partly-smooth penalties is worth mentioning and may contain some key ideas to help bridge the "isotropicity gap " in the methods developed in [START_REF] Liang | Activity Identication and Local Linear Convergence of Douglas-Rachford/ADMM under Partial Smoothness[END_REF][START_REF] Liang | Activity Identi cation and Local Linear Convergence of Inertial Forward-Backward Splitting[END_REF], concerning rates of convergence.

Lastly, the convergence rates in [START_REF] Liang | Activity Identication and Local Linear Convergence of Douglas-Rachford/ADMM under Partial Smoothness[END_REF][START_REF] Liang | Activity Identi cation and Local Linear Convergence of Inertial Forward-Backward Splitting[END_REF] are tight and given in terms of Friedrichs angles [Bauschke et al., 2014], whilst our rates are given in terms of spectral radii, and will be suboptimal in certain cases. An exception are the anisotropic cases, where we proved in part (c3) of Theorem 1 that we recover the optimal rates obtained in [START_REF] Liang | Activity Identication and Local Linear Convergence of Douglas-Rachford/ADMM under Partial Smoothness[END_REF][START_REF] Liang | Activity Identi cation and Local Linear Convergence of Inertial Forward-Backward Splitting[END_REF] in terms of Friedrichs angles. Moreover, for the Lasso, we showed that the whole algorithm converges after only nitely many iterations.

Numerical experiments and results

Here, we present results for a variety of experiments. Each experiment is an instance of problem (7.1) with an appropriate choice of the linear operators X, K, and the penalty function φ which can be the 1 -norm the 2,1 mixednorm, or a mixture of the two (as in TV-L1).

Se ing. We use a grid of 20 values of ν , evenly spaced in log-space from 10 -3 to 10 6 . For each problem model (see below), the iteration process (7.3) is started with x (0) = 0 ∈ R q×q , and iterated N = 1500 times. The nal point x (N ) is approximately a xed-point x ( * ) of the operator Λ ν . Now, the iteration process is run again (starting with the same initial x (0) ) and the distance x (k )x (N ) is recorded on each iteration k, producing a curve. This procedure is run for each value of ν from the aforementioned grid. Except otherwise stated, the n rows of design matrix X where drawn from a p-dimensional standard Gaussian. The measurements variable y is then computed as y = Xw 0 + noise, where w 0 is the true signal.

Simple models. As discussed in section 7.3, the local Q-linear convergence of ADMM on a variety of particular problems has been studied in the literature (for example [START_REF] Ghadimi | Optimal parameter selection for the alternating direction method of multipliers (ADMM): Quadratic problems[END_REF][START_REF] Nishihara | A General Analysis of the Convergence of ADMM[END_REF][START_REF] Liang | Activity Identication and Local Linear Convergence of Douglas-Rachford/ADMM under Partial Smoothness[END_REF][START_REF] Pedregosa-Izquierdo | Feature extraction and supervised learning on fMRI: from practice to theory[END_REF]). We validated empirically our linear convergence results (Theorem 1) by reproducing experiments from [START_REF] Liang | Activity Identication and Local Linear Convergence of Douglas-Rachford/ADMM under Partial Smoothness[END_REF][START_REF] Liang | Activity Identi cation and Local Linear Convergence of Inertial Forward-Backward Splitting[END_REF]. For each of these experiments the regularization parameter α was set to 1. Viz, (a) Lasso: Here the problem is an instance of (7.1) with K = I and φ = • 1 ; n = 32, q = p = 128, and w 0 is 8-sparse.

Concluding remarks

We have derived a xed-point iteration which is equivalent to the ADMM iterates for a broad class of penalized regression problems (7.1). Exploiting the formulation so obtained, we have established detailed qualitative properties of the algorithm around solution points (Theorem 1). Most importantly, under mild conditions, local Q-linear convergence is guaranteed and we have provided an explicit formula for this rate of convergence. Finally, Theorem 1 -implicitly-opens the possibility of speeding up the ADMM algorithm on problem (7.1) by selecting the tuning parameter ν so as to minimize the spectral radius (an inverted mexican-hat-shaped curve, as ν varies from 0 to +∞) of the Jacobian matrix T ν (x * ). In this chapter, we turn our attention to another key problem in neuroimaging: registration of functional brain data from di erent subjects. After a concise review of the existing literature, we present a contribution of ours, namely direct registration of fMRI data without using the subject's anatomy as a proxy.

Introduction

Registering brain images from di erent subjects in a common space (for example, the MNI space [START_REF] Collins | Automatic 3D intersubject registration of MR volumetric data in standardized Talairach space[END_REF][START_REF] Mazziotta | A probabilistic atlas of the human brain: theory and rationale for its development. The International Consortium for Brain Mapping (ICBM)[END_REF]), is an essential step in any multi-subject analysis pipeline [Friston et al., 1995]. Indeed, a voxel-to-voxel correspondence across subjects is needed for grouplevel statistics on brain maps to make sense. In addition, the use of a standard space opens the possibility to share results in a consistent fashion, hence the comparison of experiments and meta-analysis [START_REF] Wager | Meta-analysis of functional neuroimaging data: current and future directions[END_REF][START_REF] Gorgolewski | .org: a web-based repository for collecting and sharing unthresholded statistical maps of the human brain[END_REF]. This is especially true in fMRI (functional Magnetic Resonance Imaging) studies in which the activations might span just a few voxels in diameter.

Traditional indirect T1-based techniques for inter-subject registration of EPI data assume that the mismatch between a subject's T1 (i.e anatomical) image and associated EPI scan is only a ne, i.e. it includes only pose di erences related to slice orientation and eld of view selection. such methods exploit the fact that learning a deformation from the subject's T1 to a template is easier, due to the relatively high anatomical contrast in T1 images, than learning a deformation from the subject's EPI image to the template. Thus the EPI and the T1 are a nely aligned in a primary step called coregistration, then one applies the transformation T1 → template to the EPI images to warp them from subject to template space. Left: Single-band high-resolution EPI (SBRef) image of the same subject. Notice the large distortions along the Left-Right direction (inside the highlighted patches).

Center: Distortion-corrected single-band EPI image. Here, the distortioncorrection managed to undo most -but not all-of the distortions. Even after distortion correction, there are minor shape (nonlinear) di erences between the EPI and the T1-weighted image of the subject (Right). The same nativespace coordinates where used in all of the 3 plots.

A crucial assumption in these classical methods is that the T1 and EPI images of the same subject can be properly aligned to one another via an a ne transformation. Thus one assumes, for example, that distortion correction is good enough that the EPI image can be realigned to the T1 with a rigid or a ne transformation. However, high-resolution EPI sequences deviate from the aforementioned underlying assumptions of classical T1-based indirect inter-subject registration methods, namely the EPI sequences su er from distortions that push them nonlinearly out-of-match relative to the T1-weighted image of the same subject.

As an example, Fig. 8.1 illustrates this issue (distortions) on the HCP (Human Connectome Project) [van Essen et al., 2012] dataset, a reference dataset that contains high-quality EPI data acquired using state-of-the-art sequences, yet with severe distortions [START_REF] Wan | Reduction of geometric and intensity distortions in echo-planar imaging using a multireference scan[END_REF], Mangin et al., 2002[START_REF] Zeng | Image distortion correction in EPI: comparison of eld mapping with point spread function mapping[END_REF][START_REF] Andersson | How to correct susceptibility distortions in spin-echo echo-planar images: application to di usion tensor imaging[END_REF]. Indeed, as discussed in the literature (e.g [START_REF] Freire | What is the best similarity measure for motion correction in fMRI time series[END_REF]), EPI distortions and signal loss related to B0 inhomogeneities cannot be separated with registration based techniques, which are compensatory operations. Consequently, the set up of e cient distortion correction method in EPI-based imaging is an open question. Moreover, sophisticated anatomy-based methods like Freesurfer's recon-all cannot scale to huge data sets like the 5,000 participants of the initial release of the UKBioBank dataset [START_REF] Miller | Multimodal population brain imaging in the UK Biobank prospective epidemiological study[END_REF], due the long computation time that renders such approaches impractical.

The goal of this paper is to provide experimental evidence that the indirect T1-based inter-subject EPI registration explained above is no longer needed, if not sub-optimal in such settings. We also provide a computationally cheap pipeline based on publicly available tools, which bypasses the need for a T1-weighted image, and do direct inter-subject registration of the EPI images.

Methods

An important note on normalization

Let us begin by stressing that the normalization problem (i.e registration to a standard template) is not addressed in our work. We concentrate on inter-subject (nonlinear) registration, since our goal is to show the bene ts of using EPI images in place of anatomical images in pipelines. We also note that there is an increasing concern in the literature that in the future, normalization will be based on multi-modal atlases (tissue probability maps, functional parcellation maps, etc.) [START_REF] Bibliography | Interoperable atlases of the human brain[END_REF].

General preprocessing procedures

Motion correction During acquisitions, subjects move their heads in the scanner. This head movement induces an approximately a ne mismatch between di erent volumes acquired in the same acquisition run. Motion correction is done to remove this source of intra-subject variability. We used FSL's irt tool [Smith et al., 2004] for motion correction.

Distortion correction Due to inhomogeneities in the ambient B0 eld, the EPI images are distorted (i.e artifactually warped) along the phase-encoding directions (Left-Right / Right-Left in the case of HCP dataset [van Essen et al., 2012]). See Figure 8.1. In our experiments, distortion correction [START_REF] Wan | Reduction of geometric and intensity distortions in echo-planar imaging using a multireference scan[END_REF], Mangin et al., 2002[START_REF] Zeng | Image distortion correction in EPI: comparison of eld mapping with point spread function mapping[END_REF][START_REF] Andersson | How to correct susceptibility distortions in spin-echo echo-planar images: application to di usion tensor imaging[END_REF][START_REF] Jezzard | Correction for geometric distortion in echo planar images from B0 eld variations[END_REF], Wan et al., 1997b] was achieved using the methods described in [van Essen et al., 2012]. Both methods use FSL's topup tool [Smith et al., 2004] to estimate the deformation eld due to B0 inhomogeneities (the distortions) [START_REF] Glasser | The minimal preprocessing pipelines for the Human Connectome Project[END_REF].

Deformation model

We used ANTs' Symmetric Normalization (aka SyN ) deformation model [START_REF] Avants | Symmetric di eomorphic image registration with cross-correlation: evaluating automated labeling of elderly and neurodegenerative brain[END_REF][START_REF] Avants | A reproducible evaluation of ANTs similarity metric performance in brain image registration[END_REF], which has been shown to be a state-of-the art method for nonlinear registration [Klein et al., 2009]. As done usually, we initialize a nonlinear registration algorithm with an a ne (rigid-body) registration algorithm. The former is simply meant to estimate an alignment for the bounding boxes of the images (thus ensuring a suciently large region of overlap). Concretely, we stack a 2-level pyramidal 1 1 Pyramidal means multiple passes are made by a registration algorithm on the input images, with ner and ner resolution (aka pyramid). In this speedup technique, each pass of the pyramid is initialized with the solution of the previous pass (this is known as warmstarting).

a ne transformation model (as initialization) with a 3-level pyramidal SyN deformation model. Mattes mutual information [START_REF] Mattes | PET-CT image registration in the chest using free-form deformations[END_REF] is used as the loss function.

The pipelines

We now present constructions for the pipelines whose benchmark is the core of this work. All registration pipelines presented here were scripted in using command-line tools from FSL version 5.0 [Smith et al., 2004] [START_REF] Avants | Symmetric di eomorphic image registration with cross-correlation: evaluating automated labeling of elderly and neurodegenerative brain[END_REF][START_REF] Avants | A reproducible evaluation of ANTs similarity metric performance in brain image registration[END_REF]. It pools registered data from all subjects. N.B.: SBRef = singleband high-resolution 3D volume.

As in [START_REF] Glasser | The minimal preprocessing pipelines for the Human Connectome Project[END_REF], all the transformations are postponed and the original 4D EPI is resampled at the end by applying the composition of these transformations in a single step. et al., 2009] antsRegistration, antsApplyTranforms, and some custom scripts (for distortion correction) from the HCP scripts described in [START_REF] Glasser | The minimal preprocessing pipelines for the Human Connectome Project[END_REF], hosted on Github. Except stated otherwise, all a ne registrations (motion correction, coregistration) were performed using FSL's irt tool [Smith et al., 2004] with Normalized Mutual Information as the cost function (option: -cost normmi).

Classical indirect T1-based method

The classical indirect T1-based pipeline for registration of EPI images can be schematized as follows 2 :

2 The "•" symbol denotes composition of transformations. EPI → templ. = (T1 → templ.) nonlinear

• (EPI 0 BBR → T1) linear •DistCorr (8.1)
in which a deformation of the subject's T1 is estimated and this deformation is then used to warp the same subject's EPI data. We implemented the pipeline as follows. Here, EPI 0 is any single-volume EPI image previouslycoregistered with the 4D EPI sequence. Typical choices include: the middle volume of the lm or the mean volume after motion correction. In our implementations, we used the former.

For the template, a subject is xed and its T1-weighted image is used as the template. For each other subject, (a) distortion correction is used to learn a nonlinear undistorting warp eld, in a procedure already described in subsection 8.2.2 above. Then, (b) motion correction is done to realign the the subject's EPI data to the mean thereof. The subject's T1-weighted image is then aligned to this mean EPI image via coregistration (an a ne transformation). We use BBR (boundary-based registration) [START_REF] Douglas | Accurate and robust brain image alignment using boundary-based registration[END_REF] for this coregistration step, for optimal results and fair comparison. BBR is a state-of-the-art functional-to-structural registration method driven by intensity di erence across boundary (samples). It uses white-matter boundaries (via T1w segmentation). BBR need good structural images (with little contrast bias), and some anatomical contrast in the EPI image (which is the case of the single-band high-resolution reference images in the HCP dataset [van Essen et al., 2012]). The implementation we use is epi_reg script of FSL [Smith et al., 2004]. However, since BBR is an a ne correction method, it still su ers from the limitations explained in the introductory section. In particular, it is not resiliant to distortions in the input EPI image.

(c) ANTs is used to learn a deformation of the T1 image to the template (which is a xed subject). This produces a warped version of the T1-weighted image, together with the corresponding deformation (and its inverse too), for passing from the subject's space to the template space. Finally, (d) the deformation above (T1-based), and all the other postponed warp elds, a ne transformations, etc., are then applied (in respective order) to all EPI data previously aligned (rigidly, via coregistration) with the T1-weighted image of the subject; these may include EPI images acquired on the same subject during another task, for instance. This one-step resampling procedure (see subsection 8.2.2) then produces a registered, motioncorrected, undistorted version of the input EPI data.

Then mean of all the registered T1-weighted images is computed, and becomes the template henceforth. This procedure is iterated a couple of times.

Our proposed direct EPI-based non-linear inter-subject registration method

Our proposed pipeline operates just as the classical indirect T1-based pipeline described above in 8.2.3, except that the anatomical image is replaced with the single-band high-resolution EPI (the SBRef) image, which has more tissue contrast than the any volume of the 4D EPI lm being registered [START_REF] Glasser | The minimal preprocessing pipelines for the Human Connectome Project[END_REF], and also does not su er from multi-band artifacts. The anatomical image is not used anywhere in this pipeline. The pipeline can be schematized as follows:

EPI → templ. = (EPI 0 → templ.) nonlinear •DistCorr, (8.2)
where we take EPI 0 = Single-band high-resolution (SBRef) EPI image.

A note on image interpolation (resampling) To avoid degrading the images as they travel through a pipeline, we stack all intermediate transformations and postpone the resampling operations to the end of the pipeline. The transformations are then concatenated (i.e composed), and applied to the input image in a one-step resampling procedure based on the Apply-Transforms tool of the ANTs software [START_REF] Avants | Symmetric di eomorphic image registration with cross-correlation: evaluating automated labeling of elderly and neurodegenerative brain[END_REF][START_REF] Avants | Advanced normalization tools (ANTS)[END_REF]. For example, a ne transformations estimated during the motion correction step are converted to warp elds using FSL's convertwarp tool [Smith et al., 2004]. FSL's applywarp tool [Smith et al., 2004] is then used to jointly apply this a ne transformation warp elds and the warp elds corresponding to the deformations estimated by topup [Smith et al., 2004], which are then stacked with subsequent transformations. We use this strategy in both pipelines.

Relation to previous works 8.3.1 Direct EPI-to-EPI non-linear inter-subject registration

The idea of EPI-to-EPI registration has already been suggested in the literature. For example, the method in [START_REF] Grabner | A study-speci c fMRI normalization approach that operates directly on high resolution functional EPI data at 7 Tesla[END_REF] used high-resolution EPI (1.1mm isotropic) data for di erent subjects acquired at 7T to iteratively build a sequence of EPI-based study-speci c templates of increasing quality / resolution [START_REF] Grabner | Symmetric atlasing and model based segmentation: an application to the hippocampus in older adults[END_REF]. The nest of these templates shows a great deal of anatomical detail. Group-level activation patterns for a ngertapping task were also shown to be very accurately localized on the posterior bank of the central sulcus. The authors concluded that high-resolution (7T) EPI images contain enough anatomical information for inter-subject registration, and so one can e ectively by-pass the anatomical image of subjects in pipeline. This would for example allow one to avoid the classical coregistration step used to align the subject's EPI images to their anatomy. Our experiments con rm and extend the ndings of [START_REF] Grabner | A study-speci c fMRI normalization approach that operates directly on high resolution functional EPI data at 7 Tesla[END_REF]], but at an even lower resolution: 2mm resolution, obtained from 3T MRI, and on a much larger dataset. Indeed, using a much larger bail of 110 subjects, from Human Connectome Project (HCP) dataset [van Essen et al., 2012], and a variety of di erent task contrasts, we show that registration with our pipeline increases the pairwise NMI of subjects, over the classical pipeline; crucially, this leads to a decrease in residual post-registration inter-subject misalignement.

In comparison, the pipeline we propose (refer to 8.2.3) is much lighter computationally as we bypass the potentially expensive and challenging step of generating a good template from EPI data [START_REF] Grabner | Symmetric atlasing and model based segmentation: an application to the hippocampus in older adults[END_REF]. Of course, this economy is more of a compromise between complexity and accuracy, and might be potential limitation of our contribution. Finally, we note that the work in [START_REF] Grabner | A study-speci c fMRI normalization approach that operates directly on high resolution functional EPI data at 7 Tesla[END_REF] did not consider the distortion problem which are severe even at 3T [START_REF] Andersson | How to correct susceptibility distortions in spin-echo echo-planar images: application to di usion tensor imaging[END_REF], as it is the case with the HCP data.

Non-linear EPI-to-structural coregistration

A recent work [START_REF] Wang | Evaluation of Field Map and Nonlinear Registration Methods for Correction of Susceptibility Artifacts in Di usion MRI[END_REF] has considered the possibility of replacing the classical linear EPI-to-structural coregistration step with a non-linear counterpart, and then running a non-linear structural-to-template registration as usual. They show that their method outperforms the method based on distortion correction and linear EPI-to-structural coregistration followed by structural-to-template registration as usual (see 8.2.3). In contrast, our proposed method (refer to 8.2.3) does not use the structural image at all.

Experiments

We now describe benchmarks done to compare the pipelines presented in this paper (subsection 8.2.3) on the task fMRI data of 110 subjects from the HCP dataset [van Essen et al., 2012]. The task fMRI data were acquired in an attempt to assess major domains that sample the diversity of neural systems, including: 1) visual, motion, somatosensory, and motor systems; 2) language processing (semantic and phonological processing); 3) social cognition (Theory of Mind); and 4) emotion processing. Due to time constraints, our benchmarks were run only on these 4 (out of a total of 7) tasks (i.e protocols). Also, only data for LR (left-right) phase-encoding direction [START_REF] Chang | A technique for accurate magnetic resonance imaging in the presence of eld inhomogeneities[END_REF]] runs were used. In all the non T1-based pipelines, the single-band high-resolution (SBRef) image of the motor task was used to learn deformations of the subject's brain into template space (a xed subject of the same dataset). The estimated deformations were then applied to warp EPI data (previously coregistered to same subject's motor SBRef) acquired on the same subject during di erent task conditions, into template space. GLMs (General Linear Models) [Friston et al., 1994] were run using nipy [START_REF] Gorgolewski | Nipype: a exible, lightweight and extensible neuroimaging data processing framework in python[END_REF], open-source Python library for analysis of neuro-imaging data. For the purpose of reporting the results, the resulting maps were co-registered to MNI space a posteriori.

Evaluation metrics

The pipelines were evaluated using the following qualitative and quantitative metrics.

Normalized mutual information evaluation (NMI)

NMI (see Table 1 for de nition) is a popular similarity metric used to assess the quality of registration between two images, i.e how well aligned the images are to one another (for example [START_REF] Maes | Multimodality image registration by maximization of mutual information[END_REF]). It is also the loss function minimized by many optimization algorithms in image registration. A detailed overview of the use of the NMI metric in medical image registration can be found in [START_REF] Pluim | Mutual-information-based registration of medical images: a survey[END_REF]. In our experiments, FSL's irt 3 tool [Smith et al., 2004] was used to compute NMI.

3 With the "-schedule" option.

Inter-subject residual variance

In a good registration method, the residual subject-to-subject variance of the EPI image should be reduced. The aim of inter-subject registration is indeed to put subjects into spatial correspondence to facilitate later group analysis. To measure the quality of the di erent registration methods in this regards, we computed the coe cient of variation (CoV) across the di erent subjects after registration. This is de ned by CoV = variance image across subjects mean image across subjects . (8.3)

High values in this 3D image would outline regions of the brain which are not well registered across subjects.

Group-level statistics and functional brain network pa erns

Finally, in a successful inter-subject registration procedure, we expect the functional activation patterns to be more localized in space and to have higher peaks. Or could this e ect be masked by inter-subject variability in activation magnitude ? This will be discussed in detail in the discussion section 8.6.

How many (plausible) pipelines are there ?

It is worth noting that there are potentially hundreds of pipelines which could have considered for testing: should we do distortion correction ? And if yes, how ? Should we use linear or nonlinear model for the deformation eld ? What degree should we use for the interpolating splines ? In fact as noted in [START_REF] Poldrack | Scanning the Horizon: Future challenges for neuroimaging research[END_REF], there are exponentially many pipelines that can be considered, based on the answers to the above choices. Of course some of these parameters have rule-of-thumb default values (for example, there is no doubt distortion correction is a good thing to do), but others are open to preferential choice. Thus our goal is not to consider all possible pipelines, but to look at a more focal question: does direct EPI-based interregistration outperform the traditional indirect T1-based pipeline ?

Results

We now present results of experiments performed on the task fMRI protocols of the HCP dataset [van Essen et al., 2012]. Refer to section 8.4 for detailed information about the experiments we did. The di erent pipelines discussed in section 8.2.3 were used to register the data (inter-subject registration), and the quality of the registration was benchmarked using the di erent evaluation metrics discussed in Section 8.4.

Normalized Mutual Information (NMI)

The results comparing across-subject NMI for the pipelines are presented in Figure 8.3. We see that MNI is in most cases higher through our approach, which implies that our proposed direct EPI-based pipeline mildly outperforms the classical indirect T1-based pipeline.

Figure 8.3: Normalized Mutual Information -NMI (higher values are better). Each point (x, ) on the plots such that x is the NMI of a given pair of subjects registered using the pipeline on the abscissa and is the NMI of the same pair of subjects registered using the pipeline on the ordinate. From the one-sided We see that our proposed direct EPI-based pipeline signi cantly outperforms the classical indirect T1-based pipeline.

Residual inter-subject spatial variability

In Figure 8.4, we show across-subject histograms of across-subject per-voxel Coe cient of Variation (small is better). We see that our proposed direct method outperforms the classical indirect T1-based method, as the former leads to relatively more mis-aligned voxels across subjects, most concentrated on the outer edge of the cortex (see Figure 8.4 (a)).

ality of estimated EPI group template

To compare the quality of the group template produced by either pipeline, a snapshot of the resulting mean image or template is displayed in Figure 8.5. Compared to the our proposed direct method, the mean image (across all subjects) from the indirect T1-based pipeline is blurry and has "ripples" on the cortical surface, indicative of residual mismatch between subjects after registration. The across-subject mean image post-registration with our direct EPI-based pipeline is the sharpest, showing that the subjects have Coe cient-of-Variation (CoV) after registration. Top: Log10 of ratio of across-subject Coe cient of Variation (CoV) for indirect T1-based pipeline / direct EPI-based. We see that the gain of our proposed method is most pronounced along the cortical surface.

Bottom: Histograms of CoV for both pipelines. Again, we see clearly that our proposed method reduces the inter-subject variability by a much larger margin, indicative of improved subject-to-subject alignment.

been matched extremely well. Also, one notices that the mean image from the indirect T1-based pipeline still has some residual distortion (here in the left-to-right direction), even though distortion correction was done as part of both pipelines. Also the mean image from the T1-based pipeline has ripples on the cortical surface indicative of residual registration problems, which can be attributed residual EPI-distortions that could not be captured via coregistration.

Group-level statistics and Functional brain network pa erns

As regards group-level GLM scores, we see from Figure 8.6 that our proposed method does just as good as the classical indirect T1-based pipeline. Here we see that our proposed direct EPI-based registration scheme leads to slightly higher across-subject mean activation peaks. For each contrast, a cut has been made around the location of the activation peak, to display curves of the activation prole and across-subject variability thereof, in a neighborhood of this peak location. This is remarkable, as the former pipeline does not use any anatomical data. However, as noted in [Thirion et al., 2007, Thyreau et al., 2012], the intersubject variability in GLM results is not due to misregistration, but intrinsic subject di erences with a more physiological nature: the response of subjects to the same stimulus / task is modulated di erently, and is more dependent on e ect size uctuations than position. This is con rmed in the curves in Figure 8.6, where we can see that the spatial across-subject activation pro les are very similar between the compared registration methods, except for the already noted slight improvement of the peak mean activation pattern obtained by our proposed method.

Finally, Figure 8.7 comparing the functional brain networks obtained by running ICA on the images registered with each of the pipelines, shows essentially the same network patterns. The absence of a di erence between these maps can be explained by the fact that resting state networks are less focal than task-based activation-patterns, and so the former are less sensitive to the quality of the underlying registration procedure.

Discussion and concluding remarks

Classical inter-subject registration pipelines use the T1-weighted (anatomical) image of a subject to estimate the subject-to-template warp. An obvious issue is that high-quality T1-weighted images are not always available, but more generally, it is not always possible to completely align the EPI images of a subject to their T1-weighted image via coregistration. Added to this is the possibility that such an intermediate registration step is a potential source of interpolation artifacts, not to mention the added computational cost (which may exceed the rest of the computation time by many orders of magnitude, for example, in the case of surface-based methods). As shown by our experiments on the HCP dataset [van Essen et al., 2012] (Figure 8.1), this is for example the case in the presence of distortions [START_REF] Wan | Reduction of geometric and intensity distortions in echo-planar imaging using a multireference scan[END_REF], Mangin et al., 2002[START_REF] Zeng | Image distortion correction in EPI: comparison of eld mapping with point spread function mapping[END_REF][START_REF] Andersson | How to correct susceptibility distortions in spin-echo echo-planar images: application to di usion tensor imaging[END_REF]] that persist even after correction. Further, as noted in [START_REF] Yamada | E cacy of distortion correction on di usion imaging: comparison of FSL eddy and eddy_correct using 30 and 60 directions di usion encoding[END_REF], distortions cannot be separated with registration based techniques, which are compensatory operations. Consequently, the e cient distortion correction method in EPI data remains an open question. Our work proposes a direct EPI-based inter-subject registration pipeline that to some extent evades these bottlenecks.

We have proposed a computationally cheap EPI-based pipeline for direct inter-subject nonlinear registration of functional data. Our method has been empirically validated on the HCP dataset [van Essen et al., 2012], where we have shown that we obtain registered subject images with less inter-subject variability. Such direct EPI-based methods should replace the well-accepted classical T1-based strategy. Results on the HCP dataset [van Essen et al., 2012] show that the proposed pipeline outperforms the classical T1-based indirect registration strategy, according to a variety of di erent quality metrics: Normalized Mutual Information -NMI (Figure 8.3), residual inter-subject variance (Figure 8.4), and quality of estimated group template (Figure 8.5), without compromising the quality of post-registration statistical analyses results (GLM, ICA, etc.). These results replicate the ndings of [START_REF] Grabner | A study-speci c fMRI normalization approach that operates directly on high resolution functional EPI data at 7 Tesla[END_REF] on a larger dataset (110 subjects, compared to 10 in the reference paper) and in a 3T setting. [START_REF] Raichle | A default mode of brain function[END_REF], using MNI coordinates reported in Table 1 of [START_REF] Watanabe | A pairwise maximum entropy model accurately describes resting-state human brain networks[END_REF]. Remarkably, we observed that according to low-level QA metrics like NMI (Figure 8.3), residual inter-subject spatial variability (Figure 8.4) and the quality of across-subject mean registered EPI image (Figure 8.5), our proposes method outperforms the classical indirect T1-based registration. In terms of more high-level metrics like group-level GLM statistics, these gains though still present, are as not as pronounced (refer to Figure 8.6). Indeed, as noted in [Thirion et al., 2007, Thyreau et al., 2012[START_REF] Xu | Modeling inter-subject variability in fMRI activation location: A Bayesian hierarchical spatial model[END_REF], the inter-subject variability in GLM results is not due to misregistration, but intrinsic subject di erences with a more physiological nature: the size of e ects and the anatomical localization are subject-speci c. In chapter 11, we show that resting-state fMRI data can be used to predict the activation maps of a subject to a task, with an R 2 -score which can be up to 0.5 for some subjects and task. This is an enhancement on previous work by [Tavor et al., 2016], and shows di erences in task-based brain activations are largely physiological -in contrast to being driven by subjects' brain morphological di erences-and can be predicted from resting state fMRI data.

In a separate work [Dohmatob et al., 2016], also presented in chapter 9 in detail, we have considered the possibility of explicitly modeling this physiology di erences by estimating latent factors of variability across-subjects in a data-driven way using dictionary-learning. The motivating idea behind such a model, is that activation across-subjects would be governed by the same generative model (the latent model), and modulated on the subjectlevel by subject-speci c physiology.

Learning pa erns of inter-subject functional variability from data

I

, inter-subject variability is often handled as a statistical residual and discarded. Yet there is evidence that it displays structure and contains important information. Univariate models are ine ective both computationally and statistically due to the large number of voxels compared to the number of subjects. Likewise, statistical analysis of weak effects on medical images often relies on de ning regions of interests (ROIs). For instance, pharmacology with Positron Emission Tomography (PET) often studies metabolic processes in speci c organ sub-parts that are de ned from anatomy. Population-level tests of tissue properties, such as di usion, or simply their density, are performed on ROIs adapted to the spatial impact of the pathology of interest. In functional brain imaging, e.g functional magnetic resonance imaging (fMRI), ROIs must be adapted to the cognitive process under study, and are often de ned by the very activation elicited by a closely related process [START_REF] Saxe | Divide and conquer: a defense of functional localizers[END_REF]. ROIs can boost statistical power by reducing multiple comparisons that plague image-based statistical testing. If they are de ned to match spatially the di erences to detect, they can also improve the signal-to-noise ratio by averaging related signals. However, the crux of the problem is how to de ne these ROIs in a principled way. Indeed, standard approaches to region de nition imply a segmentation step. Segmenting structures on rst-level statistical maps, as in fMRI, typically yields meaningful units, but is limited by the noise inherent to these maps. Relying on a di erent imaging modality hits cross-modality correspondence problems.

Introduction and sketch of our contributions

I

, we propose to use the variability of the statistical maps across the population to de ne regions. This idea is reminiscent of clustering approaches, that have been employed to de ne spatial units for quantitative analysis of information as diverse as brain ber tracking, brain activity, brain structure, or even imaging-genetics. See [START_REF] Varol | Supervised block sparse dictionary learning for simultaneous clustering and classi cation in computational anatomy[END_REF]Davatzikos, 2014, Hibar et al., 2013] and references therein. The key idea is to group together features -voxels of an image, vertices on a mesh, ber tracts-based on the quantity of interest, to create regions -or ber bundles-for statistical analysis. However, unlike clustering that models each observation as an instance of a cluster, we use a model closer to the signal, where each observation is a linear mixture of several signals. The model is closer to mode nding, as in a principal component analysis (PCA), or an independent component analysis (ICA), often used in brain imaging to extract functional units [START_REF] Beckmann | Probabilistic independent component analysis for functional magnetic resonance imaging[END_REF]]. Yet, an important constraint is that the modes should be sparse and spatially-localized. For this purpose, the problem can be reformulated as a linear decomposition problem like ICA/PCA, with appropriate spatial and sparse penalties [Varoquaux et al., 2011[START_REF] Abraham | Extracting brain regions from rest fMRI with Total-Variation constrained dictionary learning[END_REF].

We propose a multivariate online dictionary-learning method for obtaining decompositions with structured and sparse components (aka atoms). Sparsity is to be understood in the usual sense: the atoms contain mostly zeros. This is imposed via an 1 penalty on the atoms. By "structured", we mean that the atoms are piece-wise smooth and compact, thus making up blobs, as opposed to scattered patterns of activation. We impose this type of structure via a Laplacian penalty 1 on the dictionary atoms. Combining 1 This is a slight abuse of language as we really mean the Sobolev semi-norm v → v T ∆v = (∇v) T ∇v ≥ 0 and not the Laplacian linear operator ∆ := ∇ T ∇.

the two penalties, we therefore obtain decompositions that are closer to known functional organization of the brain. This non-trivially extends the online dictionary-learning / dictionary-learning work [Mairal et al., 2010], with only a factor of 2 or 3 on the running time. By means of experiments on a large public dataset, we show the improvements brought by the spatial regularization with respect to traditional 1 -regularized dictionary learning.

We also provide a concise study of the impact of hyper-parameter selection on this problem and describe the optimality regime, based on relevant criteria (reproducibility, captured variability, explanatory power in prediction problems).

Applied to functional brain imaging, it separates successfully activation maps into localized units of brain activity. Our contribution is to frame spatial penalties as a particular case of more general Laplacian regularization and introduce an e cient online algorithm for dictionary-learning in these settings. Applied to functional brain imaging, it separates successfully activation maps into localized units of brain activity. Here we do things that are closer to probabilistic segmentations 9.2 Smooth Sparse Online Dictionary-Learning (Smooth-SODL)

Consider a stack X ∈ R n×p of n subject-level brain images x 1 , x 2 , . . . , x n each of shape n 1 × n 2 × n 3 , seen as p-dimensional row vectors -with p = n 1 × n 2 × n 3 , the number of voxels. These could be images of fMRI activity patterns like statistical parametric maps of brain activation, raw pre-registered (into a common coordinate space) fMRI time-series, PET images, etc. We would like to decompose these images as a product of k ≤ min(n, p) component maps (aka latent factors or dictionary atoms) d 1 , . . . , d k ∈ R p×1 and modulation coe cients c 1 , . . . , c n ∈ R k ×1 called codes (one k-dimensional code per sample point), i.e

x i ≈ Dc i , for i = 1, 2, . . . , n (9.1)
where

D := [d 1 | . . . |d k ] ∈ R p×k ,
an unknown dictionary to be estimated.

Figure 9.1: Dictionary-learning with smoothness and sparsity constraints on the atoms. On the right, each time point is a masked 3D brain image and corresponds to a sample, and each voxel corresponds to a feature, giving an p-byn matrix X, where n is the number of samples and p is the number of features. On the right of the equation, the sought-for p-by-k dictionary D is a low-dimensional representation these images, by means of a (non-orthonormal) basis of k min(n, p) smooth and sparse 3D brain images called atoms. In this representation, each sample point (i.e 3D brain image) X i ∈ R p is mapped onto k-dimensional vector c i , called the code of X i .

Typically, p ∼ 10 5 -10 6 (in full-brain high-resolution fMRI) and n ∼ 10 2 -10 5 (for example, in considering all the 500 subjects and all the about 15 -20 functional tasks of the Human Connectome Project dataset [van Essen et al., 2012]). Our work handles the extreme case where both n and p are large (massive-data setting). It is reasonable then to only consider undercomplete dictionaries: k ≤ min(n, p). Typically, we use k ∼ 50 or 100 components. It should be noted that online optimization is not only crucial in the case where n/p is big; it is relevant whenever n is large, leading to prohibitive memory issues irrespective of how big or small p is.

As explained in section 9.1, we want the component maps (aka dictionary atoms) d j to be sparse and spatially smooth (illustrated in Fig. 9.1). A principled way to achieve such a goal is to impose a boundedness constraint on 1 -like norms of these maps to achieve sparsity and simultaneously impose smoothness by penalizing their Laplacian. Thus, we propose the following penalized dictionary-learning model min

D∈R p×k lim n→∞ 1 n n i=1 min c i ∈R k 1 2 x i -Dc i 2 2 + 1 2 α c i 2 2 + γ k j=1
Lap(d j ).

subject to d 1 , . . . , d k ∈ C (9.
2)

The ingredients in the model can be broken down as follows:

• Each of the terms max

c i ∈R k 1 2 x i -Dc i 2
2 measures how well the current dictionary D explains data x i from subject i. The Ridge penalty term ϕ (c i ) ≡ 1 2 α c i 2 2 amounts to placing an isotropic Gaussian prior on the codes, namely

p(c i ) ∝ exp(-1 2 α c i 2
2 ). In the context of a speci c neuro-imaging problem, if there are good grounds to assume that each sample / subject should be sparsely encoded across only a few atoms of the dictionary, then we can use the 1 penalty ϕ (c i ) := α c i 1 as in [Mairal et al., 2010], corresponding to a Laplace prior on the codes, namely p(c i ) ∝ exp(-α c i 1 ). We note that in contrast to the 1 penalty, the Ridge leads to stable codes. The parameter α > 0 controls the amount of penalization on the codes.

• The constraint set C is a sparsity-inducing compact simple 2 convex sub-2 Mainly in the sense that the Euclidean projection onto C should be easy to compute. set of R p like an 1 -ball B p, 1 (τ ) or a simplex S p (τ ), de ned respectively as

B p, 1 (τ ) := d ∈ R p s.t |d 1 | + |d 2 | + . . . + |d p | ≤ τ , and S p (τ ) := B p, 1 (τ ) ∩ R p + .
Other choices (e.g ElasticNet ball) are of course possible. The radius parameter τ > 0 controls the amount of sparsity: smaller values lead to sparser atoms. The Laplacian regularization Lap (see Table 1 for de nition) is meant to impose blobs. γ ≥ 0 controls how much regularization we impose on the atoms, compared to the reconstruction error.

The above formulation, which we dub Smooth Sparse Online Dictionary-Learning (Smooth-SODL) is inspired by, and generalizes the standard dictionarylearning framework of [Mairal et al., 2010] -henceforth referred to as Sparse Online Dictionary-Learning (SODL); setting γ = 0, we recover SODL [Mairal et al., 2010].

Algorithms

The objective function in problem (9.2) is separately convex and blockseparable w.r.t each of C and D but is not jointly convex in (C, D). Also, it is continuously di erentiable on the constraint set, which is compact and convex. Thus by classical results (e.g Bertsekas [START_REF] Bertsekas | Nonlinear programming[END_REF]), the problem can be solved via Block-Coordinate Descent (BCD) [Mairal et al., 2010]. Reasoning along the lines of [START_REF] Jenatton | Structured sparse principal component analysis[END_REF], we derive that the BCD iterates are as given in Alg. 3. A crucial advantage of using a BCD scheme is that it is parameter free: there is not step size to tune. The resulting algorithm Alg. 3, is adapted from [Mairal et al., 2010]. It relies on Alg. 4 for performing the structured dictionary updates, the details of which are discussed below.

Algorithm 3: Online algorithm for the dictionary-learning problem (9.2) Require: Regularization parameters α, γ > 0; initial dictionary D ∈ R p×k , number of passes / iterations T on the data.

1: A 0 ← 0 ∈ R k×k , B 0 ← 0 ∈ R p×k (historical "su cient statistics") 2: for t = 1 to T do 3:
Empirically draw a sample point x t at random.

4:

Code update: Ridge-regression (via SVD of current dictionary D)

c t ← arg min c∈R k 1 2 x t -Dc 2 2 + 1 2 α c 2 2 . ( 9.3) 5:
Rank-1 updates:

A t ← A t -1 + c t c T t , B t ← B t -1 + x t c T t 6:
BCD dictionary update: Compute update for dictionary D using Alg. 4.

7: end for

Update of the codes: Ridge-coding. The Ridge sub-problem for updating the codes

c t = (D T D + αI) -1 D T x t (9.4)
is computed via an SVD of the current dictionary D. For α ≈ 0, c t reduces to the orthogonal projection of x t onto the image of the current dictionary D.

As in [Mairal et al., 2010], we speed up the overall algorithm by sampling mini-batches of η samples x t , . . . , x η and compute the corresponding codes c 1 , c 2 , ..., c η at once. We typically use we use mini-batches of size η ∼ 20 images.

BCD dictionary update for the dictionary atoms. Let us de ne timevarying matrices A t := t i=1 c i c T i ∈ R k ×k and B t := t i=1 x i c T i ∈ R p×k , where t = 1, 2, . . . denotes time. We x the matrix of codes C, and for each j, consider the update of the jth dictionary atom, with all the other atoms d k j kept xed. The update for the atom d j can then be written as

d j = arg min d j ∈ C 1 t t i=1 1 2 x i -Dc i 2 2 + γ Lap(d j )
= arg min

d j ∈ C t i=1 1 2 x i -Dc i 2 2 + γtLap(d j )
= arg min

d j ∈ C F γ (a j,j /t ) -1 (d j , d j old + a -1 j,j (b j -Da j ) see chapter 10 below ), (9.5)
where

F γ (d j , a) ≡ 1 2 d j -a 2 2 + γ Lap(d j ) = 1 2 d j -a 2 2 + 1 2 ∇d j 2 . (9.6)
Problem (9.5) is thus a compactly-constrained minimization of the 1-stronglyconvex quadratic functions F γ (., a) : R p → R de ned above. This problem can further be identi ed with a denoising instance (i.e in which the design Algorithm 4: BCD dictionary update with Laplacian prior

Require: D = [d 1 | . . . |d k ] ∈ R p×k (input dictionary), 1: A t = [a 1 | . . . |a k ] ∈ R k ×k , B t = [b 1 | . . . |b k ] ∈ R p×k (history) 2:
while stopping criteria not met, do 3:

for j = 1 to r do 4:

Fix the code C and all atoms k j of the dictionary D and then update d j as follows

d j ← arg min d j ∈ C F γ (a j,j /t ) -1 (d j , d j old + a -1 j,j (b j -Da j )) (9.7)
(See below for details on the derivation and the resolution of this problem) 5:

end for 6: end while matrix or deconvolution operator is the identity operator) of the GraphNet model [Grosenick et al., 2013, Hebiri andvan de Geer, 2011]. Fast rst-order methods like FISTA [Beck and Teboulle, 2009] with optimal rates O(L/ √ ϵ ) are available 3 for solving such problems to arbitrary precision ϵ > 0. One 3 For example, see [Dohmatob et al., 2014, Varoquaux et al., 2015], implemented as part of the Nilearn open-source Python library [Abraham et al., 2014].

computes the Lipschitz constant to be L F γ (.,a) ≡ 1 + γ L Lap = 1 + 4D γ , where as before, D is the number of spatial dimensions with D = 3 for volumetric images. One should also mention that under certain circumstances, it is possible to perform the dictionary updates in the Fourier domain, via FFT. This alternative approach is developed in the Appendix of [Dohmatob et al., 2016]. Finally, one notes that, since constraints in problem (9.2) are separable in the dictionary atoms d j , the BCD dictionary-update algorithm Alg. 4 is guaranteed to converge to a global optimum, at each iteration [START_REF] Bertsekas | Nonlinear programming[END_REF], Mairal et al., 2010].

How di icult is the dictionary update for our proposed model ? A favorable property of the vanilla dictionary-learning [Mairal et al., 2010] is that the BCD dictionary updates amount to Euclidean projections onto the constraint set C, which can be easily computed for a variety of choices (simplexes, closed convex balls, etc.). One may then ask: do we retain a comparable algorithmic simplicity even with the additional Laplacian terms Lap(d j ) ? The short answer is yes: empirically, we found that 1 or 2 iterations of FISTA [Beck and Teboulle, 2009] are su cient reach an accuracy of 10 -6 in problem (9.5), which is su cient to obtain a good decomposition in the overall algorithm. However, choosing γ "too large" will provably cause the dictionary updates to eventually take forever to run. Indeed, the Lipschitz constant in problem (9.5) is L t = 1 + 4Dγ (a j,j /t ) -1 , which will blow-up (leading to arbitrarily small step-sizes) unless γ is chosen so that

γ = γ t = O max 1≤j ≤k a j,j = O max 1≤j ≤k t i=1 C j 2 2 /t = O( A t ∞,∞ /t ). (9.8)
Finally, the Euclidean projections onto the 1 ball C can be computed exactly in linear-time O(p) (see for example [START_REF] Condat | Fast projection onto the simplex and the 1 ball[END_REF][START_REF] Duchi | Efcient projections onto the l 1-ball for learning in high dimensions[END_REF]). The dictionary atoms j are repeatedly cycled and problem (9.5) solved. All in all, in practice we observe that a single iteration is su cient for the dictionary update sub-routine in Alg. 4 to converge to a qualitatively good dictionary.

Implementation and practical considerations

So ware implementation. All aspects of the code where implemented in the Python programming language. For the implementation of the proposed Alg. 3, we implemented a modi ed version of scikit-learn Python library's dict_learning module, to handle more general constraint sets and more general penalties both for the codes c i and for the dictionary atoms d j . The projection onto the 1 -ball C was coded in Cython, a toolkit for writing Python code to run at the speed of the C language.

Convergence of the overall algorithm. The Convergence of our algorithm (to a local optimum) is guaranteed since all hypotheses of [Mairal et al., 2010] are satis ed. For example, assumption (A) is satis ed because fMRI data are naturally compactly supported. Assumption (C) is satis ed since the ridge-regression problem (9.3) has a unique solution. More details are provided in the Appendix of [Dohmatob et al., 2016].

Practical considerations

Hyper-parameter tuning. Parameter-selection in dictionary-learning is known to be a di cult unsolved problem [Mairal et al., 2010[START_REF] Jenatton | Structured sparse principal component analysis[END_REF], and our proposed model (9.2) is not an exception to this rule. We did an extensive study of the quality of estimated dictionary varies with the model hyper-parameters (α, γ , τ ). The data experimental setup is described in Section 11.5. The results are presented in Fig. 9.2. We make the following observations: Taking the sparsity parameter τ in (9.2) too large leads to dense atoms that perfectly explain the data but are not very intepretable. Taking it too small leads to overly sparse maps that barely explain the data. This normalized sparsity metric (small is better, ceteris paribus) is de ned as the mean ratio d j 1 / d j 2 over the dictionary atoms. Left: Percentage explained variance of the decomposition, measured on left-out data split. Right: Average normalized sparsity of the dictionary atoms.

Concerning the α parameter, inspired by [START_REF] Ying | Accelerated, parallel and proximal coordinate descent[END_REF], we have found the following time-varying data-adaptive choice for the α parameter to work very well in practice:

α = α t ∼ t -1/2 . (9.9)
Likewise, care must be taken in selecting the Laplacian regularization parameter γ . Indeed taking it too small amounts to doing vanilla dictionarylearning model [Mairal et al., 2010]. Taking it too large can lead to degenerate maps, as the spatial regularization then dominates the reconstruction error (data delity) term. We nd that there is a safe range of the parameter pair (γ , τ ) in which a good compromise between the sparsity of the dictionary (thus its intepretability) and its explanation power of the data can be reached. See Fig. 9.2. K-fold cross-validation with explained variance metric was retained as a good strategy for setting the Laplacian regularization γ parameter and the sparsity parameter τ .

Initialization of the dictionary. Problem (9.2) is non-convex jointly in (C, D), and so initialization might be a might be a crucial issue. However, in our experiments, we have observed that even randomly initialized dictionaries eventually produce sensible results that do not jitter much across di erent runs of the same experiment.

Interlude: Working in the Fourier domain (when possible)

To close this section, let us point out a few special instances cases of problem (9.6), for peculiar choices of the constraint set Q. First note that the objective in problem ( 7) can be conveniently rewritten as .11) We note that the matrix-inversion (Iγ ∆) -1 that appears in the formula above is a Laplacian lter, and can be e ciently applied in closed-form (i.e non-iteratively) in the Fourier / frequency domain. Indeed, under periodic boundary conditions, the discrete Laplacian ∆ is Block-Circulant with Circulant Blocks (BCCB) and so is diagonalizable in the Fourier domain. Precisely, ∆ = F * ΛF (9.12)

F γ A t [j,j] -1 (v, v j + A t [j, j] -1 (vA j -B j t )) = 1 2 (v -ṽj ) T (I -γ A t [j, j] -1 ∆)(v -ṽj ) = 1 2 ( v -vj ) T (I -γ A t [j, j] -1 ∆)( v -vj ), (9.10)
with ṽj := (A t [j, j]I -γ ∆) -1 v j + A t [j, j] -1 (vA j -B j t ) . ( 9
where the complex orthonormal operator F represents the fast Discrete Fourier Transform (DFT), and Λ is diagonal matrix made p eigenvalues (including multiplicities) of the Laplace operator ∆, given by

Λ(ω) := - 3 d =1 2 sin ω d π 2n d 2 = -2 3 d =1 1 -cos ω d π n d ≤ 0, for ω = (ω 1 , ω 2 , ω 3 ) ∈ [[0, n 1 -1]] × [[0, n 2 -1]] × [[0, n 3 -1]].
We note that the spectral norm of Laplace operator in D dimensions (here

D = 3) is ∆ 2 = γmax (-∆) = 2 × D × (1 + 1) = 4D.
Now, one can then harvest the closed-form solution

(I -γ ∆) -1 a = (F -1 (I -γ Λ) -1 F )(a) = F -1 (s) , (9.13)
where s ∈ R p is de ned by s(ω

) := â(ω) 1 -γ ∆(ω)
, with â := F (a). These DFT computations have complexity O(p log p).

For applying the DFTs above, one can use the FFTW 4 -or Fastest Fourier 4 FFTW is generally taught to be one of the fastest implementations of the FFT, yielding up to 3× speedup against competing libraries like LAPACK.

Transform in the West-library for computing the forward and inverse Fourier transforms needed to apply the Laplacian lter.

Pure 2 constraint. Here, the constraint set C is an L2 ball (with radius = 1, w.l.o.g) in R 2 . By the Rayleigh energy theorem (aka Parseval's identity for the DFT), one has

v 2 = p v 2 2 , ∀v ∈ R p
and so problem ( 7) can be written as

v j ← arg min v∈R p , v 2 2 ≤1 1 2 ( v -vj ) * (I -γ A t [j, j] -1 Λ)( v -vj ) = F * arg min v∈C p , v 2 2 ≤p 1 2 ( v -vj ) * (I -γ A t [j, j] -1 Λ)( v -vj ) = F * P E ( vj ) (9.14)
where

E := (I -γ A t [j, j] -1 Λ) 1 2 v s.t v ∈ C p , v 2 2 ≤ p , (9.15)
a hyper-ellipsoid in standard position (i.e 0-centered and axes-aligned). Using elementary geometric arguments, one can show that the projection P E ( vj ) can be computed e ciently using a kind of root-nding algorithm [START_REF] Dai | Fast Algorithms for Projection on an Ellipsoid[END_REF], and converges exponentially fast.

Non-negative Lasso. In case the constraint set C for the dictionary atoms is a simplex S p (τ ), the simplex (see section 9.2), then the BCD update for the jth atom becomes

v j ← arg min v∈R p , v≥0, 1 T v≤1 1 2 ( v -vj ) * (I -γ A t [j, j] -1 Λ)( v -vj ) = F * arg min v∈C p , -F * v≤0, 1T v≤1 1 2 ( v -vj ) * (I -γ A t [j, j] -1 Λ)( v -vj ) , (9.16)
which is a diagonal quadratic program with linear constraints, and can be e ectively solved via the well-known simplex method, for example.

Related works

While there exist algorithms for online sparse dictionary-learning that are very e cient in large-scale settings (for example [Mairal et al., 2010], or more recently [Mensch et al., 2016]) imposing spatial structure introduces couplings in the corresponding optimization problem [Dohmatob et al., 2014]. So far, spatially-structured decompositions have been solved by very slow alternated optimization [Varoquaux et al., 2011[START_REF] Abraham | Extracting brain regions from rest fMRI with Total-Variation constrained dictionary learning[END_REF]. Notably, structured priors such as TV-1 [Baldassarre et al., 2012] minimization, were used by [START_REF] Abraham | Extracting brain regions from rest fMRI with Total-Variation constrained dictionary learning[END_REF] to extract data-driven stateof-the-art atlases of brain function. However, alternated minimization is very slow, and large-scale medical imaging has shifted to online solvers for dictionary-learning like [Mairal et al., 2010] and [Mensch et al., 2016]. These do not readily integrate structured penalties. As a result, the use of structured decompositions has been limited so far, mostly due to the computational cost of the ensuing algorithms. Our approach instead uses a Laplacian penalty to impose spatial structure at a very minor cost and adapts the online-learning dictionary-learning framework [Mairal et al., 2010], resulting in a fast and scalable structured decomposition. Second, the approach in [START_REF] Abraham | Extracting brain regions from rest fMRI with Total-Variation constrained dictionary learning[END_REF] though very novel, is heuristic, as it does not come with theoretical guarantees. In contrast, our method enjoys the same convergence guarantees and comparable numerical complexity as the basic unstructured online dictionary-learning [Mairal et al., 2010].

Finally, one should also mention [START_REF] Varoquaux | Cohort-level brain mapping: learning cognitive atoms to single out specialized regions[END_REF] that introduced an online group-level functional brain mapping strategy for di erentiating regions re ecting the variety of brain network con gurations observed a the population, by learning a sparse representation of these in the spirit of [Mairal et al., 2010].

Experiments

Setup. Our experiments were done on task fMRI data from 500 subjects from the HCP -Human Connectome Project-dataset [van Essen et al., 2012]. These task fMRI data were acquired in an attempt to assess major domains that are thought to sample the diversity of neural systems of interest in functional connectomics. We studied the activation maps related to a task that involves language (story understanding) and mathematics (mental computation). This particular task is expected to outline number, attentional and language networks, but the variability modes observed in the population cover even wider systems. For the experiments, mass-univariate General Linear Models (GLMs) [Friston et al., 1995] for n = 500 subjects were estimated for the Math vs Story contrast (language protocol), and the corresponding full-brain Z -score maps each containing p = 2.6 × 10 5 voxels, were used as the input data X ∈ R n×p , and we sought a decomposition into a dictionary of k = 40 atoms (components). The input data X were shu ed and then split into two groups of the same size.

Models compared and metrics. We compared our proposed Smooth-SODL model (9.2) against both the Canonical ICA -CanICA [Varoquaux et al., 2010], a single-batch multi-subject PCA/ICA-based method, and the standard SODL (sparse online dictionary-learning) [Mairal et al., 2010]. While the CanICA model accounts for subject-to-subject di erences, one of its major limitations is that it does not model spatial variability across subjects. Thus we estimated the CanICA components on smoothed data: isotropic FWHM of 6mm, a necessary preprocessing step for such methods. In contrast, we did no pre-smoothing for the SODL of Smooth-SODL models. The di erent models were compared across a variety of qualitative and quantitative metrics: visual quality of the dictionaries obtained, explained variance, stability of the dictionary atoms, their reproducibility, performance of the dictionaries in predicting behavioral scores (IQ, picture vocabulary, reading pro ciency, etc.) shipped with the HCP data [van Essen et al., 2012]. For both SODL [Mairal et al., 2010] and our proposed Smooth-SODL model, the constraint set for the dictionary atoms was taken to be a simplex C := S p (τ ) (see section 9.2 for de nition). The results of these experiments are presented in Fig. 9.3 and 9.5.

Results

alitative assessment of dictionaries. As can be seen in Fig. 9.3, all methods recover dictionary atoms that represent known functional brain organization; notably the dictionaries all contain the well-known executive control and attention networks, at least in part. Vanilla dictionary-learning leverages the denoising properties of the 1 sparsity constraint, but the voxel clusters are not very structured. For, example most blobs are surrounded with a thick ring of very small nonzero values. In contrast, our proposed regularization model leverages both sparse and structured dictionary atoms, that are more spatially structured and less noisy.

In contrast to both SODL and Smooth-SODL, CanICA [Varoquaux et al., 2010] is an ICA-based method which enforces no notion of sparsity whatsoever. The result are therefore dense and noisy dictionary atoms that explain the data very well (Fig. 9.4 but which are completely unintepretable. In a futile attempt to remedy the situation, in practice such PCA/ICA-based methods (including FSL's MELODIC tool [Smith et al., 2004]) are hardthresholded in order to see information. For CanICA, the hard-thresholded version has been named tCanICA in Fig. 9.3. That notwithstanding, notice how the major structures (parietal lobes, sulci, etc.) in each atom are reproducible across the di erent models.

Stability-fidelity trade-o s. PCA/ICA-based methods like CanICA [Varoquaux et al., 2010] and MELODIC [Smith et al., 2004] are the optimal linear decomposition method to maximize explained variance on a dataset. On the training set, CanICA [Varoquaux et al., 2010] out-performs all others algorithms with about 66% (resp. 50% for SODL [Mairal et al., 2010] and 58% for Smooth-SODL) of explained variance on the training set, and 60% (resp. 49% for SODL and 55% for Smooth-SODL) on left-out (test) data. See Fig. 9.4. However, as noted in the above paragraph, such methods lead to dictionaries that are hardly intepretable and thus the user must recourse to some kind of post-processing hard-thresholding step, which destroys the estimated model. More so, assessing the stability of the dictionaries, measured by mean correlation between corresponding atoms, across di erent splits of the data, CanICA [Varoquaux et al., 2010] scores a meager 0.1, whilst the hard-thresholded version tCanICA obtains 0.2, compared to 0.4 for Smooth-SODL and 0.1 for SODL. As is to be expected, notice how the RAW model over-ts. The voxel space of worth p = 261596 voxels is reduced to k = 40 components, and then each subject Z -map x t of worth p = 261596 voxels is reduced to k = 40 coe cients via a simple ridge regression (9.3).

Prediction of behavioral variables. If Smooth-SODL captures the patterns of inter-subject variability, then it should be possible to predict cogni-Figure 9.3: Qualitative comparison of the estimated dictionaries. Each column represents an atom of the estimated dictionary, where atoms from the di erent models (the rows of the plots) have been matched via a Hungarian algorithm. Here, we only show a limited number of the most "intepretable" atoms. Notice how the major structures in each atom are reproducible across the di erent models. Maps corresponding to hard-thresholded CanICA [Varoquaux et al., 2010] components have also been included, and have been called tCanICA. In contrast, the maps from the SODL [Mairal et al., 2010] and our proposed Smooth-SODL (9.2) have not been thresholded. Right: Predicting behavioral variables of the HCP [van Essen et al., 2012] dataset using subject-level Z -maps.

tive scores y like picture vocabulary, reading pro ciency, math aptitude, etc.

(the behavioral variables are explained in the HCP wiki [hcp]) by projecting new subjects' data into this learned low-dimensional space (via solving the ridge problem (9.3) for each sample x t ), without loss of performance compared with using the raw Z -values values X. 2) over the reference sparse online dictionary-learning (SODL) model [Mairal et al., 2010]. Note the reduction in the gain in EV as more data are pooled.

                 

Concluding remarks

To extract structured functionally discriminating patterns from massive brain data (i.e data-driven atlases), we have extended the online dictionary-learning framework rst developed in [Mairal et al., 2010], to learn structured regions representative of brain organization. To this end, we have successfully augmented [Mairal et al., 2010] with a Laplacian prior on the component maps, while conserving the low numerical complexity of the latter.

Through experiments, we have shown that the resulting model -Smooth-SODL model ( 9.2)-extracts structured and denoised dictionaries that are more intepretable and better capture inter-subject variability in small medium, and large-scale regimes alike, compared to state-of-the-art models. We believe such online multivariate online methods shall become the de facto way do dimensionality reduction and ROI extraction in future.

Implementation. The authors' implementation of the proposed SSODL (9.2) model will soon be made available as part of the Nilearn package [Abraham et al., 2014].

Possible extensions

More general structure-imposing penalties. One can envisage to replace the Laplacian regularization with a general structure-inducing penalty for which the proximal operator is easy to compute. Such a framework is developed in chapter 10, and produces an entire family of models, with potentially di erent properties. Replacing the dictionary with a general neural net. One notes that the proposed model (9.2) can be seen as an auto-encoding model of brain data with linear generator 9.17) parametrized by the shared learned dictionary D, and an implicit encoder 9.18) which is by construction, the optimal encoder 5 for the the generator G D . 5 Under the prior for the codes p(c) ∝ exp(-αϕ (c)).

G D = D, . : R k → R p , c → x := Dc, (
E D : R p → R k , x → arg min c∈R k -loglik(x|c, D) = arg min c∈R k (G D (c), x) + αϕ (c), (
One could obtain much greater modeling exibility by replacing the generator (9.17) by a neural network (see Fig. 9.6)

G θ : R k → R p , c → x := G θ (c),
with parameters θ ∈ Θ. The model could be successfully trained via stochastic gradient descent (SGD) or its various enhanced variants. Such models, Proximal updates for online dictionarylearning 9.2) proposed in chapter 9, by allowing for more general block-separable penalties on the dictionary.

The power of the prox

Recall from equation (9.2) of chapter 9 that after n passes over the data, the unpenalized objective (i.e loss) function whose minimization gives the dictionary updates is E (D) := 1 2n X -DC 2 F , where C ∈ R n×k is the xed matrix of codes computed up to this point, and D ∈ R p×k is the dictionary variable. Adding a block-separable penalty term γ k j=1 j (d j ), the energy becomes

E (D) = 1 2n X -DC 2 F + γ j j (d j ). (10.1)
In particular, taking j (d j ) := 1 2 ∇d j 2 F + i B p,1 (d j ) corresponds to the model proposed in chapter 9. Now, one easily computes

∇ D 1 2n X -DC 2 F = 1 n (DC -X)C T = DA -B,
where

A := 1 n CC T = 1 n n i=1 c i c T i ∈ R k×k and B := 1 n XC T = 1 n n i=1
x i c T i ∈ R p×k . Now in BCD, the jth atom is updated whilst all the others are held constant. Selecting the jth column of (10.1), we get ∇ d j 1 2n X -DC 2 F =

Applications

Now that we have the hammer, where are the nails...

Special cases

Constraint sets

If we take j := i C j , the indicator function of a closed convex subset of R p , so that each atom d j is constrained to satisfy a set of constraints prescribed by C j , then the above updates reduce to projecting z -j onto C j , namely .5) This is interesting as long as the C j 's are su ciently "simple" to allow us compute the above projection easily (preferably in closed form).

d j ← proj C j (z -j ). ( 10
Classical choice. Taking C j to be B 2 , the unit ball for the euclidean norm on R p , we recover the updates proposed in [START_REF] Mairal | Online dictionary learning for sparse coding[END_REF][Mairal et al., , 2010]], namely .6) One notes that these constraint has no structural properties beyond preventing the dictionary atoms from becoming arbitrarily large.

d j ← z -j max(1, z -j 2 ) . ( 10
Gram-Schmidt / step-wise orthonormality constraints. Akin to ICAtype methods, one can take C j = the orthogonal complement of the linear span of the rst j -1st atoms, namely

C j = span{d l |l < j} ⊥ ∩ B 2 . (10.7)
The dictionary updates are then simply the Gram-Schmidt orthonormalization of the ordered sequence of vectors z -1 , . . . , z -j , namely 2 2 The version of the Gram-Schmidt process presented here is not to be implemented as stated, as is it known to su er from numerical instability errors in nite-precision arithmetic. There exists equivalent versions (e.g Golub & Van Loan 1996] which alleviate these instabilities.

dj ← z -j - l <j proj d l (z -j ), d j ← dj dj 2 , ( 10.8)
where

proj d l (z -j ) :=        0, if d l = 0, z -j ,d l d l ,d l d l
, otherwise is the orthogonal projection of z -j onto the line generated by the atom d l . 10.2.2 "Social" sparsity: simultaneous sparsity and smoothness via windowed group-Lasso

The rst non-trivial results of our ramblings this far is obtained by considering the social sparsity [Kowalski et al., 2013, Kowalski and[START_REF] Kowalski | Structured sparsity: From mixed norms to structured shrinkage[END_REF] prior. In this model, a weakly activated voxel in the middle of strongly activated voxels will be saved (as if rescued by the clan), whilst a strongly activated voxel in the middle of weakly activated voxels will be eliminated (as if killed by isolation). "Strongness" and "weakness" are measured with respect to a speci ed threshold α > 0, which plays a rule similar to the regularization parameter in Group-Lasso. The penalty imposes both sparsity and structure simultaneously! Formally, social sparsity corresponds to a penalty social : R p → R de ned implicitly via its proximal operator ). Outside the red group, the weights are equal to zero. When considering the blue group, coe cients are weighted by some weights γ k 2 k 0, k ∈ N (k 2 ). Adapted from [Kowalski et al., 2013].

(prox α social (z)) := z        1 -α γ •z 2 , if γ • z 2 > α, 0, otherwise = z 1 - α s ∈N ( ) (γ s ) 2 z 2
One notes the following facts

• Ez 2 2 ≡ z 2 2 , where Ez := (γ • z) ∈[[p]] ∈ R p 2
is the expansion operator associated with the weights w s . In other words, E is a linear isometry.

• social sparsity is related to Group-Lasso (GL) by the formula

prox α social = F • prox α GL • E, (10.10)
where F is the right pseudo-inverse of E.

Conclusion

These ideas have a great potential to extend the classical dictionary-learning technology providing the practitioner with a modeling framework incorporating a much larger class of constraints -namely proximable penalty functions-than is currently being done. As regards convergence of our proposed proximal online dictionary-learning scheme, direct application of [Fercoq and Richtárik, 2015] seems to su ce, since the general DL algorithm constructs unbiased estimates of ∇ j f (D t ), where f (D) := E x min c∈R k (Dc, x). However, a more careful treatment is warranted, and is left for future work.

In neuroimaging and cognitive neuroscience, it is widely believed that the functional connectivity (FC) structure at rest remains grossly unchanged during task-stimulus presentation. This makes sense by least-action principle considerations: the brain does not need to rewire the functional links between regions upon presentation of a stimulus: it conserves the same networks as during rest, except that some edges are strengthened while others are weakened, to support the cognitive load of the particular task. Pushing this even further, one can claim that the resting-state FC of the brain predictively modulates the functional responses of the brain in the presence of task. Indeed, recently, resting-state fMRI has been shown to provide relevant constraints for functional mapping, opening the possibility to capture in standardized and cheaper acquisition most of the inter-individual differences [Tavor et al., 2016[START_REF] Cole | Activity ow over resting-state networks shapes cognitive task activations[END_REF][START_REF] Bzdok | Formal Models of the Network Cooccurrence Underlying Mental Operations[END_REF]. Possible applications include the improvement of population-level analyses, e.g. by nding better imputation schemes when dealing with missing data, detecting outlier data, and clarifying between-subjects similarities in comparison with genotyping or behavioral data. An important practical question has become how to optimize information transfer across these modalities to boost the chance of capturing the essence of inter-individual di erences.

Our main contributions. In this work, we propose a general framework for the problem of predicting task fMRI activation maps from resting-stateonly features. We present 2 main contributions: (a) the stacking of data across di erent random subsets subjects to reduce model-complexity and improve the prediction on held-out subjects, and (b) a multi-target regression approach to the predictive problem which better captures the functional inter-dependencies between di erent cognitive tasks. This generalizes and improves on the ideas in [Tavor et al., 2016]. We demonstrate the empirical gains brought by this approach through experiments on real datasets.

Feature extraction

The goal is to extract from resting-state data, pertinent features that encode the functional connectivity information in each voxel. A naïve choice would be to use the adjacency vector of each voxel in the whole-brain functional connectivity matrix. This is not practical due to the large number of (noisy) voxels, as it leads to enormous feature matrices. However, due to the inherent local correlations of data from di erent voxels, all this information is captured in the a nity of each voxel to a set of brain regions or networks. One way to get such pro les is to automatically learn a lowdimensional reduction of the resting-state data X s ∈ R n s ×p of each subject s into a common latent space, of dimension k min(min s (n s ), p), as proposed in [Tavor et al., 2016]. Here, n s is the number of TRs (Repetition Time) and p is the number of voxels.

Dual regression

As before, let the n s -by-p matrix X s be the resting-state data for subject s and D be the k-by-p group-level dictionary (aka topographic basis) obtained by stacking together resting-state time-series data from all the subjects and decomposing into k components of p voxels each by running a multi-subject decomposition algorithm like PCA, ICA, or dictionary-learning, etc. (more details on obtaining D later). We assume D to be under-complete -i.e k min(min s (n s ), p)-and therefore full-rank (i.e D DT is invertible). The standard "dual-regression" procedure [Tavor et al., 2016] then proceeds as follows:

• Compute the n s -by-k matrix of subject-to-components temporal dynamics C s by regressing the group-level dictionary D onto subject data X s :

Ĉs ∈ argmin C s ∈R ns ×k X s -C s D 2 Fro (11.1)
• Compute individual dictionary Ds = ( d s,j ) 1≤j ≤k ,1≤ ≤p ∈ R k×p by regressing the subject's resting-state data X s ∈ R n s ×p onto her subject-tocomponents temporal dynamics Ĉs ∈ R n s ×k :

Ds ∈ argmin D∈R k ×p X s -Ĉs D 2 Fro (11.2)
The end result is that for each subject s and each voxel , we obtain a kdimensional encoding ds, ∈ R k of the voxel's time-series x x , ∈ R n s in a common group-level space. These are the features (see Fig. 11.1).

Using only a single regression step

For the standard "dual-regression" feature-extraction method [Tavor et al., 2016], a total of 2 regression steps are done (hence the name of the procedure). As a rst (conceptual) improvement, we note that the individual dictionary Ds = Ĉ † s X s can be rewritten as

Ds = ( ĈT s Ĉs ) -1 ĈT s X s = (( D DT ) -1 DX T s X s DT ( D DT ) -1 ) -1 ( D DT ) -1 DX T s X s = D DT ( DX T s X s DT ) -1 DX T s X s = D DT DX T s (X s DT DX T s ) † X s OLS(X s DT ,X s )
.

That is, we regress the subject's resting-state time-series data X s onto n sby-k matrix Ĉs := X s DT and then reweight the result by the componentto-component covariance matrix D DT of the group-level dictionary. All in all, only a 1 regression step is needed.

Obtaining the global dictionary D

Since the resting-state time-series data are large (for example 1200 3D volumes of 2 × 10 5 voxels in each subject in the HCP -Human Connectome Project-dataset [van Essen et al., 2012]), a decomposition method that scales well is required. We use a variant of online dictionary-learning method [Mairal et al., 2010], a very fast implementation of which has been proposed in [Mensch et al., 2016], based on random matrix sketching / sub-sampling.

Incremental PCA/ICA-based methods [START_REF] Stephen | Group-PCA for very large fMRI datasets[END_REF], Varoquaux et al., 2010] are also a competitive choice.

Relationship between dual-regression and hyper-alignment

It turns out that the shared-response "hyper-alignment" (HA) framework [START_REF] Haxby | A common, high-dimensional model of the representational space in human ventral temporal cortex[END_REF] and the "dual regression" (DR) scheme [Tavor et al., 2016] we just presented are very closely related. Indeed, [START_REF] Haxby | A common, high-dimensional model of the representational space in human ventral temporal cortex[END_REF] considers the following problem minimize 11.3) Without the orthonormality constraints "C T s C s = I k ", this problem is precisely the DR problem. (11.3) is usually solved via an alternating minimization scheme. Viz,

1 N N s=1 X s -C s D 2 Fro over D ∈ R k ×p , C s ∈ R n s ×k , subject to C T s C s = I k , ∀s ∈ [1 . . . N ]. (
• Update rotations (orthonormal Procrustes analysis):

C (t +1) s = U (t ) s V (t ) s T ∀s ∈ [1 . . . N ],
where

U (k ) s Σ s V (t )
s T is the SVD of the n s -by-k matrix X s D (t ) T .

• Update shared-dictionary:

D (t +1) s = 1 N N s=1 D (t +1) s ,
where

D (t +1) s := C (t +1) s T X s .
However,

• DR is much more attractive due to its low cost: HA performs an SVD per subject per iteration.

• HA is usually done parcel-wise (i.e locally) because, the orthonormality conditions are unreasonable globally (i.e full-brain). We now develop our model for predicting subject-speci c activation maps Y s from resting-state features D s (refer to section 11.2). Our model considers bootstraps of sub-samples of subjects instead of on a subject-bysubject basis enforces a reduction in the complexity of the model without loss in capacity. The idea is that the regression coe cients from predicting task activations from resting state should be partly shared across subjects. This re ects the hypothesis that the global cognitive organization of the brain should share some similarities across di erent subjects. Also, stacking across subjects as such corrects for covariate-shift 1 between di erent 1 Grossly speaking, covariate-shift is a situation in which the distribution of the test set is not the same as the distribution of the training set, and so the model learned on the training set may not generalize to unseen (i.e test) data.

subjects, and facilitates transfer-learning from one-subject to another at test time.

Thus, for a bootstrap sub-sample

S of b (1 ≤ b ≤ N ) subjects, let Z S = [Z 1 | . . . |Z b ] ∈ R k ×p b
be functional features masked over a parcel P of p ≤ p voxels (see Fig. 11.2) and horizontally stacked matrix of functional features. Similarly, let Y S ∈ R p b×c be corresponding activation maps to the c task contrasts, masked over the same parcel, and stacked vertically. The goal is to link these functional features with activations maps corresponding to c task-activation contrasts e.g. "Story-vs-Math", "Faces-vs-Houses", "2Back-vs-0Back", etc.

Low-rank Ridge regression

Intra-subject activation maps for so-called di erent experimental stimuli may be correlated to one another. Indeed one would expect all brain activations to any conceivable experiment to be driven by a restricted set of latent causes, which is much less than the number of possible experiments. Thus, in an experiment with a su ciently large bail of experimental conditions, one would expect that corresponding action maps would be correlated across di erent experimental conditions. Fitting a separate model per experimental condition would therefore be statistically ine cient due to model over-speci cation. We need a principled way to incorporate the covariance structure of the intra-subject activation maps into our predictive model. Low-rank linear models do just this. It produces a much smaller model (i.e few number of free parameters) which best explains the covariance structure between activation maps Y for the di erent conditions. This can be written as

Find β S ∈ R k ×c , with rank(β S ) ≤ r , s.t Y j S ≈ Z T S β j S ∀j ∈ [1 . . . c], (11.4)
for a chosen rank bound r , with 1 ≤ r ≤ min(c, k ). Here, Y j S ∈ R |S |p×1 denotes the activation maps for contrast j, for the subjects in the bootstrap sub-population S. The model (11.4) (11.5) This de nes a low-complexity linear model 11.6) for predictively linking resting-state data to individual activation maps over the parcel P. The full-rank case r = min(c, k ) together with the choice b = 1 (no bagging) corresponds to the subject-wise contrast-wise single-output linear regression model proposed in [Tavor et al., 2016]. Now, let ŶOLS S = UΣV T be the SVD (singular-value decomposition) of the least-squares prediction ŶOLS S := Z S βOLS S where βOLS S := (Z T S Z S ) † Z T S Y S is the ordinary least-squares (OLS) solution to the unconstrained version of (11.5). Of course (11.5) may fail to have a unique solution. The following elementary lemma, whose proof (Supp. Mat.) follows directly from the Eckart-Young-Mirsky theorem [START_REF] Carl | The Approximation of One Matrix by Another of Lower Rank[END_REF] and the orthogonality property of the OLS t, produces a solution for model (11.5). Viz, Lemma 3. A solution to (11.5) is given by β S = βOLS S Π S (r ) where Π S (r ) = r i=1 v i v T i is the orthogonal projector onto the subspace spanned by the rst r principal singular vectors v i ≤r of the OLS prediction ŶOLS S .

fS : Z s → Z T s β S (
It should be noted that the form of the solutions provided by the above lemma is particularly appealing: We only need to do a single t to obtain a solution to problem (11.5) from solutions to the unconstrained OLS version of [Tavor et al., 2016].

Proof. Indeed, by the orthogonality property of least-squares estimates, we have the decomposition

Y S -Z S β S 2 Fro = Y S -ŶOLS S 2 Fro + ŶOLS S -Z S β S 2 Fro ,
with the rst summand being independent of the model parameters β S . Thus (11.5) can be rewritten as

minimize 1 2 ŶOLS S -Z S β S 2 Fro w.r.t β S ∈ R k×c subject to rank(β S ) ≤ r . (11.7) It is clear that β S (r ) := βOLS S Π S (r ) is of rank at most r . One computes Z S β (r ) = Z S βOLS S Π S (r ) = Z S βOLS S r i=1 v i v T i = ŶOLS S r i=1 v i v T i ,
which, by the Eckart-Young-Mirsky theorem [START_REF] Carl | The Approximation of One Matrix by Another of Lower Rank[END_REF] for the Frobenius norm, is the best rank r approximation of ŶOLS S w.r.t the Frobenius norm.

Algorithms 11.4.1 Learning

Low-rank multi-output regression. For the estimation of the predictive model linking resting-state features to activation maps, the template model (11.5) is solved for each parcel and each bootstrap sub-sample of subjects, to obtain the coe cients for predicting individual subject activations for the di erent task contrasts, jointly. The estimation is massively parallel: it is done per bootstrap and per parcel.

Hyper-parameter tuning

The rank bound r can be selected via K-fold cross-validation: we would retain the smallest value or r = rS in the range [1, min(k, c))] which produces a cross-validation score within 1 standard deviation of the best score (the so-called 1 standard error rule), or alternatively via leave-one-out (LOO) cross validation. However, cross-validation is very costly as multiple models must be tted on di erent splits of the training data. Moreover it might not even be possible in the case limited training data. • A set of brain parcellations (de ned by sets of brain masks).

Ensure: Distributed sets of tted models {{ fS | fS ∈ F P }|P ∈ parcels}, i.e one model fS per bootstrap sub-sample S per parcel P.

1: parallel for each parcel P do 2:

parallel for each bootstrap sub-sample of subjects S do 3:

Fit a model fS from (11.5), for predicting Y S from Z S restricted on the parcel P 4:

end pararell for 5: end pararell for Generalized cross-validation. A very attractive alternative to crossvalidation is the so-called generalized cross-validation (GCV) [START_REF] Golub | Generalized Cross-Validation as a Method for Choosing a Good Ridge Parameter[END_REF], whereby one attempts to directly minimize (an unbiased estimate of) the generalization error, which in our case reduces to

GCV (r ) := Y -Ŷ(r ) 2 Fro (nc -d f (r )) 2 (11.8)
as a function of the rank parameter r . Hee, d f (r ) is an unbiased estimate of the number of degrees of freedom, that is, the number of free parameters needed to completely specify the linear prediction model given by the coe cients β (r ). One can show that GCV is the consistent asymptotic limit of leave-one out (LOO) cross-validation. However the advantage of GCV is that only one model needs to be tted per value of the hyper-parameter (cf. LOO cross-validation, where as many models as sample points need to tted per value of the hyper-parameter).

In our case of reduced rank linear regression, the approximation error term Y -Ŷ(r ) 2

Fro in (11.8) reduces to

Y -Ŷ(r ) 2 Fro = Y -ŶOLS 2 Fro + ŶOLS -Ŷ(r ) 2 Fro = Y -ŶOLS 2 Fro + r 0 l =r +1 σ 2 l , (11.9)
where σ 1 ≥ σ 2 ≥ . . . ≥ σ r 0 are nonzero singular values of ŶOLS and r 0 = rank( ŶOLS ). In [START_REF] Mukherjee | On the degrees of freedom of reduced-rank estimators in multivariate regression[END_REF], nite-sample unbiased estimates for the degrees of freedom of rank-penalized models where derived. The authors established the formula (11.11) where κ ( ŶOLS ) := σ 1 /σ r 0 is the condition number of ŶOLS . Albeit, this bound is not very interesting for even mildly ill-conditioned ŶOLS where κ ( ŶOLS ) 2 1. Finally, we note that though d f naïve (r ) under-estimates d f (r ) in (11.10) and (11.11), the former is already a very good approximation in practice, and in fact equals the latter (almost surely) in the asymptotic limit n → ∞.

Inference

At prediction time, these di erent models are queried on held-out subjects and their results are aggregated by averaging. Such a divide-and-conquer approach allows us to learn complementary aspects of the data landscape, boosting prediction scores, while reducing the variance of the individual component models of which its is made. This is a well-known statistical property of bagging ensembles. The inference be done by making a single pass in Alg. 6. Algorithm 6: Predicting activation maps from resting-state features Require: • Data from N test subjects. For each subject s, we have precomputed spatial features Xs ∈ R k ×p using their resting-state data.

• Sets of tted models fS (see Alg. 5).

Ensure: Predictions Ŷs ∈ R p×c , for each test subject s. end pararell for 8: end pararell for

Experiments

Our experiments were done on task fMRI data from 200 subjects from the HCP -Human Connectome Project-dataset [van Essen et al., 2012]. These task fMRI data were acquired in an attempt to assess major domains that sample the diversity of neural systems , including language processing (semantic and phonological processing) and working memory.

The activation maps Y to predict. This data includes task activation maps from General Linear Models (GLMs) [Friston et al., 1994] that show the activation of di erent brain voxels to di erent cognitive conditions / task contrasts, for each subject. For example of these conditions include "Math vs Story" (part of language task), and "2Back-vs-0Back" -or "2BK-vs-0BK" for short-(part of working memory task). For example, there are about 19 task contrasts activation maps -each containing p = 2 × 10 5 voxels-per subject for the working memory protocol. For each subject s, this gives an output matrix Y s ∈ R p×c , were c is the total number of contrasts considered. In our experiments, we only considered the language (3 contrasts) and working memory (19 contrasts), giving a total of c = 22 task contrasts.

Extracted resting-state only features. The data also comes shipped with resting-state fMRI data consisting of n s = 1200 3D volumes of p = 2 × 10 5 voxels each, per subject, forming an n s -by-p matrix X s . The feature extraction described in section 11.2 was then applied to transform each X s into low-dimensional functional connectivity features Xs ∈ R k ×p , with k = 100.

The setup. N train = 100 subjects were used in Alg. 5 to t an ensemble of models (section 11.3). We used parcellations in which each parcel was worth about 4000 voxels, for a total of about 60 parcels. N test = 100 subjects were held out for evaluating the models predictions, computed via Alg. 6.

Results

antitative metrics. Fig. 11.3 shows confusion matrices (via Pearson correlation) of predicted against true activation maps. We see that a subject's predicted activation maps are consistently more similar to their true activation than to other subjects', re ected by the fact that the confusion matrices are strongly diagonal-dominant. This is even more true for our proposed method. The Fig. 11.4 shows box-plots of prediction R 2 -score and Pearson correlation for the 47 distinct contrast of the HCP task fMRI dataset [van Essen et al., 2012]. We see that both the reference method [Tavor et al., 2016] and our proposed method successfully predict the subjects' activation maps well above chance, with our method doing much better. Math" task contrast. The left plot corresponds to the reference method [Tavor et al., 2016] while the right one is for our proposed method. Higher diagonal values is better. The strong diagonal dorminance of these matrices reveals that the predicted maps of the subjects are more similar to their true maps than to other subjects. Results for the reference method [Tavor et al., 2016] are also shown. Bottom: Pearson correlation for the same prediction problem. alitative metrics. In Fig. 11.5, we display level-curves of the population mean (magenta) of activation maps for the "Story-vs-Math" and "2BK-vs-0BK" task contrasts, superimposed on the true activation maps of the subjects (the background image). The population mean activation map (magenta) is shown as a baseline (dummy predictor). We see that the contours for the predicted activation maps using our proposed method (green) faithfully follow topography of the true activation maps, indicating that the model successfully predicted the topography of the subjects' activation patterns for the contrasts. .5: Level-curves of the population mean (magenta), predicted activation maps using our proposed method (green) and the reference method [Tavor et al., 2016] (cyan) for di erent contrasts. Each column represents a di erent subject (here 3), while each row represents a task contrast (here 2): rst row is for "2BK-0BK" and second row is "Story-vs-Math".

Concluding remarks

We have proposed a general framework for the problem of predicting task fMRI activation maps from resting-state-only features. Our method creates an ensemble of parcel-wise low-rank multi-target linear models, over different random sub-populations of the training subjects to leverage the full richness of the data and jointly predict activation maps to di erent cognitive hypotheses (task contrasts). This is a major improvement over the state of the art [Tavor et al., 2016], as con rmed by extensive experiments on real data.

A practical implication of our results is that, for population studies, a large amount of information can be captured solely by a T1 image + restingstate fMRI: faster, cheaper scanning OR more control on data quality (imputation, outlier control). This explores new avenues for exploring the human brain via resting-state data, in patients and healthy subjects alike.

Possible extensions with general approximators. Our work and also the previous works [Tavor et al., 2016[START_REF] Cole | Activity ow over resting-state networks shapes cognitive task activations[END_REF][START_REF] Bzdok | Formal Models of the Network Cooccurrence Underlying Mental Operations[END_REF] has shown beyond doubt, that there exists a predictive mapping from restingstate fMRI data X to task activation maps Y. That is to say, the activation patterns in a person's brain during task are pre-determined by its background functional organization at rest. These works (including ours) have been limited to linear regression models, largely due to the simplicity of the latter. A priori, there is no reason why such a presumably complex relation should be accurately captured with a straight line, since there are probably many di erent layers of abstraction between functional connectivity patterns all the way through to activation patterns observed during task.

A simple extension would therefore be to replace the linear (Ridge) regression used to predict task activations from resting-state features, with a small multi-layer perceptron (MLP) 2 , a cascade of linear transformations L l , 2 An MLP can in principle approximate any "reasonable" function up to within arbitrary precision.

merged via simple non-linear functions σ l like recti er linear units (ReLU) or sigmoids

X σ 1 •L 1 -→ . . . σ H •L H → Y.
The intermediate representations extracted by such a model would be important in their own right. There would be enough data to t such a model since voxels are the samples in this prediction problem, and there 2 × 10 5 voxels per subject. Indeed, a preliminary implementation of this generalization appears to give even much better prediction results (not presented here) than the improvements presented here over the state-of-the-art. This excursion will be continued in a future work.

So ware. The code for the models presented in this chapter will be made publicly available online soon.

Part

Conclusion

Summary of main contributions

This thesis kicked-o with the goal of developing methods for modeling inter-subject functional variability, the aim being to enhance the estimation of functional connectomes -data-driven regions of interest, connectivity matrices, etc.-across populations of subjects. Below we summarize some of my major contributions.

Scientific contributions

The quest led to the study of, and proposal of methods for, structured (sparsity, smoothness, etc.) multi-variate models for brain encoding / decoding [Dohmatob et al., 2015b, Abraham et al., 2014, Eickenberg et al., 2015, Pellé et al., 2016]. Nonlinear registration of functional brain images also came up as a natural concern, and we contributed a method for direct registration of functional brain images [Dohmatob et al., 2016a] (submitted to Neuroimage journal). We also improved the current state-of-the-art in ROI extraction and dimensionality reduction by combining techniques from online learning and structured sparsity (like TV-L1) to propose a novel scalable dictionary-learning framework for obtaining decompositions of brain images, which are closer to known neuro-biological organization of the brain: networks made of spatially localized smooth components with sharp boundaries.

The ultimate indicator for having understood a phenomenon is being able to recreate it, at least approximately. Indeed, Feynman once said, "What I cannot (re)create, I do not understand!" By combing techniques in generative modeling and ensembles, we improved state-of-the-art methods for predicting task-based activation maps (at the individual level!) from restingstate fMRI data, with accuracy well above chance. This work is being prepared for journal submission.

Conclusion and Perspectives

123

Due to the intimate relationship between modeling and optimization, the bulk of this work was made possible by development of new or improvement of existing methods of optimization, with scalability and robustness at heart [Dohmatob et al., 2014, 2015a, Varoquaux et al., 2015[START_REF] Dohmatob | A simple algorithm for computing Nash-equilibria in incomplete information games[END_REF].

Finally, some of the work done in the thesis have cross-fertilized other collaborative papers like [START_REF] Rahim | Integrating Multimodal Priors in Predictive Models for the Functional Characterization of Alzheimer's Disease[END_REF], Thirion et al., 2014].

So ware contributions

While preparing this PhD project, I have made contributions to numerous open-source projects, including:

• Nilearn http://nilearn.github.io/index.html: Python package for leveraging machine learning algorithm in neuro-imaging. For example, the multi-variate models presented in chapter 3 are implemented as part of this package.

• Pypreprocess https://github.com/neurospin/pypreprocess: Python scripts scripts for preprocessing and QA of MRI data.

• Nistats https://github.com/nistats/nistats: Python package for statistical analyis (GLM, permutation tests, etc.) on MRI data.

Ongoing work and future directions

Unified view on structured models for brain data. We are preparing journal paper synthesizing all our contributions in the regarding structured models for brain decoding and segmentation presented in chapter 3. This will bring these methods to the doorsteps of the neuroscience practitioner.

Non-linear generative models for inter-subject brain data and prediction of task-fMRI activity from resting-state data. As concerns the modelling of inter-subject variability (chapters 9 and 11), most of the work done in this thesis can be cast in a more exible framework of generative encoder-less models (see Fig. 9.6, for example) 1 , with the space of 1 Since the encoding representation is gotten by simply minimizing a reconstruction loss between the generated and the true brain image.

parameters carefully constrained to ensure tractability and intepretability.

Synthèse en français

La thèse à démarrée avec l'objectif de développer des nouvelles méthodes pour la modélisation de la variabilité inter-sujet fonctionnelle, le but ultime étant l'amélioration de l'estimation de connectômes fonctionnelles sur des populations de sujets (chapitre 2). Cette quête à conduit à la proposition des méthodes de pénalisation structurée (parcimonie, variation totale, etc.) multi-variées pour l'encoding / decoding [Dohmatob et al., 2015b, Abraham et al., 2014, Eickenberg et al., 2015, Pellé et al., 2016]. Le récalage fonctionnel est aussi souvenu naturellement, est nous avons contribué une méthode pour le récalage directe des images fonctionnelles (EPI) vers un cerveau standard (template) [Dohmatob et al., 2016a] (soumit au Frontiers). Voir chapitres 3, 4, 5, 6, 8, et 7. Nous avons aussi amélioré l'état de l'art sur l'extraction de régions d'intérêt et la réduction de dimension en neuro-imagérie, par des méthodes combinant des techniques d'apprentissage en ligne et de parcimonie structurée (par exemple avec les pénalités TV-L1). Notre proposition est une nouvelle technique d'apprentissage de dictionnaire pour la décomposition d'images de cerveau plus conformes avec des a priori neurobiologique sur l'organisation fonctionnelle du cerveau: des réseaux spatialement localisés avec des contours bien délimités. Il s'agit d'un modèle génératif de base dimension, encodant succinctement la variabilité inter-sujet. Nous referons le lecteur aux chapitres 9 et 10.

Finalement, nous nous sommes intéressés à l'utilisation de techniques d'apprentissage supervisé pour expliquer la relation entre l'activité spontanée (activité au repos) et les enregistrements avec stimulations (activité évoquée dans des conditions précises). Nous avons proposé une méthode couplant un apprentissage non-supervisé (de type réduction de dimension par apprentissage de dictionnaire partagé) et des modèles prédictifs de faible rang pour exploiter les interdépendances entre les di érentes fonctions cognitives. Les expériences numériques réalisées (200 sujets du projet HCP [van Essen et al., 2012]) montrent que nous apportons une amélioration considérable à l'état de l'art. Voir chapitre 11.

Les travaux réalisés on donné lieu à des nombreuses publication à des conférences et journaux tels que NIPS, ICASSP, MICCAI, Frontiers in Neurosciences, etc. Une liste complète des publications peut être consulter ma page Google scholar https://scholar.google.fr/citations? user=FDWgJY8AAAAJ&hl=fr. En chi res • Nombre total de publications ≥ 15.

• h index ≥ 4.

• 110 index ≥ 3, dont • Parcimonie et régularisation spatiale: [Dohmatob et al., 2014], [Dohmatob et al., 2015b], [Abraham et al., 2014], [Eickenberg et al., 2015], [Pellé et al., 2016] • Récalage: [Dohmatob et al., 2016a] • Optimisation: [Dohmatob et al., 2015a], [Varoquaux et al., 2015], [START_REF] Dohmatob | A simple algorithm for computing Nash-equilibria in incomplete information games[END_REF] • Modelisation de variabilité fonctionnelle inter-sujet: [Dohmatob et al., 2016b] • Neurosciences: [START_REF] Rahim | Integrating Multimodal Priors in Predictive Models for the Functional Characterization of Alzheimer's Disease[END_REF], [Thirion et al., 2014] Il y a aussi des manuscrits en cours de préparation pour être publier dans des journaux:

• Vue globale sur la parcimonie et régularisation spatiale en neuroimagérie: "Structured penalties for brain decomposition and decoding: a uni ed view", pour Neuroimage

• "Inter-subject registration of functional images: do we need anatomical images ?" , pour Frontiers

• "Enhanced prediction of task-based activation maps from resting-state data", pour Neuroimage Logicielles contribuées. Pendant la préparation de la thèse, des nombreuses contributions dans de projets open-source on été réalisées. Pour en citer quelques unes:

• Nilearn http://nilearn.github.io/index.html: Librairie Python pour l'apprentissage statistique pour la neuro-imagerie. Par exemple les méthodes multi-variées présentées au chapitre 3 font partir des modules de cette librairie.

• Pypreprocess https://github.com/neurospin/pypreprocess: Des scripts Python pour le pré-traitement d'images IRMf.

  1 = j (∇ x w) 2 j + (∇ w) 2 j + (∇ z w) 2 j w SV Sparse Variation regularization ∇ ρ w 2,1 = j ρ 2 w j + (1 -ρ) 2 ∇w) ∇w 1,1 = j |(∇ x w) j | + |(∇ w) j | + |(∇ z w) . . . .. . . . . . . . . . . . . . . . . . . . . 11 1.2 Sketch of contributions . . . . . . . . . . . . . . 13

Figure 2 . 2 :

 22 Figure 2.2: Motion-correction and coregistration. The top plot shows an overlay of a subjects anatomy onto their mean functional image (the background). Here the red contours are well aligned with the background image, indicative of a successful co-registration.Typical things that can go wrong include: lesions (missing brain tissue), bad orientation headers in the images, non-brain tissue in the images (e.g skull), etc.The bottom plot show estimated motion parameters for the subjects-head motion.Here all movements are well below 0.5mm, which is generally considered as ne. The preprocessing and plots displayed here were done using Pypreprocess https://github.com/

Figure 2 . 3 :

 23 Figure 2.3: Tissue segmentation and normalization. Showing (top) outlines of a subject's anatomical / T1 image (foreground) projected onto an MNI template image (background) and also the tissue probability maps (TPMs), after registration to the latter. The contours should match the background image well. Typically impediments to correct registration include: lesions (missing brain tissue), corrupted image headers, non-brain tissue in anatomical image (i.e needs brain extraction), etc.

Figure 2 . 4 :

 24 Figure 2.4: Subject-level Activation maps for auditory versus visual and visual versus auditory conditions. Here, we show an axial via (z-coordinate) of the Zvalues corresponding the size of the e ect, in each voxel. Values range from -13 (light blue) to +13 (light red). Analysis was done using the Nistats open-source Python library https://github.com/ nistats/nistats.

Figure 2 . 5 :

 25 Functional connectivity patterns extracted from resting state data. The nodes are regions of the brain, and the thickness of the edges represent the relative strength average signal between the two corresponding regions.

Figure 2

 2 Figure 2.6: Inter-subject functional variability. Showing Z maps across di erent subjects for activation in Story vs Math language condition of the HCP -Human Connectome Project-dataset [van Essen et al., 2012]. The acrosssubject mean activation (top row, middle column) is also shown. Notice how the activations di er across subjects both in magnitude and spatial location.

Figure 3 . 1 :

 31 Figure 3.1: Decoding models mine patterns of activity to discriminate between cognitive states [Dehaeneet al., 1998]. Di erent activation patterns re ect di erent mental states. For example, those associated with di erent images viewed by the subject. In a training phase, the classi er will learn to discriminate between brain activity measured under di erent cognitive states. In the testing phase the generalization performance of the trained model is quanti ed by evaluating the classi er on the testing set and comparing the output of the classi er with the true labels associated with the stimuli. The prediction accuracy of the model is used as a measure of the quantity of information about the cognitive task shared by the voxels. Adapted from[START_REF] Pedregosa-Izquierdo | Feature extraction and supervised learning on fMRI: from practice to theory[END_REF] .

Figure 3

 3 Figure 3.3: p unit ball for various values of p. The kinks in the cases 0 ≤ p ≤ 1 impose sparsity.One notes however that, the cases 0 ≤ p < 1 lead to non-convex intractable optimization problems.

  Figure 3.4: A cartoon showing a sparse and blobby (step-wise constant / cartoon-like) brain map, as would be sought for by Total-Variation regularization (9.2).

  (a) 1 iteration of 5-Fold cross-validation

Figure 3 . 7 :

 37 Figure 3.7: Bar chart showing percentage classi cation on leftout, for one-vrs-one classi cation on the visual recognition dataset[Haxby et al., 2001]. We see that the highly structured maps produced by SpaceNet models (3.1) (e.g seeFig. 3.6) are not at the expense of model accuracy.

Figure 4 . 2 :

 42 Figure 4.2: Benchmarking solvers for TV-1 penalized models. Top: TV-1 penalized Logistic Regression on the visual recognition face-house discrimination task. Top Left: excess energy E(w t ) -E (w * ) as a function of time. Top Right: convergence time of the various solvers for di erent choice of regularization parameters. Broken lines correspond to a tolerance of 10 0 , whilst full-lines correspond to 10 -2 . The thick vertical line indicates the best model selected by cross-validation. Bottom: TV-1 penalized Least-Squares Regression. Bottom Left: on the visual recognition face-house discrimination task; Bottom Right: on the Mixed gambles dataset. The thick vertical line indicates the best model selected by cross-validation.

Figure 5

 5 Figure 5.1: Univariate featurescreening for the Graph-Net[Hebiri and van de Geer, 2011, Grosenick et al., 2013] problem (3.1) on di erent datasets. This gure shows spatial maps of X T j y, thresholded so that only voxels j with (from left to rightmost column)|X T j y| ≥ p 10% (|X T y|), |X T j y| ≥ p 20% (|X T y|), |X T j y| ≥ p 50% (|X T y|), and |X T j y| ≥ p 100% (|X T y|) (fullbrain) respectively, survive. The green contours enclose the elite voxels which are selected by the screening procedure at the respective threshold levels. (a): Mixed Gambles dataset[Jimura and Poldrack, 2012].Weights maps obtained for the GraphNet model (3.1) with these di erent screening-percentiles are shown in Figure5.4. (c): OASIS dataset[START_REF] Daniel S Marcus | Open Access Series of Imaging Studies (OASIS): cross-sectional MRI data in young, middle aged, nondemented, and demented older adults[END_REF] with VBM. See Figure5.2 for weights maps and age predictions obtained using these di erent screeningpercentiles.

Figure 5 .

 5 Figure 5.2: Predicting age from gray-matter concentration maps from the OASIS dataset [Marcus et al., 2007]. Top: Weights maps (solutions to problem (3.1)). Bottom-left: Mean Square Error (MSE) in age prediction, for di erent subjects of the validation set, for varying levels of the early-stopping tolerance ("es tol" for short), with the screeningpercentile (sp) held constant at 100 (full-brain). Bottom-right: MSE in age prediction, for varying levels of the screening-percentile (sp).

Figure 5 .

 5 Figure 5.3: Predicting age from gray-matter concentration maps from the OASIS dataset[START_REF] Daniel S Marcus | Open Access Series of Imaging Studies (OASIS): cross-sectional MRI data in young, middle aged, nondemented, and demented older adults[END_REF]. Top: Weights maps (solutions to problem (3.1)). Bottom-left: Mean Square Error (MSE) in age prediction, for di erent subjects of the validation set, for varying levels of the early-stopping tolerance ("es tol" for short), with the screeningpercentile (sp) held constant at 100 (full-brain). Bottom-right: MSE in age prediction, for varying levels of the screening-percentile (sp). Running times: Increasing est tol (from -10 -4 to 10): 100.2m, 171.4m, 188.8m, 289.6m. For increasing sp (10 to 100): 44.2m, 81.3m, 186.5m, 341.3m Figure 5.4: Visual recognition dataset [Haxby et al., 2001]. (a): Weights maps for the Face vs House contrast, for di erent early-stopping and univariate feature-screening thresholds. One can see that the supports of these maps for di erent values of the thresholds are quite similar to cases involving no heuristic at all (the case where est = 10 and the where case sp = 100%). (b), top-left: Prediction scores as a function of the early-stopping tolerance (est), for di erent task contrasts. (b), top-right: Prediction scores as a function of the screening-percentile (sp), for di erent task contrasts. (b), bottom-row: Running times in minutes for the di erent thresholds of the heuristics.

  Contents 6.1 Derivation . . . . . . . . . . . . . . . . . . . . . . . 54 6.2 The algorithm: iGraphNet . . . . . . . . . . . . 55 6.3 Experimental results . . . . . . . . . . . . . . . . 56

Figure 6 . 1 :

 61 Figure 6.1: Estimated coe cients ŵ on the Face vs House condition of the visual recognition dataset[Haxby et al., 2001]. Classi cation accuracies on held-out data are shown in the legends. We monitor the evolution of the model as a function of the number of iGraph-Net iterations k = 0, 1, 2, . . .. We see that as more and more iterations of iGraphNet are run, the coe cients become more and more spatially denoised and localized (and therefore more intepretable), without deteriorating of the model accuracy.

( b )

 b Group-Lasso: Here K = I and φ = • 2,1 , n = 48, p = 128, number of blocks d = 32, block size = c = 4, q = d × c = 128, w 0 is has 2 non-zero blocks. (c) Sparse spikes deconvolution: Here, K = I, X is a projector onto low Fourier frequencies (Dirichlet kernel) and the penalty φ is the 1 -norm; n = p = 200 (with rank X = 40). The true signal w 0 is a 20-sparse vector (of length p), containing randomly distributed spikes with Gaussian values at a minimum pairwise distance of 5.

Figure 7 . 2 :

 72 Figure 7.2: Experimental results from ICASSP paper[Dohmatob et al., 2015]. showing local Q-linear convergence for ADMM on problem(7.1). The "theoretical" line is the exponential curve t → x (0)x * r (Λ (x * )) t . The red broken vertical line marks the instant the support of the xed-point x * is identi ed.

Figure 8 . 1 :

 81 Figure 8.1: Nonlinear mismatch between EPI and T1-weighted image of the same subject of the HCP dataset [van Essen et al., 2012], before and after distortioncorrection.Left: Single-band high-resolution EPI (SBRef) image of the same subject. Notice the large distortions along the Left-Right direction (inside the highlighted patches).Center: Distortion-corrected single-band EPI image. Here, the distortioncorrection managed to undo most -but not all-of the distortions. Even after distortion correction, there are minor shape (nonlinear) di erences between the EPI and the T1-weighted image of the subject (Right). The same nativespace coordinates where used in all of the 3 plots.

Figure

  Figure 8.4:Coe cient-of-Variation (CoV) after registration. Top: Log10 of ratio of across-subject Coe cient of Variation (CoV) for indirect T1-based pipeline / direct EPI-based. We see that the gain of our proposed method is most pronounced along the cortical surface.Bottom: Histograms of CoV for both pipelines. Again, we see clearly that our proposed method reduces the inter-subject variability by a much larger margin, indicative of improved subject-to-subject alignment.

Figure 8 . 5 :

 85 Figure 8.5: Mean EPI image across all subjects after registration (aka estimated population templates). Patches on the images have been zoomed to highlight details. The mean image from the indirect T1-based pipeline (Left) is more blurry (as seen here in the cerebellum), compared to our direct EPI-based pipeline post-registration across-subject mean image (Right) which is much sharper, indicative of a better inter-subject registration.Also the mean image from the T1-based pipeline has ripples on the cortical surface indicative of residual registration problems, which can be attributed residual EPI-distortions that could not be captured via coregistration.

Figure 8 . 6 :

 86 Figure 8.6: Qualitative comparison of pipelines via GLM results. Across-subject mean activation maps of Z -scores for di erent contrasts.Here we see that our proposed direct EPI-based registration scheme leads to slightly higher across-subject mean activation peaks. For each contrast, a cut has been made around the location of the activation peak, to display curves of the activation prole and across-subject variability thereof, in a neighborhood of this peak location.

Figure 8 . 7 :

 87 Figure 8.7: Comparing functional brain networks from subject fMRI images registered with both pipelines, namely the classical indirect T1-based method, and our proposed direct EPI-based method). Shown here are group-level unthresholded sub-component maps of the Default Mode Network (DMN)[START_REF] Raichle | A default mode of brain function[END_REF], using MNI coordinates reported in Table1of[START_REF] Watanabe | A pairwise maximum entropy model accurately describes resting-state human brain networks[END_REF].

Figure 9 . 2 :

 92 Figure 9.2: In uence of model parameters. In the experiments, α was chosen according to (9.9).

Figure 9 . 4 :

 94 Figure 9.4: Left: Mean explained variance of the di erent models on both training data and test (left-out) data. N.B.: Bold bars represent performance on test set while faint bars in the background represent performance on train set.

Figure 9 . 6 :

 96 Figure 9.6: General neural net schema for structured online dictionary learning. The model is trained via stochastic gradient descent -SGD. For an incoming image x ∈ R p (or mini-batch of images), a code z ∈ R k is sampled and fed into th generator G θ to produce a reconstructed image x recon ∈ R p , which is compared with the original via a loss function. The model(9.2) corresponds to a linear generator G θ = G D = D, . , and can be solved via Alg. 4.

Contents 10 . 1

 101 The power of the prox . . . . . . . . . . . . . . 103 10.2 Applications . . . . . . . . . . . . . . . . . . . . . . 105 10.2.1 Special cases . . . . . . . . . . . . . . . . . 105 10.2.2 "Social" sparsity: simultaneous sparsity and smoothness via windowed group-Lasso . . . . . . . . . . . . . . . . . 105 10.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . 106 I , we consider possible generalization of model (

  γ • z := (γ 1 z 1 , γ 2 z 2 , . . . , γ p z p ) ∈ R p for weights (γ s ) ,s ∈[[p]] satisfying |γ s | 2 = 1 for all s, and N ( ) := {s ∈ [[p]]|γ s 0} is the neighborhood of the th voxel, assumed to be non-empty. Thus each γ can be thought of as (normalized) mean-lter supported on a patch N ( ) around the voxel . Examples include rectangular lters, truncated Gaussians, etc.

Figure 10 . 1 :

 101 Figure 10.1: Social sparsity illustrated in 2D. The neighborhood of the coe cient k 1 is given by the red window, and the neighborhood of the coe cient k 2 by the blue one. These two neighborhoods share one coe cient. When considering the red group, coe cients are weighted by some weights γ k 1 k 0, k ∈ N (k 1 ). Outside the red group, the weights are equal to zero. When considering the blue group, coe cients are weighted by some weights γ k 2

Figure 11 . 1 :

 111 Figure 11.1: FC feature-extraction. Left: A component of the group dictionary. Right: Corresponding component for an individual subject's dictionary estimated using the proposed formula (11.2.2).

11. 3 Figure 11

 311 Figure 11.2: A parcellation is simply a collection of contiguous / simply-connected masks P called parcels (the colored patches) which cover the brain. Each voxel of the brain belongs / is assigned to exactly one parcel. In the parcellation shown here, each parcel contains approximately 1000 voxels.

Algorithm 5 :

 5 Training model for predicting activation maps from restingstate Require: • Data from N train subjects. For each subject s we have precomputed spatial features Z s ∈ R k ×p .

d

  f (r ) = d f naïve (r ) any r ∈ [1, r 0 ] with σ r > σ r +1 if r < r 0 . The naïve estimate d f naïve (r ) := kc -(kr )(cr ) is simply the number of free parameters needed to com-pletely specify a matrix of rank r ∈ [1, r 0 ]. trivial bound d f (r ) ≥ d f naïve (r ) + 2r (r 0 -r ) κ ( ŶOLS ) 2 ,

1 : 5 :

 15 Ŷ ← 0 ∈ R N test ×p×c 2: parallel for each parcel P do 3: parallel for each trained model fS on P do 4: parallel for each test subject s do Predict the activation maps of subject s with model fS : ŶS | P ← ŶS | P + fS ( Xs | P )

Figure 11 . 3 :

 113 Figure 11.3: Confusion matrices for predicted versus true activation maps for the "Story vsMath" task contrast. The left plot corresponds to the reference method[Tavor et al., 2016] while the right one is for our proposed method. Higher diagonal values is better. The strong diagonal dorminance of these matrices reveals that the predicted maps of the subjects are more similar to their true maps than to other subjects.

Figure 11

 11 Figure 11.4: Top: R 2 -score for predicting subject-speci c activation maps for di erent task contrasts, from their resting data. These results are for the di erent contrasts of the HCP dataset [van Essen et al., 2012] are shown.Results for the reference method[Tavor et al., 2016] are also shown. Bottom: Pearson correlation for the same prediction problem.

Figure 11

 11 Figure 11.5: Level-curves of the population mean (magenta), predicted activation maps using our proposed method (green) and the reference method[Tavor et al., 2016] (cyan) for di erent contrasts. Each column represents a di erent subject (here 3), while each row represents a task contrast (here 2): rst row is for "2BK-0BK" and second row is "Story-vs-Math".

  Summary of main contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12.1.1 Scientific contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12.1.2 So ware contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12.2 Ongoing work and future directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	Conclusion
	12 Concluding remarks
	12.1 13 Synthèse en français

1 Introduction to brain decoding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.2 Sparsity and structure-inducing priors: towards intepretable multi-variate models . . .

3.3 SpaceNet: sparse structured models for brain data . . . . . . . . . . . . . . . . . . . . . . .

  Contents 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . 58 7.1.1 The ADMM algorithms . . . . . . . . . 59 7.1.2 Examples . . . . . . . . . . . . . . . . . . . 59 7.2 Our contributions . . . . . . . . . . . . . . . . . . 60 7.2.1 Preliminaries . . . . . . . . . . . . . . . . . 60 7.2.2 Behavior of ADMM around fixed-points 61 7.3 Relation to prior work . . . . . . . . . . . . . . 62

	7.3.1 Ridge, QP, and nonnegative Lasso . . 62
	7.3.2 Fréchet-di erentiable nonlinear sys-
	tems

. . . . . . . . . . . . . . . . . . . . . . 63 7.3.3 Partly-smooth functions and Friedrichs angles . . . . . . . . . . . . . . . . . . . . . 63 7.4 Numerical experiments and results . . . . . 64 7.5 Concluding remarks . . . . . . . . . . . . . . . . 65

  Introduction . . . . . . . . . . . . . . . . . . . . . . 69 8.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . 71 8.2.1 An important note on normalization 71 8.2.2 General preprocessing procedures . . 71 8.2.3 The pipelines . . . . . . . . . . . . . . . . 71 8.3 Relation to previous works . . . . . . . . . . . 74 8.3.1 Direct EPI-to-EPI non-linear inter-subject registration . . . . . . . . . . . . . . . . . . 74 8.3.2 Non-linear EPI-to-structural coregistration

	Direct EPI-to-EPI inter-subject non-
	linear registration
	Contents
	8.1	Part
	III -Functional
	inter-subject variability
	Bibliography
	Luca Baldassarre, Janaina Mourao-Miranda, and Massimiliano Pontil. Struc-
	tured sparsity models for brain decoding from fMRI data. In PRNI, 2012.
	Heinz H Bauschke, JY Cruz, Tran TA Nghia, Hung M Phan, and Xianfu
	Wang. Optimal rates of convergence of matrices with applications.

arXiv:1407.0671, 2014. Ilker Bayram and Ivan W Selesnick. A Subband Adaptive Iterative Shrinkage/Thresholding Algorithm. Signal Processing, IEEE Transactions on, 58, 2010. . . . . . . . . . . . . . . . . . . . . . 75 8.4 Experiments . . . . . . . . . . . . . . . . . . . . . . 75 8.4.1 Evaluation metrics . . . . . . . . . . . . . 76 8.4.2 How many (plausible) pipelines are there ? . . . . . . . . . . . . . . . . . . . . . 76 8.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . 77 8.6 Discussion and concluding remarks . . . . . 80

  Summary of main contributions . . . . . . . 122 12.1.1 Scientific contributions . . . . . . . . . . 122 12.1.2 So ware contributions . . . . . . . . . . 123 12.2 Ongoing work and future directions . . . . . 123

	Concluding remarks
	Contents
	12.1

Acknowledgements

Many thanks to my PhD advisors Bertrand Thirion and Gael Varoquaux for their help and support during these 3 years preparing my PhD project. Special thanks to the Parietal team assistants Regine Bricquet and Ti any Caristan for their generous help with administrative matters. I would also want to thank the rest of the Parietal team -past and present-for ensuring such a nice work environment. Finally, I am thankful to all the members of the jury: John Ashburner, Gabriel Peyré, Marc Schoenauer, and Moritz Grosse-Wentrup, for kindly accepting to evaluate my PhD project.

Contents

I -Introduction

Fréchet-di erentiable nonlinear systems

In the SISTA algorithm [Bayram and Selesnick, 2010], the authors linked the rate of convergence of their multi-band ISTA (refer to [Daubechies et al., 2004] and the references therein, for the original ISTA algorithm) scheme to the spectral radius of a certain Jacobian matrix related to the problem data and dependent on the xed-point [Bayram and Selesnick, 2010, Propositions 6 and 7], provided this spectral radius is less than 1. Most importantly, the authors show [Bayram and Selesnick, 2010, Proposition 8] how their algorithm can be made as fast as possible by choosing the shrinkage parameter per sub-band to be "as large as possible". Finally, analogous to our Theorem 1(a), Lemma 2 of [Bayram and Selesnick, 2010] shows that the SISTA iteration projects points su ciently close to xed-points onto the support of these xed-points.

Partly-smooth functions and Friedrichs angles

In the recent work [START_REF] Liang | Activity Identication and Local Linear Convergence of Douglas-Rachford/ADMM under Partial Smoothness[END_REF] which focuses on Douglas-Rachford/ADMM, and [START_REF] Liang | Activity Identi cation and Local Linear Convergence of Inertial Forward-Backward Splitting[END_REF] which uses the same ideas as in [START_REF] Liang | Activity Identication and Local Linear Convergence of Douglas-Rachford/ADMM under Partial Smoothness[END_REF] but with a forward-backward scheme [START_REF] Patrick | Signal recovery by proximal forward-backward splitting[END_REF], the authors consider a subclass PSS (refer to de nition 2.2 of [START_REF] Liang | Activity Identi cation and Local Linear Convergence of Inertial Forward-Backward Splitting[END_REF]) of the class of so-called partly-smooth (PS) penalties and general C 2 loss functions with Lipschitz gradient. Under nonlinear complementarity requirements analogous to the non-degeneracy assumption "ϵ (x * ) > 0" of Theorem 1(b), and rank constraints analogous to the requirement that the Jacobian matrix Λ (x * ) have spectral radius less than 1 (in Theorem 1(c2)), the authors of [START_REF] Liang | Activity Identication and Local Linear Convergence of Douglas-Rachford/ADMM under Partial Smoothness[END_REF][START_REF] Liang | Activity Identi cation and Local Linear Convergence of Inertial Forward-Backward Splitting[END_REF] prove nite-time activity identi cation and local Q-linear convergence at a rate given in terms of Friedrichs angles, via direct application of [Bauschke et al., 2014, Theorem 3.10]. The authors show that their arguments are valid for a broad variety of problems, for example the anisotropic TV penalty. Still in the framework of partly-smooth penalties, [START_REF] Demanet | Eventual linear convergence of the Douglas-Rachford iteration for basis pursuit[END_REF] showed local Q-linear convergence of the Douglas-Rachford algorithm on the Basis Pursuit problem. [START_REF] Liang | Activity Identication and Local Linear Convergence of Douglas-Rachford/ADMM under Partial Smoothness[END_REF][START_REF] Liang | Activity Identi cation and Local Linear Convergence of Inertial Forward-Backward Splitting[END_REF]. The works which are most comparable to ours are [START_REF] Liang | Activity Identication and Local Linear Convergence of Douglas-Rachford/ADMM under Partial Smoothness[END_REF] and [START_REF] Liang | Activity Identi cation and Local Linear Convergence of Inertial Forward-Backward Splitting[END_REF], already presented above. Let us point out some similarities and di erences between these papers and ours. First, though our constructions are entirely di erent from the techniques developed in [START_REF] Liang | Activity Identication and Local Linear Convergence of Douglas-Rachford/ADMM under Partial Smoothness[END_REF][START_REF] Liang | Activity Identi cation and Local Linear Convergence of Inertial Forward-Backward Splitting[END_REF], one notes that both approaches are ultimately rooted in the same idea, namely the work of B. Holmes [START_REF] Holmes | Smoothness of certain metric projections of Hilbert space[END_REF] on the smoothness of the euclidean projection onto convex sets, and other related functionals (Minkowski gauges, etc.). Indeed, Theorem 1 builds directly upon [START_REF] Holmes | Smoothness of certain metric projections of Hilbert space[END_REF], whilst, [START_REF] Liang | Activity Identi cation and Local Linear Convergence of Inertial Forward-Backward Splitting[END_REF] and [START_REF] Liang | Activity Identication and Local Linear Convergence of Douglas-Rachford/ADMM under Partial Smoothness[END_REF] are linked to [START_REF] Holmes | Smoothness of certain metric projections of Hilbert space[END_REF] via [START_REF] Wright | Identi able surfaces in constrained optimization[END_REF], which builds on [START_REF] Fitzpatrick | Di erentiability of the metric projection in Hilbert space[END_REF], and the latter builds on [START_REF] Holmes | Smoothness of certain metric projections of Hilbert space[END_REF].

Detailed comparison with

Second, part (c1) of Theorem 1 ( nite-time identi cation of active set) of the theorem can be recovered as a consequence of the results established in [START_REF] Liang | Activity Identication and Local Linear Convergence of Douglas-Rachford/ADMM under Partial Smoothness[END_REF][START_REF] Liang | Activity Identi cation and Local Linear Convergence of Inertial Forward-Backward Splitting[END_REF]. However, the rest of our results, notably part (c2) (Q-linear convergence) cannot be recovered from the aforementioned      loadings (C) prediction IQ, age, etc.

Let RAW refer to the direct prediction of targets y from X, using the top 2000 most voxels most correlated with the target variable. Results of the comparison are shown in Fig. 9.4. Only variables predicted with a positive mean (across the di erent methods and across subjects) R-score are reported. We see that the RAW model, as expected over-ts drastically, scoring an R 2 of 0.3 on training data and only 0.14 on test data. Overall, for this metric CanICA performs best than all the other models in predicting the di erent behavioral variables on test data. However, our proposed Smooth-SODL model outperforms both SODL [Mairal et al., 2010] and tCanICA, the thresholded version of CanICA.

Running time. On the computational side, the vanilla dictionary-learning SODL algorithm [Mairal et al., 2010] with a batch size of η = 20 took about 110s (≈ 1.7 minutes) to run, whilst with the same batch size, our proposed Smooth-SODL model (9.2) implemented in Alg. 3 took 340s (≈ 5.6 minutes), which is slightly less than 3 times slower than SODL. Finally, CanICA [Varoquaux et al., 2010] for this experiment took 530s (≈ 8.8 minutes) to run, which is about 5 times slower than the SODL model and 1.6 times slower than our proposed Smooth-SODL (9.2) model. All experiments were run on a single CPU of a modern laptop.

Is spatial regularization really needed ? Ideally, one does not need spatial regularization if data are abundant (like in the HCP). So we computed learning curves of mean explained variance (EV) on test data, as a function of the amount training data seen by both Smooth-SODL and SODL [Mairal et al., 2010] (Fig. 9.5). In the beginning of the curve, our proposed spatially regularized Smooth-SODL model starts o with more than 31% explained variance (computed on 241 subjects), after having pooled only 17 subjects. In contrast, the vanilla SODL model [Mairal et al., 2010] scores a meager 2% explained variance; this corresponds to a 14-fold gain of Smooth-SODL over SODL. As more and more that are pooled, both models explain more variance, and the gap between Smooth-SODL and SODL reduces, and both models perform comparably asymptotically.

Da jb j . Putting things together, we get p = arg min

where z -j := a -1 j,j b jl j a j,l d l , and the last equivalence results from the following elementary lemma which reveals that the prox of a function at a point can be seen as an implicit gradient step. Viz

Lemma 2. For a function f : H → (-∞, +∞] (convex or not), recall the de nition of its subdi ferential at a point p ∈ H , namely

We have the following characterization of the prox

Putting everything together, we have the close-form BCD update formula for regularized online dictionary learning 1

1 We assumed W.L.O.G that a j,j 0 (i.e a j,j > 0), meaning that the jth atom is active in the representation of at least one sample. Otherwise, we can update the jth atom with a random vector, or even skip it altogether.

BCD updates DL with general separable penalties. The BCD updates for a penalized DL model (10.1) is d j ← prox γ a -1 j,j j (z -j ), (10.3) where z -j := a -1 j,j r j and r j := b jl j a j,l d l .

For the purposes of practical implementation, one notes that r j is precisely the jth column of the p-by-k matrix R := DA -B + d j • a j . (10.4) Thus at the beging of the BCD updates, we precompute the di erence DA -B, and substract (resp. add) the rank-1 term d j • a j before and after updating the j atom d according to (10.3). Such rank-1 updates are natively optimized in all linear algebra scienti c computation packages.