.

(2.7)

Titre : Analyse non-supervisée en composante de données d'imageries neuronales Mots clés : neuro-imagerie, analyse en composantes, apprentissage automatique Résumé : Cette thèse d'informatique et de mathématiques s'applique au domaine des neurosciences, et plus particulièrement aux recherches sur la modélisation de l'activité cérébrale humaine par électrophysiologie et imagerie. Dans ce champ, la tendance est actuellement d'expérimenter avec des stimuli naturels, comme le visionnage d'un film ou l'écoute d'une piste audio, et non plus avec des stimuli étroitement contrôlés mais outrageusement simples. L'analyse de ces stimuli « naturels » et de leurs effets demande toutefois de disposer d'une immense quantité d'images, par ailleurs très coûteuses. Sans outils mathématique, identifier l'activité neuronale à partir des données est quasi impossible. Toutefois, ces stimuli sont compliqués à modéliser et à analyser, car l'utilisation de méthodes fondées sur des régressions est limitée par la difficulté de modéliser les stimuli. C'est ce qui motive l'utilisation de méthodes non-supervisées qui ne font pas d'hypothèses sur ce qui déclenche les activations neuronales.

Dans cette thèse, nous considérons d'abord le cas du modèle de réponse partagée (MRP), dans lequel les sujets sont supposés partager une réponse commune. Ce modèle est utile pour réduire la dimension des données, mais son entraînement est coûteux pour les données d'imagerie fonctionnelle (IRMf) dont la dimension peut être immense. Nous présentons une version bien plus rapide et beaucoup plus économe en mémoire. Mais le MRP fait des hypothèses irréalistes sur les données d'imagerie.

Des hypothèses plus réalistes sont utilisées dans l'analyse en composantes indépendantes (ACI) mais cette méthode est difficile à généraliser aux jeux de données qui contiennent plusieurs sujets. Nous proposons alors une extension de l'ACI appelée ACI multi-vue, fondée sur le principe de maximum de vraisemblance et qui convient à des jeux de données multi-sujets. L'ACI multi-vue a une vraisemblance en forme fermée qui peut être maximisée efficacement. Toutefois, cette méthode suppose la même quantité de bruit pour tous les sujets.

Nous présentons donc l'ACI partagée, une généralisation de l'ACI multi-vue qui s'accompagne d'un modèle de bruit plus général. Contrairement à presque tous les modèles fondés sur l'ACI, l'ACI partagée peut séparer des sources gaussiennes et non gaussiennes et propose une estimation optimale des sources communes (au sens des moindres carrés), qui pondère chaque sujet en fonction de son niveau de bruit estimé. En pratique, l'ACI partagée et l'ACI multi-vue permettent d'obtenir, en magnéto-encéphalographie et en IRMf, une estimation plus fiable de la réponse commune que leurs concurrents.

Enfin, nous utilisons l'ACI comme base pour faire de l'augmentation de données. Plus précisément, nous présentons l'ACI conditionnelle, une méthode d'augmentation de données qui exploite la grande quantité de données d'IRMf non étiquetées pour construire un modèle génératif en utilisant seulement un petit nombre de données étiquetées. L'ACI conditionnelle permet d'augmenter de façon appréciable la précision du décodage sur huit grands jeux de données d'IRMf.

Nos principaux apports nous semblent consister dans l'accélération de l'entraînement du MRP ainsi que dans l'introduction de deux modèles plus réalistes pour l'analyse de l'activité cérébrale de sujets exposés à des stimuli naturels : l'ACI multi-vue et l'ACI partagée. Enfin, nos résultats sont prometteurs concernant l'utilisation de l'ACI pour faire de l'augmentation de données.

Nous présentons pour finir quelques pistes qui pourraient guider des travaux ultérieurs. D'un point de vue pratique, des modifications mineures de nos méthodes pourraient permettre l'analyse des données d'imagerie obtenues sur des sujets au repos en faisant l'hypothèse d'une organisation spatiale partagée à la place d'une réponse partagée. D'un point de vue théorique, les travaux futurs pourraient se concentrer sur la compréhension de la façon dont la réduction de dimensions et l'identification de la réponse partagée peuvent être réalisées conjointement.

Title : Unsupervised component analysis for neuroimaging data Keywords : neuro-imaging, component analysis, machine learning

Abstract :

This thesis in computer science and mathematics is applied to the field of neuroscience, and more particularly to the mapping of brain activity based on imaging electrophysiology. In this field, a rising trend is to experiment with naturalistic stimuli such as movie watching or audio track listening, rather than tightly controlled but outrageously simple stimuli. However, the analysis of these "naturalistic" stimuli and their effects requires a huge amount of images that remain hard and costly to acquire. Without mathematical modeling, the identification of neural signal from the measurements is very hard if not impossible. However, the stimulations that elicit neural activity are challenging to model in this context, and therefore, the statistical analysis of the data using regressionbased approaches is difficult. This has motivated the use of unsupervised learning methods that do not make assumptions about what triggers brain activations in the presented stimuli.

In this thesis, we first consider the case of the shared response model (SRM), where subjects are assumed to share a common response. While this algorithm is useful to perform dimension reduction, it is particularly costly on functional magnetic resonance imaging (fMRI) data where the dimension can be very large. We considerably speed up the algorithm and reduce its memory usage. However, SRM relies on assumptions that are not biologically plausible.

In contrast, independent component analysis (ICA) is more realistic but not suited to multisubject datasets. In this thesis, we present a wellprincipled method called MultiViewICA that extends ICA to datasets containing multiple subjects. MultiViewICA is a maximum likelihood estimator. It comes with a closed-form likelihood that can be efficiently optimized. However, it assumes the same amount of noise for all subjects.

We therefore introduce ShICA, a generalization of MultiViewICA that comes with a more general noise model. In contrast to almost all ICA-based models, ShICA can separate Gaussian and non-Gaussian sources and comes with a minimum mean square error estimate of the common sources that weights each subject according to its estimated noise level. In practice, MultiViewICA and ShICA yield on magnetoencephalography and functional magnetic resonance imaging a more reliable estimate of the shared response than competitors.

Lastly, we use independent component analysis as a basis to perform data augmentation. More precisely, we introduce CondICA, a data augmentation method that leverages a large amount of unlabeled fMRI data to build a generative model for labeled data using only a few labeled samples. CondICA yields an increase in decoding accuracy on eight large fMRI datasets.

Our main contributions consist in the reduction of SRM's training time as well as in the introduction of two more realistic models for the analysis of brain activity of subjects exposed to naturalistic stimuli : MultiViewICA and ShICA. Lastly, our results showing that ICA can be used for data augmentation are promising.

In conclusion, we present some directions that could guide future work. From a practical point of view, minor modifications of our methods could allow the analysis of resting state data assuming a shared spatial organization instead of a shared response. From a theoretical perspective, future work could focus on understanding how dimension reduction and shared response identification can be achieved jointly.

U N S U P E R V I S E D C O M P O N E N T A N A LY S I S F O R N E U R O I M A G I N G D ATA hugo richard

PhD thesis S Y N T H È S E

Cette thèse d'informatique et de mathématiques s'applique au domaine des neurosciences, et plus particulièrement aux recherches sur la modélisation de l'activité cérébrale humaine par électrophysiologie et imagerie. Dans ce champ, la tendance est actuellement d'expérimenter avec des stimuli naturels, comme le visionnage d'un film ou l'écoute d'une piste audio, et non plus avec des stimuli étroitement contrôlés mais outrageusement simples. L'analyse de ces stimuli « naturels » et de leurs effets demande toutefois de disposer d'une immense quantité d'images, par ailleurs très coûteuses. Sans outils mathématique, identifier l'activité neuronale à partir des données est quasi impossible. Toutefois, ces stimuli sont compliqués à modéliser et à analyser, car l'utilisation de méthodes fondées sur des régressions est limitée par la difficulté de modéliser les stimuli. C'est ce qui motive l'utilisation de méthodes non-supervisées qui ne font pas d'hypothèses sur ce qui déclenche les activations neuronales. Dans cette thèse, nous considérons d'abord le cas du modèle de réponse partagée (MRP), dans lequel les sujets sont supposés partager une réponse commune. Ce modèle est utile pour réduire la dimension des données, mais son entraînement est coûteux pour les données d'imagerie fonctionnelle (IRMf) dont la dimension peut être immense. Nous présentons une version bien plus rapide et beaucoup plus économe en mémoire. Mais le MRP fait des hypothèses irréalistes sur les données d'imagerie.

Des hypothèses plus réalistes sont utilisées dans l'analyse en composantes indépendantes (ACI) mais cette méthode est difficile à généraliser aux jeux de données qui contiennent plusieurs sujets. Nous proposons alors une extension de l'ACI appelée ACI multi-vue, fondée sur le principe de maximum de vraisemblance et qui convient à des jeux de données multi-sujets. L'ACI multi-vue a une vraisemblance en forme fermée qui peut être maximisée efficacement. Toutefois, cette méthode suppose la même quantité de bruit pour tous les sujets.

Nous présentons donc l'ACI partagée, une généralisation de l'ACI multi-vue qui s'accompagne d'un modèle de bruit plus général. Contrairement à presque tous les modèles fondés sur l'ACI, l'ACI partagée peut séparer des sources gaussiennes et non gaussiennes et propose une estimation optimale des sources communes (au sens des moindres carrés), qui pondère chaque sujet en fonction de son niveau de bruit estimé. En pratique, l'ACI partagée et l'ACI multi-vue permettent d'obtenir, en magnéto-encéphalographie et en IRMf, une estimation plus fiable de la réponse commune que leurs concurrents.

iii Enfin, nous utilisons l'ACI comme base pour faire de l'augmentation de données. Plus précisément, nous présentons l'ACI conditionnelle, une méthode d'augmentation de données qui exploite la grande quantité de données d'IRMf non étiquetées pour construire un modèle génératif en utilisant seulement un petit nombre de données étiquetées. L'ACI conditionnelle permet d'augmenter de façon appréciable la précision du décodage sur huit grands jeux de données d'IRMf.

Nos principaux apports nous semblent consister dans l'accélération de l'entraînement du MRP ainsi que dans l'introduction de deux modèles plus réalistes pour l'analyse de l'activité cérébrale de sujets exposés à des stimuli naturels : l'ACI multi-vue et l'ACI partagée. Enfin, nos résultats sont prometteurs concernant l'utilisation de l'ACI pour faire de l'augmentation de données.

Nous présentons pour finir quelques pistes qui pourraient guider des travaux ultérieurs. D'un point de vue pratique, des modifications de nos méthodes pourraient permettre l'analyse des données d'imagerie obtenues sur des sujets au repos en faisant l'hypothèse d'une organisation spatiale partagée à la place d'une réponse partagée. Une autre extension possible serait de prendre en compte le décalage temporel entre sujets que l'on observe en magnéto-encéphalographie. D'un point de vue théorique, les travaux futurs pourraient se concentrer sur la compréhension de la façon dont la réduction de dimensions et l'identification de la réponse partagée peuvent être réalisées conjointement.

iv A C K N O W L E D G M E N T S

Thank you to Pr. Moritz Große-Wentrup and Pr. Tülay Adali for reviewing my thesis and to Pr. Sylvain Chevalier, Pr Mathieu Kowalski, Pr Aapo Hyvärinen and Pr Alexandre Gramfort for being part of the jury. Their insightful remarks and comments have helped me improve the quality of this thesis.

Je remercie chaleuresement mon directeur de thèse Bertrand Thirion pour avoir largement contributé à mon épanouissement durant ces années de thèse. Merci Bertrand de m'avoir appris à écrire des papiers scientifiques, à présenter des résultats et à aborder des problèmes complexes avec méthode. Merci aussi pour ta très grande disponibilité et ton engagement dans mon travail.

Merci du fond du coeur Pierre, pour m'avoir transmis ta grande aisance à formaliser les problèmes, tes connaissances sur l'ICA, sur python et pour avoir investi du temps sur mes problèmes. I would also like to give special thanks to Aapo, who carefully reviewed our proofs and provided deep insights and always well-reasoned comments. It has been an honour to work with you. Merci Alex, la base de code et la pipeline d'intégration continue ne serait pas aussi propre sans ton aide. I deeply thank all other co-authors Luigi, Ana-Luisa, Jonathan, Lucas, Guillaume for your help in writing and great scientific discussions.

Merci à tous mes collègues et amis pour leur soutien. Merci à ceux qui ont partagé avec moi les joies du périple vers Neurospin: Thomas B., Jérome Alexis, Hamza, Binh, Antonia, Zac, Jerome, Arthur, Kamalaker. Merci aux volleyeurs Alexis T., Thomas M., Raphael, Louis, Maeliss, Demian qui m'ont permis de pratiquer au labo le sport que j'aime le plus au monde. Merci aussi à Olivier et Gaël pour ces discussions stimulantes. Merci à toute l'équipe. Merci à Corinne et Stephanie pour m'avoir guidé dans le labyrinthe administratif.

Merci à ma famille. Merci à ma mère, à mon père et à mon frère que j'aime fort et qui m'ont permi d'arriver jusque-là. Un grand merci à mon grand-père maternel Papou, qui malgré son grand âge a contribué significativement à l'écriture du résumé de cette thèse. Merci aussi à mon grand-père paternel Papet pour ses encouragements et son soutient. Enfin, merci à Laurie ma fiancée que j'admire et que j'aime qui m'a toujours accompagné, inspiré et soutenu dans les moments faciles comme dans les plus durs. We write vectors in bold letter v and scalars in lower case a. Upper case letters M are used to denote matrices. We denote |W| the absolute value of the determinant of W. x ∼ N(µ, Σ) means that x ∈ R k follows a multivariate normal distribution of mean µ ∈ R k and covariance Σ ∈ R k×k . When x ∼ N(µ, Σ), its density is given by x → N(x; µ, Σ). The j, j entry of a diagonal matrix Σ i is denoted Σ ij , the j entry of y i is denoted y ij . δ is the Kronecker delta. We use the usual scalar product for matrices A, B = tr(A B) and the associated norm is denoted A = A, A . Vectors can be seen as tall matrices and therefore the scalar product and the norm are the same as for matrices. The gradient of a real function f(x) ∈ R is denoted ∂f(x) ∂x and is seen as a column vector. The Jacobian of a vector valued function f(x) is denoted ∂f(x) ∂x and is a matrix such that the line j is given by ∂f(x j ) ∂x where x j is the j-th coordinate of x. When a subject is reading a sentence or when she is listening to a sentence, her brain activity is expected to differ. In order to measure where it differs and by which amount, one can perform a controlled experiment. The controlled experiment comes with a design matrix that describes the features driving her brain activation across time. In the above case, the occurrences of the subject listening to a sentence or reading a sentence are encoded in a design matrix. Then, a model is introduced to explain how the design matrix relates to brain activity. A simple model can consider that brain activity only differs when the task differs and therefore have a stereotypical representation of brain activity when she is reading a sentence and another one when she is listening to a sentence. While this model may simplify the reality (two repetitions of the same task would never yield exactly the same pattern or amount of brain activity), it is easy to interpret: the two stereotypical representations of each task can easily be compared and analyzed.

This procedure naturally extends to multiple subjects. Indeed, a model can be fit independently for each subject using the same design matrix.

Naturalistic stimuli

While controlled experiments give some insights about brain functionality, the subjects' experience is far from their every-day life. Naturalistic stimuli are meant to overcome this issue. Example of naturalistic stimuli include movie watching, music listening or resting (subjects are just asked to lie still in the scanner without further instruction). While there is a broad interest in understanding how the brain reacts in such ecological conditions, the recorded brain activity is difficult to analyze. In particular, design matrices are notoriously difficult to construct for naturalistic stimuli.

Component analysis

A possible solution is to learn the design matrix as part of the model. The widely used independent component analysis applied on the 1 data of one subject extracts a set of components that are maximally independent. Such components give a plausible design matrix as each component may be seen as a different set of features driving brain activity. However, many questions remain. How do we efficiently generalize such methods to multi-subject data ? How to measure their performance ? Why is it useful to construct a generative model for neuroimaging data ? In this thesis, we develop well principled unsupervised methods for component analysis of neuroimaging data.

organization of the manuscript

In part I, we give some background on statistical learning and optimization (chapter 2), neuroscience (chapter 3) and unsupervised methods popular in neuroimaging (chapter 4).

The parts II, III, IV and V highlight four different contributions that led to different publications.

Fast shared response model for fMRI data

When subjects are exposed to the same stimuli, their brain activity likely exhibits some common, or shared response. Our goal is to recover this shared response. The shared response model is one possible solution to this problem. However, the algorithm used in the shared response model does not scale well with the size of input . This is problematic because fMRI data have a very large size. Indeed, the fMRI data of each subject have a dimension on the order of p = 10 5 and a number of samples on the order of n = 10 3 . Therefore, there is a need for faster algorithms. This will be presented in part II.

MultiViewICA for neuroimaging data

In order to obtain a meaningful shared response, we need to impose constraints on the model. In the shared response model, the linear combination of the shared sources is done under orthogonality constraints which are often deemed not biologically plausible. A more biologically plausible constraint is to assume independence of the recovered responses. However, most popular methods using this assumption are partially heuristic when multiple subjects are involved. In part III, we introduce MultiViewICA: a method based on the maximum likelihood principle that can be efficiently optimized. In practice our method is able to recover more reliable estimate of the shared response on fMRI and MEG data than competitors. The model used in MultiViewICA allows the response of each subject to differ from the stereotypical response. However, it is assumed that the deviation is on the same order of magnitude for all subjects and all components. In Shared ICA, we use a more general model that allows to model different deviations from the stereotypical response depending on the subject considered and/or on the response. Importantly, the difference between subjects can be used as an additional source of information to recover an even better estimate of the shared response. In practice, we observe that Shared ICA improves upon MultiViewICA. This work will be presented in part IV. Our results suggest that ICA is a good generative model for fMRI data. Can it be used to perform data augmentation ? We design a data augmentation method, that leverages the large amount of (unlabeled) resting state data to generate realistic fake task data from a small set of images. Our data augmentation method yields an improvement in classification accuracy on eight large datasets. This will be discussed in part V.

chapter ordering

The appendices can be skipped at first read. Chapters 2, 3 and 4 present the necessary background. The reader more interested in the theory should first read chapters 5, 7, 9 and 11. Chapters 6, 8, 10 and 12 focus on practical results.

probabilistic generative models

In a probabilistic generative modeling framework, learning from experiments means identifying the underlying process that generates the data we observe.

More formally, consider X = x (1) , . . . , x (n) ∈ R k , n random vectors with joint density ν * X and consider the available data X ∈ R k,n as an observation of X. We also say that data X are generated from ν * X . The broad goal of probabilistic generative modeling is to recover ν * X from X. When samples x (1) , . . . , x (n) are independent and identically distributed, ν * X = ⊗ n i=1 ν * x so that X can be seen as n observations of the random variable x with density ν *

x . In the remaining of the thesis, we use the same notation x (i) for the random variable associated to the i-th sample and its observation (the i-th column of X). In addition, we use ν * to denote both ν *

x and ν * X depending on whether samples are independent and identically distributed or not.

In practice, we assume a model for the true density ν * meaning that we assume that ν * belongs to a family of densities F. When it is indeed true that ν * ∈ F, we say that the model holds.

Often, we assume F to be a set of parametric densities so that any density in F can be written as ν θ where θ ∈ Θ is a set of parameters and Θ is the set of all possible θ. When the model holds, there exists an optimal set of parameters θ * such that ν * = ν θ * . The goal is then to find θ * .

Before we even try to find θ * , it is instructive to wonder whether the problem is well defined. Ideally we would like our model to be such that if ν θ 1 = ν θ 2 then θ 1 = θ 2 . When this is the case, we say that the model is identifiable.

In all the proofs in this section, we assume that the integration and differentiation operators can always be exchanged and that all quantities introduced are well defined. The set of parameters θ is viewed as a vector.

Note that our definition of a model includes the latent variable model (which specifies the dependencies across the different random variables) and the densities of each random variable. However it does not include the estimation procedure. Including the densities in the model allows us to state claims like : "this model is identifiable" which would not always be possible if the densities were not included. Yet, a model is most useful when it comes with an efficient estimation procedure so we often present together the model and the algorithm used to estimate its parameters.

Desirable properties of estimators

An estimator θn of θ * is a function of the observations X = x (1) , . . . , x (n) that aims at finding θ * . An estimator is a random variable, as it depends on X. Therefore, it can almost never be perfectly accurate and that is why we need a criterion to measure its inaccuracy.

A common choice is the mean squared error criterion given by:

E[ θn -θ * 2 ] = E[ θn -E[ θn ] + E[ θn ] -θ * 2 ] (2.1) = E[ θn -E[ θn ] 2 ] + E[ θn ] -θ * 2 (2.2) = tr(V( θn )) + E[ θn ] -θ * 2 (2.3)
where the right hand side in equation 2.3 gives the bias-variance decomposition. The left term is the trace of the covariance V( θn

) = Cov( θn -E[ θn ], θn -E[ θn ]) where Cov(a, b) = E[(a -E[a])(b - E[b])
] and the right term is the squared norm of the bias given by

E[ θn ] -θ * .
The norm of the bias can be minimized exactly and such estimators that achieve E[ θn ] -θ * = 0 are called unbiased.

In Example 1, we study the bias of the sample mean and sample variance. The sample mean is shown to be unbiased. In contrast, the sample variance is biased: we show how to correct the estimator of the variance so that it becomes unbiased.

Example 1 (Biased and unbiased estimate of the parameters of a 1D Gaussian). Consider n observations x (1) , . . . x (n) of x with mean µ * and variance σ 2 * . Consider the sample mean: μe = 1 n n i=1 x (i) . This estimate is unbiased as

E[ μe ] = 1 n n i=1 E[x (i) ] = µ * . Consider the sample variance: σ2 e = 1 n n i=1 (x (i) -n z=1 1 n x (z) ) 2 .
We have

E[ σ2 e ] = 1 n n i=1 E[(x (i) - n z=1 1 n x (z) ) 2 ] (2.4) = 1 n n i=1 E[(x (i) ) 2 -2x (i) ( n z=1 1 n x (z) ) + ( n z=1 1 n x (z) ) 2 ] (2.5) = 1 n n i=1 E[(x (i) ) 2 -2x (i) ( n z=1 1 n x (z) ) + 1 n 2 n y=1,z=1
x (y) x (z) ]

(2.6)

= 1 n n i=1 (σ 2 * + µ 2 * - 2 n σ 2 1 n σ 2 * (2.8) = n -1 n σ 2 * (2.9)
so the sample covariance is biased. In contrast, if we consider the estimator

σ2 u = n n -1 σe 2 (2.10)
we can see that it is unbiased.

However the variance of an estimator cannot be arbitrarily low as shown by Proposition 1 (Cramer-Rao bound).

Proposition 1 (Cramer-Rao bound). Let θn be an estimator of θ * . Then,

V( θn ) ∂E[ θn ] ∂θ I n (θ * ) -1 ( ∂E[ θn ] ∂θ )
where A B is understood as A -B is positive semi-definite. We introduced I n (θ), the Fisher information matrix given by

I n (θ) = E[ ∂ log(ν θ (X)) ∂θ ∂ log(ν θ (X)) ∂θ ] (2.11)
Lastly, the quantity l(x, θ) = log(ν θ (x)) is called log-likelihood of x and its derivative ψ(x, θ)

= ∂ log(ν θ (x)) ∂θ is called the score function of x.
Proof of Cramer-Rao bound. First let us show that when the optimal parameter is used, the expected score function cancels:

E X [ψ(θ * )] = E X [ ∂ log(ν θ * (X)) ∂θ ] (2.12) = E X [ 1 ν θ * (X) ∂ν θ * (X) ∂θ ] (2.13) = X 1 ν θ * (X) ∂ν θ * (X) ∂θ ν θ * (X)dX (2.14) = X ∂ν θ * (X) ∂θ dX (2.15) = ∂ X ν θ * (X)dX ∂θ (2.16) = ∂1 ∂θ (2.17) = 0 (2.18) so that V[ψ(θ * )] = I n (θ * ).
We also have that:

Cov( θn , ψ(θ * )) = E[( θn -E[ θn ])(ψ(θ * ) -E[ψ(θ * )]) ] (2.19) = E[ θn ψ(θ * ) ] (2.20) = X θn ∂ν θ * (X) ∂θ dX (2.21) = ∂ X θn ν θ * (X)dX ∂θ (2.22) = ∂E[ θn ] ∂θ (2.23)
and similarly Cov(ψ(θ * ), θn ) = ( ∂E[ θn ] ∂θ ) . Then we apply the following Cauchy Schwartz inequality for random vectors:

∀x, y Var(y) Cov(y, x) Var(x) -1 Cov(x, y) (2.24)
(see a proof in [START_REF] Tripathi | A matrix extension of the Cauchy-Schwarz inequality[END_REF]) and therefore get the expected result:

Var( θn ) ∂E[ θn ] ∂θ I n (θ * ) -1 ( ∂E[ θn ] ∂θ ) (2.25)
When an estimator is both unbiased and reaches the Cramer-Rao bound we call the estimator efficient. If we assume that samples are independent and identically distributed we have: ψ(X, θ) = n i=1 ψ(x (i) , θ) and therefore

I n (θ * ) = E[ψ(X, θ * )ψ(X, θ * ) ] (2.26) = n i=1 E[ψ(x (i) , θ * )ψ(x (i) , θ * ) ] (2.27) = nI(θ * ) (2.28)
where

I(θ * ) = E x∼ν θ * [ψ(x, θ * )ψ(x, θ * ) ] (2.29)
When samples are independent and identically distributed and θn is unbiased the Cramer-Rao bound is given by: Var( θn )

1 n I(θ * ) -1 (2.30)
The Example 2 hows that the sample mean is efficient.

Example 2 (Sample mean is efficient). Consider n observations x (1) , . . . x (n) of x generated from ν µ * (x) = N(x; µ * , σ 2 * ) and consider the sample mean μe = 1 n n i=1 x (i) . The variance is given by:

V[ μe ] = V[ 1 n n i=1 x (i) ] (2.31) = σ 2 * n (2.32)
The Fisher information matrix is given by

I = E[( ∂ log(ν µ * (x)) ∂µ ) 2 ] (2.33) = E[( ∂ -1 2σ 2 * (x -µ * ) 2 -1 2 log(2πσ 2 * ) ∂µ ) 2 ] (2.34) = E[(- 1 σ 2 * (µ * -x)) 2 ] (2.35) = 1 σ 4 * E[(µ * -x) 2 ] (2.36) = 1 σ 2 *
(2.37)

and therefore we have

V[ μe ] = 1 n I -1 (2.38)
so the sample mean is efficient.

In practice, efficient estimators are extremely rare. In the next section, we introduce the maximum likelihood estimator which is not always unbiased nor efficient in the finite sample case, but satisfies these properties asymptotically.

Maximum likelihood

The maximum likelihood estimates the parameters θ such that the density ν θ at x (1) , . . . , x (n) is the highest among all possible values for θ ∈ Θ:

θn = argmax θ∈Θ ν θ (x (1) , . . . , x (n) ) (2.39)
The quantity ν θ (x (1) , . . . , x (n) ) is called the likelihood and we call θn the maximum likelihood estimator.

It is often assumed that samples are independent and identically distributed so that the joint density can be written: ν θ (x (1) , . . . ,

x (n) ) = n i=1 ν θ (x (i) )
Let us define the empirical expected log-likelihood:

l n (θ) = 1 n n i=1 l(x (i) , θ) (2.40) = E n l(x, θ) (2.41)
where E n is the empirical expectation operator defined by

E n f(x) = 1 n n i=1 f(x (i) ) (2.42)
The next lines show that finding the maximum likelihood is done by optimizing the empirical expected log-likelihood.

θn = argmax θ∈Θ ν θ (x (1) , . . . , x (n) ) (2.43) = argmax θ∈Θ n i=1 ν θ (x (i) ) (2.44) = argmax θ∈Θ n i=1 log(ν θ (x (i) )) (2.45) = argmax θ∈Θ 1 n n i=1 l(x (i) , θ) (2.46) = argmax θ∈Θ l n (θ) (2.47) (2.48)
By the law of large numbers, as the number of samples increases, the empirical expected likelihood converges almost surely to the expected log-likelihood:

l(θ) = E x [l(x, θ)].
Example 3 shows that, under a Gaussian model, the maximum likelihood estimator of the mean and the variance is the sample mean and sample variance respectively. Example 3. Consider n observations x (1) , . . . , x (n) of an unknown random variable x and consider the model x ∼ N(µ, σ 2 ). The empirical expected log-likelihood is:

l n (µ, σ 2 ) = - 1 n n i=1 1 2σ 2 (x (i) -µ) 2 - 1 2 log(2πσ 2 ) (2.49)
First order conditions yield:

n i=1 (µ -x (i) ) = 0 ⇐⇒ µ = 1 n n i=1 x (i) (2.50) 1 n n i=1 1 (σ 2 ) 2 (x (i) -µ) 2 - 1 σ 2 = 0 ⇐⇒ σ 2 = 1 n n i=1 (x (i) -µ) 2 (2.51)
So the maximum likelihood estimators of the mean and the variance are the sample mean and sample variance respectively.

The maximum likelihood estimator θn is a random variable as it depends on X = x (1) , . . . , x (n) . As the number of samples n increases, we expect θn to get closer to the optimal set of parameters θ * . In Proposition 2, it is shown that the maximum likelihood estimator is consistent meaning θn

P ----→ n→∞ θ * where θn P ----→ n→∞ θ * denotes the convergence in probability defined by ∀ > 0, p( θn -θ * < ) ----→ n→∞ 1.
Proposition 2 (Consistency of the maximum likelihood estimator). Assume Θ is compact, Assume l n converge uniformly in probability to l and assume l is continuous. Lastly, assume that the model is identifiable. Then, the maximum likelihood estimator θn is consistent meaning that θn

P ----→ n→∞ θ * .
where uniform convergence in probability means sup θ l n (θ) -

l(θ) P ----→ n→∞ 0.
Proof. We first show that l is maximum at θ * :

l(θ * ) -l(θ) = E x∼ν θ * [log ν θ * (x) -log(ν θ (x)]
(2.52)

= D KL (ν θ * , ν θ ) (2.53) 0 (2.54)
where D KL is the Kullback-Leibler divergence that is always positive. The maximum is unique. This comes from the identifiability of the model that implies

l(θ * ) = l(θ) =⇒ θ * = θ (2.55) Let ε > 0 and define V ε = {θ, θ -θ * < ε} an open neighborhood of θ * . Because Θ is compact and V ε open, Θ ∩ V C ε is compact and since l is continuous, max θ∈Θ∩V C ε l(θ) is reached for a value θ 0 ∈ Θ ∩ V C ε .
Let us define δ = l(θ * ) -l(θ 0 ) > 0 and consider the events

A n = sup θ∈Θ∩V C ε l n (θ) -l(θ) < δ 2 (2.56) B n = sup θ∈V ε l n (θ) -l(θ) < δ 2 (2.

57)

We have

A n =⇒ ∀θ ∈ Θ ∩ V C ε , l n (θ) -l(θ) < δ 2 (2.58) =⇒ ∀θ ∈ Θ ∩ V C ε , l n (θ) < δ 2 + l(θ 0 ) (2.59) =⇒ ∀θ ∈ Θ ∩ V C ε , l n (θ) < - δ 2 + l(θ * ) (2.60)
and

B n =⇒ ∀θ ∈ V ε , l n (θ) > - δ 2 + l(θ) (2.61) =⇒ l n (θ * ) > - δ 2 + l(θ * ) (2.62) (2.63)
Then consider θn = argmax θ∈Θ (l n (θ)) and assume θn ∈ V C ε ∩ Θ. We have

A n ∩ B n =⇒ l n ( θn ) < - δ 2 + l(θ * ) < l n (θ * ) (2.64) So if θn ∈ V C ε ∩ Θ, A n ∩ B n contradicts the fact that θn = argmax θ∈Θ (l n (θ)). Therefore A n ∩ B n =⇒ θn ∈ V ε . From the uniform convergence in probability of l n to l we have p(A n ∩ B n ) ----→ n→∞ 1 and therefore p( θn ∈ V ε ) ----→ n→∞ 1 which means θn P ----→ n→∞ θ * .
As we have ween in Example 1, the sample variance is biased. However, in the large sample limit the bias disappears. Is this a consequence of consistency ? Intuitively, consistency seems to be a stronger condition than asymptotic unbiasedness since consistency implies convergence of the random variable θn whereas asymptotic unbiasedness only implies convergence of the mean E[ θn ]. Proposition 3 shows that this intuition is correct as long as the variance is bounded.

Proposition 3 (Consistency implies asymptotic unbiasedness). Assume

θ n P ----→ n→∞ θ * and assume ∃M ∈ R, E[ θ n -θ * 2 ] < M. Then E[ θ n -θ * ] ----→ n→∞ 0
Proof. Set ε > 0, we have

θ n P ----→ n→∞ θ * =⇒ p( θ n -θ * > ε 2 ) ----→ n→∞ 0 (2.65) =⇒ ∃N ∈ N, ∀n > N, p( θ n -θ * > ε 2 ) < ε M (2.66)
So ∀n > N we get: Proposition 4 (Asymptotic normality of maximum likelihood estimators). We assume the same hypothesis as in Proposition 2. We further assume that θ * is in the interior of Θ. Lastly, denoting

E[ θ n -θ * ] = E[ θ n -θ * 1 θ n -θ * ε 2 ] + E[ θ n -θ * 1 θ n -θ * > ε 2 ] (2.67) ε 2 + E[ θ n -θ * 2 ]p( θ n -θ * > ε 2 ) (2.
J n (θ n ) = ∂ 2 l n ∂θ 2 (θ n ) and J(θ * ) = ∂ 2 l ∂θ 2 (θ * ) we assume that θ n P ----→ n→∞ θ * =⇒ J n (θ n ) P ----→ n→∞ J(θ * ) Then, √ n θn -θ * D ----→ n→∞ N(0, I(θ * ) -1 )
where the convergence in distribution

x (n) D ----→ n→∞ x means ∀t, F x (n) (t) ----→ n→∞ F x (t) where F x is the distribution function of x.
Proof. First order conditions give ∂l n ( θn ) ∂θ = 0. From the mean value theorem there exists θ 0 such that

∂l n ( θn ) ∂θ - ∂l n (θ * ) ∂θ = J n (θ 0 )( θn -θ * ) (2.70) ⇐⇒ - 1 √ n n i=1 ψ θ * (x (i) ) = J n (θ 0 ) √ n( θn -θ * ) (2.71)
We know from the proof of the Cramer-Rao bound (Proposition 1) that E x [ψ θ * (x)] = 0. We can use the central limit theorem and write:

1 √ n n i=1 ψ θ * (x (i) ) D -→ N(0, I(θ * )) (2.72)
Then, we have: From Slutsky's theorem we have

J(θ * ) = E x [ ∂ψ θ * (x) ∂θ ] (2.73) = E x [ ∂ 1 ν θ * (x) ∂ν θ * (x) ∂θ ∂θ ] (2.74) = E x [ 1 ν θ * (x) ∂ 2 ν θ * (x) ∂θ 2 - ∂ν θ * (x) ∂θ ν θ * (x) ∂ν θ * (x) ∂θ ν θ * (x) ] (2.75) = -I(θ * ) (2.76) (2.
I(θ * ) √ n( θn -θ * ) D -→ N(0, I(θ * )) (2.78) ⇐⇒ √ n( θn -θ * ) D -→ N(0, I(θ * ) -1 ) (2.79)
Looking at the variance of θnθ * , we see that it behaves like 1 n I(θ * ) -1 which is the quantity in the Cramer-Rao bound for unbiased estimators (see Proposition 1). Because of this, the maximum likelihood estimator is called asymptotically efficient.

Some shortcomings

It may look like an efficient and unbiased estimator should always yield the best possible mean squared error. However, this is not the case. A striking example is given by the Stein paradox [START_REF] Stein | Inadmissibility of the Usual Estimator for the Mean of a Multivariate Normal Distribution[END_REF]. As can be seen in Example 4, there exists an estimator of the mean that achieves a strictly better mean squared error than the sample mean. However, we have shown that the sample mean (which is also the maximum likelihood estimator) is unbiased and efficient. Stein's estimate shows that biased estimators can sometimes achieve lower mean squared error than unbiased ones.

Example 4 (Stein's estimate of the mean). Let us consider a single observation x generated from a multivariate Gaussian of dimension v 3: N(µ, I) where I is the identity matrix and µ ∈ R v is the unknown mean that we want to estimate. As in the one-dimensional case, the sample mean e(x) = x is unbiased and efficient. The mean squared error is given by:

E[(e(x) -µ) 2 ] = tr(V(e(x))) = v (2.80)
Now consider the estimator s(

x) = x -(v -2) x x 2 .
The mean squared error is given by:

E[(s(x) -µ) 2 ] (2.81) = E[(x -(v -2) x x 2 -µ) 2 ] (2.82) = E[(x -µ) 2 ] -2(v -2)E[ (x -µ) x x 2 ] + (v -2) 2 E[ 1 x 2 ] (2.83) = E[(x -µ) 2 ] -2(v -2) v j=1 E[ (x j -µ j )x j x 2 ] + (v -2) 2 E[ 1 x 2 ]
(2.84) where x j is the j-th coordinate of x By integration by part

E[ (x j -µ j )x j x 2 ] = E[ ∂ x j x 2 ∂x j ] (2.85) = E[ 1 x 2 - 2x 2 j x 4 ] (2.86) so that E[(s(x) -µ) 2 ] (2.87) = E[(x -µ) 2 ] -2(v -2) v j=1 E[ 1 x 2 - 2x 2 j x 4 ] + (v -2) 2 E[ 1 x 2 ]
(2.88)

= E[(x -µ) 2 ] -2(v -2)E[ v x 2 - 2 v j=1 x 2 j x 4 ] + (v -2) 2 E[ 1 x 2 ] (2.89) = E[(x -µ) 2 ] -2(v -2) 2 E[ 1 x 2 ] + (v -2) 2 E[ 1 x 2 ] (2.90) = v -(v -2) 2 E[ 1 x 2 ] (2.91)
Therefore s(x) always yields lower expected mean squared error than e(x).

Example 4 shows that shrinking the sample mean towards the origin decreases the mean squared error of the estimate in high dimensions. A similar technique can also be used for estimating the covariance of variables. For instance, Ledoit and Wolf showed in [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF] that shrinking the sample covariance towards identity decreases lower expected mean squared error in high dimensions.

optimization

The maximum likelihood estimator gives a natural way to obtain estimators that are consistent and asymptotically efficient. However, it requires finding the maximum of the empirical expected log-likelihood. This can rarely be done using a closed form formula and one almost always have to resort to iterative methods.

Some iterative optimization algorithms

Let us consider a function f : R v → R that we want to minimize. f(θ), where θ ∈ R v , can be seen as the negative expected likelihood.

Optimization algorithms that only use first order derivatives to make a step are called first order methods. We will begin by presenting the famous gradient descent. Then, we move on to second order methods with Newton and quasi-Newton methods. This section closely follows the chapter 9 of [START_REF] Boyd | Convex optimization[END_REF].

Gradient descent

We assume that f is differentiable everywhere. From a point θ 0 , assuming a small step size α and a direction d such that d = 1, a Taylor decomposition at first order gives:

f(θ 0 + αd) = f(θ 0 ) + ∂f(θ 0 ) ∂θ , αd + o(α) (2.92)
The best direction is the one that minimizes f(θ 0 + αd). If we neglect the terms in o(α), it is given by argmin

d, d =1 f(θ 0 ) + ∂f(θ 0 ) ∂θ , αd = -∂f(θ 0 ) ∂θ ∂f(θ 0 ) ∂θ (2.93)
Therefore, gradient descent updates are given by:

θ k+1 = θ k -α ∂f(θ k ) ∂θ (2.94)
where α is a small quantity.

Under certain conditions that we specify in Proposition 5, gradient descent converges to the minimum.

Proposition 5 (Convergence of gradient descent). Assume that f is twice differentiable and µ-strongly convex:

∂ 2 f(θ) ∂θ 2 µI (2.95)
In addition, assume that f is -smooth:

∂ 2 f(θ) ∂θ 2 I (2.96)
where I is the identity matrix. From θ 0 ∈ R v and given α ∈ R such that

1 α > 0, the iterates θ k+1 = θ k -α ∂f(θ k )
∂θ converge to the minimum θ * according to

θ k+1 -θ * 2 (1 -αµ) k+1 θ 0 -θ * 2 (2.97)
Proof. We have

θ k+1 -θ * 2 = θ k -α ∂f(θ k ) ∂θ -θ * 2 (2.98) = θ k -θ * 2 -2 θ k -θ * , α ∂f(θ k ) ∂θ + α ∂f(θ k ) ∂θ 2 (2.99)
From µ-strong convexity and Lagrange inequality we get:

f(θ * ) f(θ k ) + ∂f(θ k ) ∂θ k , θ * -θ k + µ 2 θ * -θ k 2 (2.100) so that -2α ∂f(θ k ) ∂θ k , θ k -θ * -2α(f(θ k ) -f(θ * )) -µα θ * -θ k 2 (2.101)
and using the fact that f is -smooth we have:

∀θ, y, f(y) f(θ) + ∂f(θ) ∂θ (y -θ) + 2 y -θ 2 (2.102) =⇒ ∀θ, f(θ - 1 ∂f(θ) ∂θ ) -f(θ) - 1 2 ∂f(θ) ∂θ 2 (2.103) =⇒ ∀θ, f(θ * ) -f(θ) - 1 2 ∂f(θ) ∂θ 2 (2.104) =⇒ α ∂f(θ k ) ∂θ 2 2α 2 l(f(θ k ) -f(θ * )) (2.105)
So from (2.99) using the inequalities (2.101) and (2.105) we get

θ k+1 -θ * 2 (1 -αµ) θ k -θ * 2 -2α(1 -α )(f(θ k ) -f(θ * ))
(2.106) and since 0 < α < 1 we get that

θ k+1 -θ * 2 (1 -αµ) θ k -θ * 2 (2.107)
which by induction yields the desired result.

The type of convergence we get in Proposition 5 is called a linear convergence (because θ k+1 -θ * is bounded by a linear function of θ kθ * ) .

Newton method and quasi-Newton methods

We assume that f is twice differentiable everywhere. As in gradient descent, we depart from a point θ 0 , assume a small step size α and consider a direction d such that d = 1. A Taylor decomposition at the second order gives:

f(θ 0 + αd) = f(θ 0 ) + ∂f(θ 0 ) ∂θ , αd + αd, ∂ 2 f(θ 0 ) ∂θ 2 αd + o(α 2 ) (2.108)
If we neglect the terms in o(α 2 ), the best direction is given by argmin

d, d =1 f(θ 0 ) + ∂f(θ 0 ) ∂θ , αd + αd, ∂ 2 f(θ 0 ) ∂θ 2 αd (2.109) = - ( ∂ 2 f(θ 0 ) ∂θ 2 ) -1 ∂f(θ 0 ) ∂θ ( ∂ 2 f(θ 0 ) ∂θ 2 ) -1 ∂f(θ 0 ) ∂θ (2.110)
Therefore, Newton updates are given by:

θ k+1 = θ k -α( ∂ 2 f(θ k ) ∂θ 2 ) -1 ∂f(θ k ) ∂θ (2.111)
where α is a small quantity. Unfortunately, Newton's method is not guaranteed to converge. This contrasts with the gradient descent method that is guaranteed to converge when the step-size is small enough. In order to converge, Newton's method needs that each step yields a sufficient decrease of the loss. This is done by a line search.

The exact line search

sets α = argmin t∈[0,1] f(θ k -t( ∂ 2 f(θ 0 )
∂θ 2 ) -1 ∂f(θ k ) ∂θ ) while the backtracking line search is an iterative procedure where one starts with α = 1 and repeatedly halves α until f(θ k+1 ) < f(θ k ).

As shown by Proposition 6, Newton's method is guaranteed to converge when an exact line search is used.

Proposition 6 (Convergence of Newton's method). Assume that f is twice differentiable, µ-strongly convex and -smooth:

µI ∂ 2 f(θ) ∂θ 2 I (2.112)
where I is the identity matrix. Assume that the Hessian is Lipschitz with constant h:

∂ 2 f(θ) ∂θ 2 - ∂ 2 f(θ ) ∂θ 2 h θ -θ (2.113) From θ 0 ∈ R v the iterates θ k+1 = θ k -α( ∂ 2 f(θ 0 ) ∂θ 2 ) -1 ∂f(θ k ) ∂θ
where α is chosen using an exact line search converge to the minimum θ * . Depending on the norm of the gradient two phases exist:

• The damped phase where ∂f(θ) ∂θ µ 2 h . In this phase we have f(θ k+1 ) -f(θ k ) µ 5 2h 2 2 . Therefore the number of iterations in this phase i cannot be larger than 2h 2 2 µ 5 (f(θ 0 ) -f(θ * )).

• The pure Newton phase where ∂f(θ) ∂θ

< µ 2 h . In this phase f(θ k ) - f(θ * ) 2µ 3 h 2 ( 1 2 ) 2 k-i+1
Once the Newton phase is reached, it continues until convergence of the algorithm.

Proof. In the damped phase we have ∂f(θ) ∂θ µ 2 h . From Lagrange inequality and -smoothness, denoting

d k = -( ∂ 2 f(θ) ∂θ 2 ) -1 ∂f(θ)
∂θ we have:

f(θ k + td k ) f(θ) + t ∂f(θ) ∂θ , d k + t 2 2 d k 2 (2.114) f(θ) + t ∂f(θ) ∂θ , d k - t 2 2µ ∂f(θ) ∂θ , d k (2.115) (2.116)
We then have

f(θ k + αd k ) f(θ) + t ∂f(θ) ∂θ , d k - t 2 2µ ∂f(θ) ∂θ , d k (2.117) (2.118) setting t = µ gives f(θ k + αd k ) f(θ) + µ 2 ∂f(θ) ∂θ , d k (2.119) f(θ) - µ 2 2 ∂f(θ) ∂θ 2 (2.120) f(θ) - µ 5 2h 2 (2.121) (2.122)
which gives the desired result.

In the pure Newton phase, we assume α = 1 (the exact line search can only be better than this). We have

∂f(θ k+1 ) ∂θ = ∂f(θ k + d k ) ∂θ - ∂f(θ k ) ∂θ - ∂ 2 f(θ k ) ∂θ 2 d k (2.123) = t∈[0,1] ∂ 2 f(θ k + td k ) ∂θ 2 d k dt - ∂ 2 f(θ k ) ∂θ 2 d k (2.124) = t∈[0,1] ( ∂ 2 f(θ k + td k ) ∂θ 2 - ∂ 2 f(θ k ) ∂θ 2 )d k dt (2.125) (2.126) Then ∂f(θ k+1 ) ∂θ t∈[0,1] ( ∂ 2 f(θ k + td k ) ∂θ 2 - ∂ 2 f(θ k ) ∂θ 2 )d k dt (2.127) t∈[0,1] ( ∂ 2 f(θ k + td k ) ∂θ 2 - ∂ 2 f(θ k ) ∂θ 2 ) d k dt (2.128) t∈[0,1] th d k 2 dt (2.129) h 2µ 2 ∂f(θ k ) ∂θ 2 (2.130) (2.131) Since ∂f(θ k ) ∂θ µ 2
h the sequence of gradients converges to zero. The convergence is quadratic. At iteration k we therefore have:

h 2µ 2 ∂f(θ k ) ∂θ ( 1 2 ) 2 k-i (2.132) (2.133)
Using strong convexity we conclude that

f(θ k ) -f(θ * ) 1 2µ ∂f(θ k ) ∂θ 2 2 µ 3 h 2 ( 1 2 ) 2 k-i (2.134)
In the pure Newton phase, the convergence is quadratic. This is a very strong advantage of Newton methods. A second advantage of Newton method is its equivariance as demonstrated in Proposition 7. Equivariance means that for any invertible matrix A, working with parameters Aθ instead of parameters θ has no impact on the result given by the algorithm.

Proposition 7 (Newton algorithms are equivariant). Let us consider θ k the current estimate of the minimum of f obtained by the Newton method after k iterations starting from θ 0 . For any invertible matrix A, consider y k the current estimate of the minimum of g(y) = f(Ay) obtained by the Newton method after k iterations starting from y 0 = A -1 θ 0 . Then,

y k = A -1 θ k .
Proof. The relation holds for k = 0 and assuming y k = A -1 θ k we have:

y k+1 = y k -α( ∂ 2 g(y k ) ∂y 2 ) -1 ∂g(y k ) ∂y (2.135) = y k -αA -2 ( ∂ 2 f(θ k ) ∂θ 2 ) -1 A ∂f(θ k ) ∂θ (2.136) = A -1 (θ k -α( ∂ 2 f(θ k ) ∂θ 2 ) -1 ∂f(θ k ) ∂θ ) (2.137) = A -1 θ k+1 (2.138) (2.139)
In contrast, gradient descent is not equivariant. However, despite its attractive properties, Newton's method is rarely used in practice because it is often intractable. Indeed inverting the Hessian ∂ 2 f ∂θ 2 is difficult as it is an operator in dimension R v×v×v×v . In some cases, it is possible to construct an approximation of the inverse of the Hessian in a reasonable time. This can be done iteratively by building a sequence of matrices B k that approaches the inverse of the Hessian as k grows. Methods that use an approximation of the inverse instead of the true inverse are called quasi-Newton methods [START_REF] Philip | Quasi-Newton methods for unconstrained optimization[END_REF].

EM and generalized EM

In the maximum likelihood framework, our goal is to find the parameters θ * that maximize the expected likelihood l(θ) of observations X = x (1) , . . . , x (n) assuming a model ν θ (X). In the previous section, we have presented some iterative optimization methods that can be used to directly maximize the log-likelihood. In this section, we present an alternative technique suited to latent variable models. Latent variable models include, in addition to observed data x, unobserved data z called latent variables. The observed data x is assumed to depend on the unobserved data z. With some abuse of notation we denote ν θ (z) the likelihood of z, ν θ (x|z) the density of x|z evaluated at x, z and ν θ (x, z) = ν θ (z)ν θ (x|z) the joint likelihood of x, z called the complete likelihood.

We can relate the likelihood and the completed likelihood by:

ν θ (x) = z ν θ (x, z)dz (2.140)
The EM algorithm maximizes an expression that depends on the complete log-likelihood rather than the log-likelihood so it is useful when the former is much simpler to maximize than the later. This section follows the work in [START_REF] Radford | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF].

In order to make this more concrete, let us focus on the Gaussian mixture model in Example 5.

Example 5 (Gaussian mixture models). Let us consider a latent variable z sampled from a Bernoulli distribution with parameter φ, and x is given by x|z = 0 ∼ N(µ 0 , 1) and and x|z = 1 ∼ N(µ 1 , 1).

We call x (1) , . . . , x (n) the observed samples and z (1) , . . . , z (n) the corresponding unobserved latent variables.

The log-likelihood of x (i) is given by:

l(x (i) , (µ 0 , µ 1 , φ)) (2.141) = log(φN(x (i) , µ 1 , 1) + (1 -φ)N(x (i) , µ 0 , 1)) (2.142) = log φ exp(-(x (i) -µ 1 ) 2 2 ) √ 2π + (1 -φ) exp(-(x (i) -µ 0 ) 2 2 ) √ 2π (2.143)
Whereas the completed log-likelihood of (x (i) , z (i) ) is given by

l((x (i) , z (i) ), (µ 0 , µ 1 , φ)) (2.144) = log((φN(x (i) , µ 1 , 1)) z (i) ((1 -φ)N(x (i) , µ 0 , 1)) 1-z (i) ) (2.145) = z (i) (- (x (i) -µ 1 ) 2 2 + log(φ)) (2.146) + (1 -z (i) )(- (x (i) -µ 0 ) 2 2 + log(1 -φ)) + c (2.147)
where c is a constant that does not depend on µ 1 or µ 2 .

It is easy to see that the expected completed log-likelihood will be much easier to maximize than the expected log-likelihood.

This example shows that unsurprisingly, it would be easier to find the parameters of a Gaussian mixture model if we knew the component from which are generated each sample.

Instead of optimizing the expected log-likelihood l(θ) directly, the EM algorithm optimizes the following function:

F(q, θ) = E q [log(ν θ (x, z))] + H q (2.148)
where q is a density and

E q [log(ν θ (x, z))] = z log(ν θ (x, z))q(z)dz (2.149) H q = z -log(q(z))q(z)dz (2.150)
H q is called the entropy of q and is always positive. Let us introduce ν θ (z|x) the density of z|x, we have:

ν θ (x, z) = ν θ (z|x)ν θ (x)
and therefore

F(q, θ) = E q [log(ν θ (z|x)) + log(ν θ (x))] + E q [-log(q)] (2.151) = log(ν θ (x)) -E q [log(q) -log(ν θ (z|x))] (2.152) = l(x, θ) -D KL (q, ν θ (z|x)) (2.153)
From the fact that D KL is positive and

D KL (a, b) = 0 ⇐⇒ a = b we have that F(q, θ) l(x, θ) (2.154)
F(ν θ (z|x), θ) = l(x, θ) (2.155)
and therefore max q,θ F(q, θ) = max θ (max

q F(q, θ)) = max θ l(x, θ) (2.156)
Any EM algorithm maximizes l by maximizing F. The most common practice is to maximize alternatively F with respect to q and θ. At iteration k, let us call θ k the current estimate of θ. According to equation (2.155), the maximum with respect to q is given by q = ν θ k (z|x). Then, we have to maximize equation (2.148) with respect to θ. As the entropy H does not depend on θ the function to maximize is given by:

Q(θ) = E z∼ν θ k (z|x) [log(ν θ (x, z))] (2.157)
Computing Q is called the E-step. Then we maximize Q and set

θ k+1 = argmax θ Q(θ) (2.158)
This step is called the M-step.

In Example 6 we us use the EM algorithm to optimize the Gaussian mixture model introduced in Example 5.

Example 6 (Optimizing the Gaussian mixture via EM). Let us take the same data as in Example 5. We call µ 0 k , µ 1 k , φ k the estimates of µ 0 , µ 1 , φ at iteration k and define:

γ (i) k = p(z (i) = 1|x (i) ) = φ k N(x (i) ; µ 0 k , 1) φ k N(x (i) ; µ 0 k , 1) + (1 -φ k )N(x (i) ; µ 1 k , 1) (2.159)
The E-step is given by:

Q(µ 0 , µ 1 , φ) = n i=1 γ (i) k log(ν µ 0 ,µ 1 ,φ (x (i) , 1))
(2.160)

+ (1 -γ (i) k ) log(ν µ 0 ,µ 1 ,φ (x (i) , 0)) (2.161) = n i=1 γ (i) k (- (x (i) -µ 1 ) 2 2 + log(φ)) (2.162) + (1 -γ (i) k )(- (x (i) -µ 0 ) 2 2 + log(1 -φ)) + c (2.163)
The M-step is given by:

∂Q ∂φ (φ k+1 ) = 0 ⇐⇒ φ k+1 = n i=1 γ (i) k n (2.164) ∂Q ∂φ (µ k+1 0 ) = 0 ⇐⇒ µ k+1 0 = n i=1 (1 -γ (i) k )x (i) n (2.165) ∂Q ∂φ (µ k+1 1 ) = 0 ⇐⇒ µ k+1 1 = n i=1 γ (i) k x (i) n (2.166)
The maximization of F via the E-step and M-step, like in Example 6, is the historical version of the EM algorithm. However, other optimization techniques referred to as EM variants are possible. For instance, one could replace the maximization in equation (2.158) by just one step of an iterative optimization algorithm giving a set of parameters θ k+1 that verifies:

Q(θ k+1 ) > Q(θ k ) (2.167)
When the M-step is performed in such an approximated way, the EM algorithm is renamed generalized EM [START_REF] Arthur P Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. Similarly, one could replace the exact E-step by just improving the current density estimate with respect to the value function q → F(q, θ) [START_REF] Harmeling | Multiframe blind deconvolution, super-resolution, and saturation correction via incremental EM[END_REF].

conclusion

In this chapter, we have introduced the maximum-likelihood estimator and have shown that it is consistent and asymptotically efficient. This estimator is the solution of an optimization problem. We have presented gradient descent and Newton's method, two popular iterative optimization methods that can be used to maximize the likelihood directly. Lastly, we explored the special case of latent variable models, where the maximum likelihood estimator can be estimated via the EM algorithm.

In the next chapter, we introduce the necessary background on neuro-imaging.

In the previous chapter, we have introduced some background on statistics and optimization. In this chapter, we present the necessary background on the data used in this thesis. In our experiments, we use functional magnetic resonance imaging (fMRI) or magnetoencephalograpy (MEG) data. As will be seen in this chapter, these two modalities have very different characteristics.

functional magnetic resonance imaging (fmri)

This section relies on the handbook of fMRI data analysis [START_REF] Russell A Poldrack | Handbook of functional MRI data analysis[END_REF].

The BOLD signal and hemodynamic response

The story of fMRI begins with the discovery of the blood-oxygenlevel-dependent (BOLD) imaging contrast in 1990 by Ogawa [START_REF] Ogawa | Brain magnetic resonance imaging with contrast dependent on blood oxygenation[END_REF]. When neurons fire, they consume oxygen and therefore, the level of oxygen in blood changes (it is actually over-compensated by blood supply). Since oxygenated and deoxygenated blood (oxyhemoglobin and deoxyhemoglobin) do not have the same magnetic susceptibility, changes in the oxygen level in blood can be tracked by a magnetic resonance imaging (MRI) scanner.

When a short stimuli occurs, the relative change in the MRI signal (BOLD signal) is not instantaneous. The typical response to a short stimuli follows the curve displayed in figure 3.1 and is called the hemodynamic response. The BOLD signal can be seen as the convolution of neural activation with the hemodynamic response.

Spatial and temporal resolution of fMRI data

The success of fMRI is due to its advantages compared to positron emission topography (PET) that uses radioactive tracers to follow glucose or water levels. PET and fMRI both have a spatial resolution of a few millimeters (which still encompasses several hundred thousands of neurons). However, while taking an image with PET takes about a minute, it takes on the order of 2s to produce with fMRI. In addition, fMRI is non-invasive and, unlike PET, it does not involve being exposed to radioactive tracers.

Experimental designs

In resting state fMRI, subjects are instructed to lie still in the scanner without further instruction. This kind of data can be used to extract networks that highlight regions that tend to co-activate. This contrasts with task fMRI, in which subjects have to comply with specific instructions. In task fMRI, two types of experimental designs are widespread. Block designs are suited to functional imaging modalities with low temporal resolutions. In block designs, subjects are exposed to a continuous stimulus that lasts of a rather long period of time before it switches to another. This contrasts with event related designs that are only suited to modalities with a finer time resolution. In event related designs, a sequence of short events are presented and separated by a time window (inter-stimulus interval). When the inter-stimulus interval is shorter than the length of the hemodynamic response, we talk about "fast" even related designs. While block designs or event related designs are controlled designs, a rising trend is to use more naturalistic paradigms that are unconstrained from behavioral manipulations and thus, more ecological with respect to real-world conditions. In such naturalistic paradigms, referred to as naturalistic task fMRI in this thesis, subjects are exposed to naturalistic stimuli such as movie watching or audio track listening. These kinds of stimuli constitute a middle ground between resting state and task fMRI data as activations are time-locked but the environment is not as controlled as in a typical task fMRI setting. In this thesis, we use mostly naturalistic task fMRI, as one of our motivations is to propose a suitable framework to analyse such data.

Preprocessing

The fMRI signal is particularly noisy. Several techniques are used to reduce the noise and enhance the quality of the data.

Distortion correction

In the MRI scanner, a constant magnetic field is applied. However, in areas where there is an air/tissue interface, the magnetic field is distorted leading to errors in the localization of voxels near the sinus or ears. However, most scanners come with a map quantifying the distance that each voxel has been shifted. By inverting the map, one can try and recover the correct location of voxels.

Slice timing correction

Images are acquired in slices, a few slices at a time, so that different voxels are acquired at different times. In order to correct for this effect, a reference slice is chosen and other slices are interpolated so that we can consider that all voxels are acquired at the same time.

Motion correction

Head motion causes massive distortion of the signal. Such effects are corrected by chosing a reference image and applying a rigid body transformations so that other images match the orientation and localization of the reference image. One also records estimated head movements parameters so that they can be regressed out. However, if the task is correlated with head motion, this regression can lead to a loss of information.

Spatial smoothing

Spatial smoothing blurs the image through the application of a Gaussian kernel with a given width (in mm). When the signal of interest has a large spatial extend, smoothing increases the signal to noise ratio. Smoothing is also used as a way to decrease between-subject variability. In our experiments, we usually don't apply smoothing in order to measure how our methods are able to handle fined grained details. While there is currently no consensus on smoothing, we observe that a smoothing of 3 mm generally improves most analysis on fMRI data.

Frequency filtering

Heating of the scanner causes a low temporal frequency noise to appear. In order to remove such noise, a high-pass filtering is applied with a cut-off frequency of 0.01Hz. Sometimes a low-pass filter with a cut-off frequency of 0.1 Hz is also used since artifacts induced by motion are usually of higher frequency than the signal of interest.

Spatial normalization

Each subject has a brain of different size and shape. The goal of spatial normalization is to align brain images of different subjects in order to reduce anatomical variability.

The most widespread technique is to register all images into the Montreal Neurological Institute template (MNI template). The MNI template is built from 305 anatomical images that are aligned based on anatomical features via a non-linear registration model. The nonlinear registration model uses rigid body transformations followed by non-linear diffeomorphic deformations to better match brain shape.

The projection of fMRI data to the MNI space proceeds in two steps. First, a high dimensional anatomical image of each subject is aligned on the MNI template. Then the fMRI images of each subjects are aligned on the anatomical image of the same subject (this step is called co-registration). The fMRI images are mapped to the MNI template by composing the two transformations.

Masking

It is often the case that only a subpart of the brain is of interest. When this is the case, we use a mask to only keep the set of voxels corresponding to particular locations defined by the mask. In the case of our experiments, we use a gray-matter mask, since this is the only part that matters when one wants to capture functional activations in the brain that reflect behavioral responses.

Runs

For the comfort of participants and to ensure their active participation, it is in general not possible to use the scanner without interruptions for more than 15 minutes. When long acquisitions must be performed, they are split into short runs of approximately 15 minutes.

Example of fMRI Datasets

In this subsection we provide examples of fMRI datasets. These datasets will be used to evaluate the methods developed in this thesis. Datasets are preprocessed with FSL http://fsl.fmrib.ox. ac.uk/fsl and SPM https://www.fil.ion.ucl.ac.uk/spm/software using slice timing correction, distortion correction spatial realignment, co-registration to the T1 image and affine transformation of the functional volumes to a template brain (MNI). Using Nilearn [START_REF] Abraham | Machine learning for neuroimaging with scikit-learn[END_REF], preprocessed data are masked (using a full brain mask available at http://cogspaces.github.io/assets/data/hcp _ mask.nii.gz), detrended and standardized (so that any voxel's timecourse has zero mean and unit variance). We also apply a high-pass filter and a lowpass filter with cut-off frequencies of 0.01 Hz and 0.1 Hz respectively).

Sherlock

In the SHERLOCK dataset, 17 participants are watching "Sherlock" BBC TV show (episode 1). These data are downloaded from http: //arks.princeton.edu/ark:/88435/dsp01nz8062179. Data were acquired using a 3T scanner with an isotropic spatial resolution of 3 mm. More information including the preprocessing pipeline is available in [START_REF] Chen | Shared experience, shared memory: a common structure for brain activity during naturalistic recall[END_REF]. Subject 5 is removed because of missing data, leaving us with 16 participants. Although SHERLOCK data contains originally only 1 run, we split it into 4 blocks of 395 timeframes and one block of 396 timeframes for the needs of our experiments.

Forrest

In the FORREST dataset, 20 participants are listening to an audio version of the movie Forrest Gump. FORREST data are downloaded from OpenNeuro [START_REF] Russell A Poldrack | Toward open sharing of task-based fMRI data: the OpenfMRI project[END_REF]. Data were acquired using a 7T scanner with an isotropic spatial resolution of 1 mm (see more details in [START_REF] Hanke | A high-resolution 7-Tesla fMRI dataset from complex natural stimulation with an audio movie[END_REF]. More information about the forrest project can be found at http://studyforrest.org. Subject 10 and run 8 are discarded because of missing data. We therefore use full brain data of 19 subjects split in 7 runs of respectively 451, 441, 438, 488, 462, 439 and 542 timeframes.

CamCAN

In the fMRI CamCAN dataset, 647 participants aged from 18 to 88 years are watching Alfred Hitchcock's "Bang! You're Dead" (edited so that it lasts only 8 minutes). CamCAN consists of data obtained from the CamCAN repository (available at http://www.mrc-cbu.cam.ac.uk/ datasets/camcan/) (see [START_REF] Taylor | The Cambridge Centre for Ageing and Neuroscience (Cam-CAN) data repository: structural and functional MRI, MEG, and cognitive data from a cross-sectional adult lifespan sample[END_REF] and [START_REF] Meredith A Shafto | The Cambridge Centre for Ageing and Neuroscience (Cam-CAN) study protocol: a cross-sectional, lifespan, multidisciplinary examination of healthy cognitive ageing[END_REF]). We use all available subjects and runs yielding 647 participants and 1 run of 193 timeframes.

Raiders

The RAIDERS dataset reproduces the protocol described in [START_REF] James V Haxby | A common, highdimensional model of the representational space in human ventral temporal cortex[END_REF]. 10 participants are watching the movie "Raiders of the Lost Ark". This dataset pertains to the Individual Brain Charting dataset [START_REF] Luísa Pinho | Individual Brain Charting dataset extension, second release of high-resolution fMRI data for cognitive mapping[END_REF][START_REF] Luísa Pinho | Individual Brain Charting, a high-resolution fMRI dataset for cognitive mapping[END_REF]. We use full brain data of 10 subjects split in 9 runs of respectively 374, 297, 314, 379, 347, 346, 350, 353 and 211 timeframes. Note that the Raiders dataset is different from the one used in [START_REF] Chen | A reduced-dimension fMRI shared response model[END_REF], as it involves different subjects, and because data were acquired at NeuroSpin using a 3T scanner with an isotropic spatial resolution of 1. dataset [START_REF] Luísa Pinho | Individual Brain Charting dataset extension, second release of high-resolution fMRI data for cognitive mapping[END_REF][START_REF] Luísa Pinho | Individual Brain Charting, a high-resolution fMRI dataset for cognitive mapping[END_REF] is an extension of the raiders dataset where the first two scenes of the movie are shown twice (130 mins).

CLIPS

The CLIPS dataset reproduces the protocol of original studies described in [START_REF] Nishimoto | Reconstructing visual experiences from brain activity evoked by natural movies[END_REF] and [START_REF] Alexander G Huth | A continuous semantic space describes the representation of thousands of object and action categories across the human brain[END_REF]. 10 participants are exposed to short clips. The data were acquired in 17 runs of 325 timeframes. The CLIPS dataset also pertains to the Individual Brain Charting dataset ( [START_REF] Luísa Pinho | Individual Brain Charting dataset extension, second release of high-resolution fMRI data for cognitive mapping[END_REF][START_REF] Luísa Pinho | Individual Brain Charting, a high-resolution fMRI dataset for cognitive mapping[END_REF]). The CLIPS and RAIDERS data are available in OpenNeuro under the identification number: ds002685. Protocols on the visual stimuli presented are available in a dedicated repository on Github:

https://github.com/hbp-brain-charting/public _ protocols.

Unless stated otherwise we use spatially unsmoothed data, except for the sherlock dataset, for which the available data are already preprocessed with a 6 mm spatial smoothing. The temporal resolution or repetition time (TR) is 2s for all datasets except for the Sherlock dataset where the TR is 1.5s.

A summary about the size of each dataset is available in Table 3.1. Note in particular that all datasets have been resampled to 2mm isotropic resolution, leading to 212,445 voxels in the brain mask.

An example of univariate analysis: analyzing task fMRI data via the general linear model

In this section, we describe a framework successfully applied to identify brain regions involved in specific tasks: the general linear model (GLM) [START_REF] Karl J Friston | Analysis of fMRI time-series revisited[END_REF]. We refer the reader to Poline [START_REF] Poline | The general linear model and fMRI: does love last forever?[END_REF] for a more detailed description of the model and only cover here what we consider to be the most important parts. We associate to each image t a set of numbers s(t) that describe the experiment and nuisance parameters. For instance, in an experiment where the subject is either asked to listening a sentence or reading a sentence we could have the association: s(t) = (φ(t), τ(t), x(t), y(t), z(t)) where φ(t) = 1 if the subject is reading a sentence and 0 otherwise, τ(t) = 1 if the subject is listening to a sentence and 0 otherwise and x(t), y(t), z(t) describe the position of the head. While φ and τ describe the experiment, x, y, z describe nuisance parameters. Note that in practice we use many more nuisance regressors, such as parameters describing the rotation of the head, heartbeats, respiration rhythm and really any other parameter that might introduce artifacts. We also usually convolve the regressors describing the experiment with the hemodynamic function so that they are closer to the actual brain response and sometimes also its derivative to account for small temporal deviations from our model of hemodynamic response. We call S ∈ R k,n the design matrix such that column t of S is given by s(t). Then, the general linear model sees the data of one subject X ∈ R v×n as a linear combination of the design matrix S with additive noise N.

X = AS + N (3.1)
where A ∈ R v×k contains k spatial maps in the sense that column j of A can be plotted as a brain image representing the localization of activity described by the row j of S. The noise is often assumed to be Gaussian with unit variance so that A can be obtained as the result of a linear regression. Another common hypothesis is to assume temporal correlation of the noise in the form of an AR-1 process. In this case, we consider instead the data XB where B ∈ R n×n is chosen so that the noise BN has no temporal correlation allowing us to perform linear regression.

The result of the GLM procedure is often given by contrasts maps where the spatial maps corresponding to conditions of interest are subtracted. In our reading versus listening example, we would typically display the first column of A (that corresponds to φ which describes the activation related to reading a sentence) minus the second column of A (that corresponds to τ which describes the activation related to listening to a sentence). An example of "reading versus listening contrast is displayed in Figure 3.2.

In this thesis, we always work with unthresholded contrast maps like the one in Figure 3.2. For interpretation purposes, contrast maps are usually z-scored and thresholded to keep only activity that is significantly different from zero. The GLM model is called univariate since interactions between voxels are not modeled. Univariate methods are often criticized for their inability to capture well correlations and interactions between brain-wide measurements. This makes it difficult to precisely locate brain functions from brain maps. An approach to overcome this issue is to train classifiers to decode brain maps, i.e to discriminate between different stimulus of task types [95, 138, Contrast "sentence listening versus sentence reading" computed from the fMRI data of one subject in the IBC dataset [START_REF] Luísa Pinho | Individual Brain Charting dataset extension, second release of high-resolution fMRI data for cognitive mapping[END_REF][START_REF] Luísa Pinho | Individual Brain Charting, a high-resolution fMRI dataset for cognitive mapping[END_REF]. This contrast was downloaded from neurovault [START_REF] Krzysztof | NeuroVault. org: a web-based repository for collecting and sharing unthresholded statistical maps of the human brain[END_REF] (collection 2138 subject 1 session 0) .

151].

When linear classifiers are used, the weights of this classifiers localize the brain functions. This approach can be more powerful, as classifiers take correlations between voxels into account. Decoding is also popular for individual imaging-based diagnosis.

magneto electro encephalography (meg)

Most material in this section is inspired from the book of Riitta Hari and Aina Puce [START_REF] Hari | MEG-EEG Primer[END_REF].

Principle of MEG

When a neuron fires, it emits a current that generates a magnetic field. By recording the magnetic field close to the skull, we gain insight on neural activity.

A limiting factor is that the magnetic fields induced by neuronal currents are extremely weak (on the order of 10fT ) which is much lower than the ambient magnetic noise (10 8 fT ). Therefore recordings are performed in a shielded room and extremely sensitive magnetometers are used. The best current tools can measure the magnetic field generated by approximately 50 000 neurons oriented in the same direction.

MEG recording device also include gradiometers in addition to magnetometers. These are sensitive to the local variations of the magnetic field. However in our experiments, we only use the magnetometers data.

MEG has a temporal resolution on the order of the ms which is essentially the same as electroencephalography (EEG). This is much better than fMRI. The spatial resolution is similar to that of EEG and is on the order of the centimeter. In general the brain sources are slightly better localized in MEG since the magnetic field is not affected by changes in conductivity in the head.

Solving the inverse problem

In order to locate activation from the brain measurements, we need to solve an inverse problem: what kind of sources can generate the magnetic field we observe ? This section gives an overview of the sLORETA method [START_REF] Domingo | Standardized lowresolution brain electromagnetic tomography (sLORETA): technical details[END_REF] and is strongly inspired from the tutorial of the authors [START_REF] Roberto D Pascual-Marqui | Discrete, 3D distributed, linear imaging methods of electric neuronal activity. Part 1: exact, zero error localization[END_REF].

From a vector describing the 3D current density inside the brain j ∈ R 3×v (v is the number of voxels) we can predict the magnetic field recorded by the sensors b ∈ R p (p is the number of sensors) by solving Maxwell's equations. This is called the forward model. Models describing the geometry of the head (head models) are used to provide an approximate solution. We have therefore

b = Kj (3.2)
where K ∈ R p×3v is the solution to the forward model. K can be seen as

K = [K 1 . . . K v ]
where K i ∈ R p×3 describes how voxel i contributes to the measured magnetic field b.

Note that if b and K are given there is an infinite number of possible j that satisfy equation (3.2), so the inverse problem must be regularized in order to select a solution among all possible ones. A common parametrization is to assume that the current at voxel i, j i ∈ R 3 verifies:

j i = (K i CK i ) -1 2 K i Cb (3.3)
for a symmetric matrix C ∈ R p×p . In sLORETA, we have C = (KK + αI) † where † represents Moore's pseudo inverse, I is the identity matrix and α is an hyperparameter. Taking a point wise source

a ∈ R 3 at voxel i * equation (3.2) yields b = K i * a. Then, j i = (K i CK i )K i * a is such that j i
2 is maximum at i = i * showing that the method can correctly recover point sources.

Preprocessing steps

This section is inspired from the preprocessing recommandations in [START_REF] Jas | A reproducible MEG/EEG group study with the MNE software: recommendations, quality assessments, and good practices[END_REF] and the book of Riitta Hari and Aina Puce [START_REF] Hari | MEG-EEG Primer[END_REF].

Temporal filtering

Temporal filtering is performed in two steps. First an analogous pass band filter selects a large band of frequencies between 0.01Hz and 200Hz. Then, a digital filtering is applied. Typically, a low-pass filtering with cut-off frequency at 40Hz is used as it is removes the line frequency at 50Hz while keeping most of the brain's signals.

Independent component analysis

Independent component analysis (ICA) can be used to isolate artifacts due to heart beating or muscle contraction. As will be seen in section 4.1. ICA extracts independent components from the data. The data can be cleaned by removing components corresponding to artifacts [START_REF] Jung | Extended ICA removes artifacts from electroencephalographic recordings[END_REF].

Maxwell filtering

Maxwell filtering also called signal space separation (SSS) [START_REF] Taulu | Spatiotemporal signal space separation method for rejecting nearby interference in MEG measurements[END_REF] identifies the contributions to the magnetic field from brain sources outside the brain and removes it. This decreases the effects of outside noise and even allows to scan subjects with implanted stimulators.

Experimental designs

In magneto-encephalography, we measure event related fields (ERF) related to an evoked response in the brain. Experiments are often repeated several times yielding a number of trials. By averaging the event related fields over trials, the signal to noise ratio increases.

MEG datasets

In this section, we give examples of MEG datasets. These examples are used in the rest of the thesis to evaluate the methods we develop.

The Sinusoidal Phantom MEG dataset uses data collected with a realistic head phantom, which is a plastic device mimicking real electrical brain components. Eight current dipoles positioned at different locations can be switched on or off. We only consider the 102 magnetometers. An epoch corresponds to 3 s of MEG signals where a dipole is switched on for 0.4 s with an oscillation at 20 Hz and a peak-topeak amplitude of 200 nAm. We have access to 100 epochs per dipole. Maxwell filtering is applied on the data as well as low-pass filtering with a cut-off frequency at 40 Hz.

The Elekta Phantom MEG dataset also uses data collected with a realistic head phantom and is available as part of the Brainstorm application [START_REF] Tadel | Brainstorm: a user-friendly application for MEG/EEG analysis[END_REF]. This dataset uses 32 dipoles positioned at different locations. Like in the Sinusoidal Phantom MEG dataset, we only consider the 102 magnetometers and apply Maxwell filtering and low-pass filtering with a cut-off frequency at 40 Hz. The dipoles emit a signal at either a strong, medium or low amplitude yielding 3 different datasets. We use the dataset where the emitted signal is very strong to recover the true signal by performing a PCA with 1 component. Then we work with the dataset where the emitted signal is the weakest. Each epoch corresponds to 301 samples and 20 epochs are available in total.

The CamCAN dataset [START_REF] Meredith A Shafto | The Cambridge Centre for Ageing and Neuroscience (Cam-CAN) study protocol: a cross-sectional, lifespan, multidisciplinary examination of healthy cognitive ageing[END_REF] contains the MEG data of 647 different subjects exposed to an audio-visual stimuli. More precisely, subjects are presented simultaneously an auditory stimuli lasting 300ms at frequency 300, 600 or 1200 Hz and a checkerboard pattern lasting 34ms. 120 trials are available.

conclusion

In this chapter, we described the principles on which fMRI and MEG imaging are based as well as the standard preprocessing pipelines and the datasets used in this thesis.

In the next chapter, we review several methods to perform unsupervised analysis of neuroimaging data.

R E V I E W O F S E L E C T E D U N S U P E R V I S E D M E T H O D S P O P U L A R I N N E U R O I M A G I N G S T U D I E S

When exposed to naturalistic stimuli (e.g. movie watching or simulated driving), subjects' experience is closer to their every-day life than with classical psychological experiments. This makes naturalistic paradigms an attractive class of stimulation protocols for brain imaging. However, such stimulations are difficult to model, therefore the statistical analysis of the data using supervised regression-based approaches is challenging. This has motivated the use of unsupervised learning methods that do not make assumptions about what triggers brain activations in the presented stimuli.

In this chapter, we first present independent component analysis (ICA), a widely used unsupervised method for neuroimaging studies routinely applied on individual subject electroencephalography (EEG) [START_REF] Makeig | Independent component analysis of electroencephalographic data[END_REF], magnetoencephalography (MEG) [START_REF] Vigário | Independent component analysis for identification of artifacts in magnetoencephalographic recordings[END_REF] or functional MRI (fMRI) [START_REF] Martin | Independent component analysis of fMRI data: examining the assumptions[END_REF] data. Then, we review multiview unsupervised techniques that leverage the availability of data from multiple subjects performing the same experiments.

independent component analysis

Independent component analysis (ICA) models a set of signals as the product of a mixing matrix and a component matrix containing independent components. As will be seen in this section, the required assumptions on the independent components to guarantee identifiability are rather weak, making ICA a method of choice to analyze the data of subjects exposed to a stimulus that is difficult to quantify.

ICA is applied to fMRI data to analyze resting state data [START_REF] Christian F Beckmann | Investigations into resting-state connectivity using independent component analysis[END_REF] or when subjects are exposed to natural [98] [14] or complex stimuli [START_REF] Vince D Calhoun | Different activation dynamics in multiple neural systems during simulated driving[END_REF]. In M/EEG processing, it is widely used to isolate acquisitions artifacts from neural signal [START_REF] Jung | Extended ICA removes artifacts from electroencephalographic recordings[END_REF], and to identify brain components of interest [START_REF] Delorme | Independent EEG sources are dipolar[END_REF][START_REF] Vigário | Independent component approach to the analysis of EEG and MEG recordings[END_REF].

Mainly for computational reasons, it is often assumed that the number of components k is much lower than the dimensionality of the data v. However in ICA, the dimensionality of the data must be equal to the number of components. We therefore first present principal component analysis, a standard method to perform dimension reduction, then present non-Gaussian ICA and non-stationary ICA.

Principal component analysis (PCA)

Let us assume our data are given by n observations of a random vector x ∈ R v that we stack into a matrix X ∈ R v×n . The data are centered (E[x] = 0). Principal component analysis (PCA) yields an orthonormal familly of p vectors R ∈ R v×p , R R = I p , such that the projected data RR x does not yield a large mean reconstruction error in Frobenius norm. The corresponding optimization problem is given by: argmin

R,R R=I E x x -RR x 2 (4.1) = argmin R,R R=I E x trace(x (I -RR )(I -RR )x) (4.2) = argmin R,R R=I E x ( x 2 -R x 2 ) (4.3) = argmax R,R R=I E x R x 2 (4.4) = argmax R,R R=I trace(R V[x]R) (4.5)
Therefore the matrix R is just given by the first p eigenvectors of V[x].

In practice, we can use the sample variance 1 n XX to estimate V[x]. However, such an estimate is not practical in high dimension: when v is very large XX is a prohibitively large matrix.

In such case, assuming n is small, we rely instead of the singular value decomposition (SVD) of X:

X = UDV (4.6)
where U ∈ R v×n is an orthogonal matrix (U U = I t ) of left-singular vectors, D ∈ R n is a positive diagonal matrix of singular values and V ∈ R n×n is an orthogonal matrix of right-singular vectors. Such a decomposition can be computed in Õ(vn 2 ) operations which is not prohibitive when n is small enough. Every matrix has a singular value decomposition and, provided that all singular values are distinct, this decomposition is unique up to a permutation and sign indeterminacy. More precisely, if UDV and U D V are two singular decomposition of X and if D has only distinct values, then U = UΠΞ, D = Π DΠ and V = ΞΠ V where Π is a permutation matrix and Ξ a diagonal matrix with diagonal values in {-1, 1}.

Let X = UDV be a singular value decomposition of X, we have XX = UD 2 U and therefore the left-singular vectors U p corresponding to the largest p singular values of X yield the p largest eigenvectors of XX . Therefore R is given by the left-singular vectors corresponding to the largest p singular values: R = U p .

Sometimes the PCA includes whitening R = U p D -1 p where D p contains the p largest singular values of X so that the components of R x are uncorrelated (but R is no longer orthonormal). In this thesis what we call PCA does not include signal whitening.

After the data are reduced, we can perform independent component analysis (ICA). ICA can exploit a number of different properties of the signal. We focus in the next section on non-Gaussian ICA.

Non Gaussian ICA

ICA sees the data x ∈ R p as a linear mixtures of components s ∈ R p . Data are therefore modeled as:

x = As (4.7)
where A is the unmixing matrix and s is a vector with independent components, of which at most one is Gaussian and whose densities are not reduced to a point-like mass. Without loss of generality, the data x are assumed centered (E[x] = 0). This section relies heavily on the ICA review [START_REF] Hyvärinen | Independent component analysis: algorithms and applications[END_REF] and on [START_REF] Cardoso | Infomax and maximum likelihood for blind source separation[END_REF].

Identifiability

A matrix P that can be written P = ΞΠ, where Ξ is a diagonal matrix and Π is a permutation matrix, is called a scale and permutation matrix.

If furthermore Ξ only has diagonal values in {-1, 1}, then P is a sign and permutation matrix.

It is easily seen that if s is a vector with independent components, of which at most one is Gaussian and whose densities are not reduced to a point-like mass, so is Ps where P is a scale and permutation matrix. Therefore, if x = As, take A = AP and s = Ps and we have x = A s . In other words, there exists a permutation and scaling indeterminacy.

To fix the scaling indeterminacy, we assume that s has unit variance (V[s] = I). In addition, we assume that the data are whitened (V[x] = I). Note that the whitening assumption does not imply any loss of generality. If a matrix H is used to whiten the data, the mixing matrix of the unwhitened data is given by H -1 A where A is the unmixing matrix obtained on the whitened data. With these assumptions, A is orthogonal. Indeed,

V[x] = A A = I.
Are there any other indeterminacies ? Assume that there exists two mixing matrices A 1 and A 2 such that x = A 1 s 1 and x = A 2 s 2 with A 1 and A 2 orthogonal. Then s 1 = Os 2 where O = A 2 A 1 is an orthogonal matrix. The following theorem in [START_REF] Comon | Independent component analysis, a new concept?[END_REF] shows that O is necessarily a sign and permutation matrix. Theorem 8. Let s 2 be a vector with independent components, of which at most one is Gaussian, and whose densities are not reduced to a point-like mass. Let O be an orthogonal matrix and s such that s = Os 2 . Then, s has independent component if and only if O is a sign and permutation matrix.

As a result, if we assume that at most one component is Gaussian and that components densities are not reduced to a point-like mass, ICA is identifiable up to a scale and permutation indeterminacy.

Infomax: Maximum likelihood estimation

Infomax is introduced in [START_REF] Anthony | An informationmaximization approach to blind separation and blind deconvolution[END_REF] and although initially formulated as a maximum entropy problem, it has been shown in [START_REF] Cois | Infomax and maximum likelihood for blind source separation[END_REF] to be equivalent to maximum likelihood.

Assume for simplicity that the components have the same density δ, denoting f(x) =log(δ(x)) and f(x) = m i=1 f(x i ) the expected negative log-likelihood is given by:

L(W) = E[-l(x, θ)] = E[-log(p(x))] (4.8) = -log(|W|) -E[log(δ(Wx))] (4.9) = -log(|W|) -E[log(δ(y))] (4.10) = -log(|W|) + E[f(y)] (4.11) = -log(|W|) + E[ m i=1 f(y i )] (4.12)
where y = Wx and y i is the i-th component of y. At line (4.9) we used the change of variable s = Wx.

Optimizing the maximum likelihood: relative gradient, Hessian and approximations

The matrix W needs to be invertible. In practice, if the constraint is not enforced, we could see numerical instabilities appear. A simple rule that preserves the invertibility of W is to use updates of the form:

W ← (I + αD)W (4.13)
where α is a small step-size and D is a direction to be found. Following the same reasonning as in section 2.2.1.1, we assume a small step-size and write at first order: argmin

D, D =1 L((I + αD)W) = argmin D, D =1 L(W) + ∂L(W) ∂W , αDW (4.14)
= argmin

D, D =1 L(W) + ∂L(W) ∂W W , αD (4.15)
= -

∂L(W) ∂W W ∂L(W) ∂W W (4.16)
Therefore the steepest direction for that update rule is D = ∂L(W) ∂W W and therefore updates are given by: W ← (I -αG)W where G = ∂L(W) ∂W W is called the relative gradient [START_REF] Cardoso | Equivariant adaptive source separation[END_REF]. Second order extensions of the method are obtained by following the steps in section 2.2.1.2. The corresponding Hessian is called the relative Hessian.

Let us follow [START_REF] Ablin | Faster independent component analysis by preconditioning with Hessian approximations[END_REF] and use a quasi-Newton algorithm to minimize the likelihood. The relative gradient and Hessian of L are given by:

G = E[f (y)y ] -I p , (4.17)
where y = Wx and f (y) = ∂f(y) ∂y and

H abcd = δ ad δ bc + δ ac E[f (y a )y b y d ] (4.18)
Following [START_REF] Ablin | Faster independent component analysis by preconditioning with Hessian approximations[END_REF], we can approximate the Hessian by

Habcd = δ ad δ bc + δ ac δ bd Γ ab (4.19)
where

Γ ab = E[f (y a )y 2 b ].
The approximation is exact when the true unmixing matrix is found since E[s b s d ] = s 2 b δ bd . The updates are then given by:

W ← (I -α H-1 G)W (4.20)
The approximated Hessian H is block diagonal. Indeed we have for any matrix M:

( HM) ab ( HM) ba = Γ ab 1 1 Γ ba M ab M ba if a = b ( HM) aa = (1 + Γ aa )M aa (4.21)
So that each block corresponds to a pair (a, b) and is of size 2 if a = b and of size 1 when a = b. Therefore H can be easily regularized and inverted.

Robustness to density mismatch

In practice, we observe that sometimes, the quasi-Newton algorithm described in the previous section fails to recover the true mixing matrix. This happens when the sources used in the model are too far from the actual generating sources. In this section, we explain why this happens and derive stability conditions. This follows the work done in [START_REF] Cardoso | Blind signal separation: statistical principles[END_REF].

When the density of the sources corresponds to the one used in the model, one recovers the correct unmixing matrices via the maximum likelihood estimate in the limit of large samples. This comes from the consistency of maximum likelihood estimators. When this is not the case, a mismatch appears which can be quantified. Let us denote x * the true data. As highlighted in (2.54) the expected negative loglikelihood coincides, up to a constant, with the KL divergence between the true distribution of the data and the distribution of the data as hypothesized in the model:

L(W) (4.22) = D KL (p(x * ), p(W -1 s)) (4.23) = D KL (p(Wx * ), p(s)) (4.24) = D KL (p(Wx * ), i p i (w i x * )) + D KL ( i p i (w i x * ), p(s)) (4.25) = D KL (p(Wx * ), i p i (w i x * )) + D KL ( i p i (w i x * ), i p i (s i )) (4.26) = D KL (p(Wx * ), i p i (w i x * )) + i D KL (p i (w i x * ), p i (s i )) (4.27)
where equation (4.24) comes from the invariance of the KL divergence and we denote i p i (x i ) the product of the marginal densities of x.

The first term in equation (4.27) is the mutual information. It quantifies the independence of unmixed data Wx. The second term quantifies the mismatch between the assumed distribution of the components and the marginals of unmixed data. Therefore if the assumed components are too far from the true components, the recovered unmixing matrices may be far from the true ones.

Looking at the relative gradient G, we see that if components are truly independent, the true unmixing matrix will be a stationary point of the loss up to a scaling: G(ΛA -1 ) = 0 where Λ = diag(λ 1 . . . λ p ). However, in order for the quasi-Newton to reach this point, it must be a local minimum (the Hessian must be definite positive).

Let us therefore consider the Hessian H at point W = ΛA -1 where Λ is chosen such that G = 0. The unmixed data y = Wx = Λs are independent and therefore, we have H = H. As already mentioned, H is block diagonal with blocks described by equation (4.21). The condition for stability is that the Hessian is positive definite. Therefore all the blocks need to be positive definite so we get the conditions:

∀a, 1 + Γ aa > 0 for blocks of size 1 (4.28) ∀a = b, Γ ab Γ ba > 1
for blocks of size 2 (4.29)

When the conditions are satisfied, the quasi-Newton algorithms will recover the true unmixing matrices if initialized close to a solution. However, we have no theoretical guarantees of convergence because of the non-convexity of the problem.

Non-stationary ICA and joint diagonalization

Up to now we have assumed that samples were independent an identically distributed. In non-stationary ICA, samples are no longer identically distributed as the distribution can vary over time. This section follows the work in [START_REF] Pham | Blind separation of instantaneous mixtures of nonstationary sources[END_REF]. The components are assumed to be Gaussian with a variance that varies between samples but is assumed to be piece-wise constant Σ t = Σ k for t ∈ T k where (T k ) k is a partition of [1, 2, . . . , n] and Σ k is a positive diagonal matrix.

Denote X[i] to be sample i of the observed data X. The negative empirical expected log-likelihood is given by:

L = -log(|W|) + 1 n k 1 2 i∈T k Σ -1 2 k WX[i] 2 + 1 2 log(|Σ k |) (4.30) Denoting C k = i∈T k X[i]X[i] we have L = -log(|W|) + 1 n k 1 2 trace(Σ -1 k WC k W ) + 1 2 log(|Σ k |) (4.31)
Minimizing L with respect to Σ k yields Σ k = diag(WC k W ) and therefore up to a constant:

L = -log(|W|) + 1 n k 1 2 log(| diag(WC k W )|) (4.32)
which up to a constant can be rewritten:

L = 1 2n k log(| diag(WC k W )|) log(|WC k W |) (4.33)
L is a joint diagonalization criterion. The optimal W can be found via a quasi-Newton method very similar to the one we used for non-Gaussian ICA [START_REF] Ablin | Beyond Pham's algorithm for joint diagonalization[END_REF]. However, the updates only depend on the covariance matrices C k and no longer on the number of samples making this approach very fast when the number of samples is large.

Other approaches and extensions

non-white ica While non-stationary ICA relaxes the assumption that samples are identically distributed, non-white ICA relaxes the assumption that samples are independent. Following the work in [START_REF] Pham | Blind separation of instantaneous mixture of sources via the Gaussian mutual information criterion[END_REF], the entropy of a stationary Gaussian process y(t) is given by:

h(y(t)) = 1 4π π -π log(det(4π 2 eF(λ)))dλ (4.34)
where F is the spectral density matrix of the process. Then the authors define the Gaussian mutual information between stationary processes y 1 (t), . . . , y k (t) by:

I g (y 1 , . . . , y k ) = 1 4π π -π log(det(diag(F(λ)))) -log(det(F(λ)))dλ (4.35)
which is a joint diagonalization criterion similar to (4.33).

extended infomax: more flexible densities In Infomax, the densities are fixed once and for all. However, it can be useful to make them more flexible. Note that the updates in equation (4.20) and the stability criterion in equation (4.28) only depend on the score function f (y) = ∂f(y) ∂y . Choosing f (x) = x + tanh(x) would only separate super-Gaussian sources while choosing f (x) = xtanh(x) would only separate sub-Gaussian sources. In [START_REF] Lee | Independent component analysis using an extended infomax algorithm for mixed subgaussian and supergaussian sources[END_REF], the authors therefore propose an extended Infomax algorithm, where the score function can change based on observed statistics so that the stability criterion is fulfilled in any case.

ica via mutual information minimization In the limit of large samples, we have seen in equation (4.27) that maximum likelihood minimizes the sum of two terms: the mutual information and a second term that quantifies the mismatch between the assumed distribution of the components and the marginals of unmixed data. It might be desirable to only minimize the mutual information and forget the second term. Since the mutual information depends on the entropy, different works such as [START_REF] Hyvärinen | Analysis and projection pursuit[END_REF] [88] [START_REF] Li | Blind spatiotemporal separation of second and/or higher-order correlated sources by entropy rate minimization[END_REF] or [START_REF] Li | Independent component analysis by entropy bound minimization[END_REF] studied how the entropy can be estimated accurately and efficiently.

analysis of multiview data

In this section, we present multiview unsupervised techniques suited to analyze the data of multiple subjects exposed to the same complex stimuli. Such techniques assume some similarity between the data of different subjects. This assumption can be justified by the findings of [START_REF] Hasson | Intersubject synchronization of cortical activity during natural vision[END_REF] showing that brains exposed to the same natural stimuli exhibit synchronous activity. The task of finding common patterns or responses that are shared between subjects is called shared response modeling.

In the general linear model presented in section 3.1.6, the shared response is assumed to be known. Therefore, multiple subjects can be studied separately assuming that the data of different subjects are independent given the shared response. In the unsupervised setting it may not be so straightforward to deal with multiple subjects and therefore many different methods for data-driven multivariate analysis of neuroimaging group studies have been proposed. We summarize the characteristics of some of the most commonly used ones.

Multiset canonical correlation analysis

Canonical correlation analysis is initially designed to find a linear combination that maximizes the correlation between two datasets. The extension to more than two datasets is ambiguous, and many different generalized CCA methods have been proposed. [START_REF] Kettenring | Canonical analysis of several sets of variables[END_REF] introduces 6 objective functions that reduce to CCA when m = 2 and [107] considered 4 different possible constraints leading to 24 different formulations of Multiset CCA.

In [START_REF] Kettenring | Canonical analysis of several sets of variables[END_REF], the different formulations of multiset CCA are not derived from a probabilistic model. However, later works (see [START_REF] Li | Joint blind source separation by multiset canonical correlation analysis[END_REF] or [START_REF] Adali | Diversity in independent component and vector analyses: Identifiability, algorithms, and applications in medical imaging[END_REF]) have shown that some formulations of multiset CCA can be related to probabilistic models.

In this section, we present the formulation refered to in [START_REF] Aasbjerg | Multiset canonical correlations analysis and multispectral, truly multitemporal remote sensing data[END_REF] as "SUMCORR with constraint 4" which is one of the fastest to fit.

Let us consider X 1 , . . . , X m ∈ R p,n , m datasets and consider the following (SUMCORR) objective: max

a 1 ∈R v ,...,a m ∈R v m i=1 m j=1 a i , X i X j a j (4.36)
This objective can be arbitrarily large if not constrained. Constraint 4 is given by:

m i=1 a i , X i X i a i = 1 (4.37)
The Lagrangian is given by:

m i=1 m j=1 a i , X i X j a j -λ( m i=1 a i , X i X i a i -1) (4.38)
Taking the gradient with respect to a i we obtain

m j=1 X i X j a j = λX i X i a i (4.39)
This is a generalized eigenvalue problem of the form Ca = λDa where C is a block matrix where block i, j is given by X i X j , D is the block diagonal matrix formed by the block i, i of C and a ∈ R m×v yields the dataset specific projections vectors:

a =     a 1 • • • a m     .
The leading eigenvector correspond to the first canonical vectors. The second canonical vectors is given by the second eigenvalues and so on. They are orthogonal for the scalar product: a, b D = a, Db .

Group independent component analysis

Given the success of ICA in analyzing the data of one subject. It is natural to look for extensions of ICA in a multiview setting. Several works assume that the subjects share a common mixing matrix, but with different components [START_REF] Pfister | Robustifying independent component analysis by adjusting for group-wise stationary noise[END_REF] [START_REF] Svensén | ICA of fMRI group study data[END_REF]. Instead, we focus on models where the subjects share a common components matrix, but have different mixing matrices.

CanICA and ConcatICA

In the single subject setting, we reduce the data (for example using PCA) and apply ICA on reduced data. Therefore a natural framework to perform group ICA is to first aggregate the data of individual subjects into a single dataset, often resorting to dimension reduction technique and then apply off-the-shelf ICA on the aggregated dataset. When PCA is used to aggregate the data, the method is referred to as ConcatICA [START_REF] Vince D Calhoun | A method for making group inferences from functional MRI data using independent component analysis[END_REF]. An alternative is to use multiset canonical correlation analysis (CCA) leading to a method called CanICA [START_REF] Varoquaux | CanICA: Model-based extraction of reproducible group-level ICA patterns from fMRI time series[END_REF].

This framework has the advantage of being simple and straightforward to implement since it resorts to customary single-subject ICA method.

When datasets are high-dimensional, a three steps procedure is often used: first dimensionality reduction is performed on data of each subject separately; then the reduced data are merged into a common representation; finally, an ICA algorithm is applied for shared components extraction.

CanICA and ConcatICA are popular methods for fMRI [START_REF] Vince D Calhoun | A review of group ICA for fMRI data and ICA for joint inference of imaging, genetic, and ERP data[END_REF] and EEG [START_REF] Eichele | EEGIFT: group independent component analysis for event-related EEG data[END_REF] group studies. These methods directly recover only group level, shared components; when individual components are needed, a back-reconstruction step is required (such as GICA1 [START_REF] Vince D Calhoun | A method for making group inferences from functional MRI data using independent component analysis[END_REF] GICA2, GICA3 [START_REF] Barry Erhardt | Comparison of multi-subject ICA methods for analysis of fMRI data[END_REF] or dual-regression [START_REF] Christian F Beckmann | Group comparison of resting-state fMRI data using multi-subject ICA and dual regression[END_REF]).

In our experiments, when we fit CanICA or ConcatICA, we use Infomax and the picard solver [START_REF] Ablin | Faster independent component analysis by preconditioning with Hessian approximations[END_REF] with the tanh non-linearity. To obtain k components, we first apply a subject specific PCA with k components and then aggregate the data using either a Group PCA (in ConcatICA) or a multiset CCA (in CanICA) with k components. When individual sources are needed, we use dual-regression. As already noted, many different procedure to perform single subject ICA exists (such as Infomax, EBM [START_REF] Li | Blind spatiotemporal separation of second and/or higher-order correlated sources by entropy rate minimization[END_REF], ERBM [START_REF] Li | Independent component analysis by entropy bound minimization[END_REF]), the number of components in the first PCA can be chosen in many different ways and so can the back-reconstruction method. This leads to a number of different CanICA / ConcatICA algorithms. Many of these different versions of CanICA / ConcatICA are implemented in the GIFT toolbox https://trendscenter.org/software/gift/. We highlight that GIFT uses a slightly different terminology as us: dual-regression is referred to as spatial-temporal regression and ConcatICA is simply referred to as GroupICA. In this work, GroupICA describes instead the problem of recovering common sources from multiple subjects.

In the rest of the thesis, we often make the point that CanICA and ConcatICA are not maximum likelihood estimators. This is because the group PCA / multiset CCA used to fuse the data of different subjects does not come from a maximum likelihood approach.

Likelihood based method

While CanICA and ConcatICA are simple to implement and very fast to fit, they do not rely on maximum likelihood estimators. Therefore they do not benefit of advantages of such estimators such as asymptotic efficiency.

The model of [START_REF] Guo | A unified framework for group independent component analysis for multi-subject fMRI data[END_REF] considers the very general model x i = A i s + n i where A i represent mixing matrices, s the common sources and n i are view specific noise. Therefore this model is an instance of noisy ICA as defined in [START_REF] Hyvarinen | Gaussian moments for noisy independent component analysis[END_REF]. The noise covariance is learned from the data and each source is assumed to be a mixture of Gaussians p(s j ) = q z=1 N(µ zj , σ zj ). The parameters of the Gaussian mixtures are learned which makes the E-step impossible to compute in closed form. In order to solve this issue, an approximate E-step is introduced. Unfortunately, it leads to an update rule involving a sum over q p terms making their algorithm intractable when the number of components p is larger than 20.

In this section, we show why a sum of an exponential number of terms appears. We start with the same model as in [START_REF] Guo | A unified framework for group independent component analysis for multi-subject fMRI data[END_REF] but to make the computations more tractable, we assume that the Gaussian mixture is given by: p(s j ) = 1 q

α j ∈A p(s j |α j ) (4.40) p(s j |α j ) = N(s j ; 0, α j ), (4.41)
where α j takes its value in a known discrete set A with equal probability 1 q where q is the cardinal of A. We call α the random vector with independent coordinates such that coordinate j is given by α j . We further assume that the noise distribution is the same for all components and all subjects leading to the formulation:

x = As + n (4.42)
where

A =     A 1 . . . A m     , x =     x 1 . . . x m     and n =     n 1 . . . n m     with n ∼ N(0, σ 2 I).
This formulation is a special case of [START_REF] Guo | A unified framework for group independent component analysis for multi-subject fMRI data[END_REF] and constitues a single subject noisy ICA problem almost identical to [START_REF] Moulines | Maximum likelihood for blind separation and deconvolution of noisy signals using mixture models[END_REF].

We now follow [START_REF] Moulines | Maximum likelihood for blind separation and deconvolution of noisy signals using mixture models[END_REF] and write:

p(x, s, α) (4.43) = p(x|s)p(s|α)p(α) (4.44) = N(x; As, σ 2 I p )N(s; 0, diag(α)) 1 2 p (4.45) ∝ exp -1 2 ( 1 σ 2 x -As 2 + s, diag(α) -1 s ) (|σ 2 I p || diag(α)|) 1 2 (4.46) = exp -1 2 ( 1 σ 2 ( x 2 -2 x, As + As 2 ) + s, diag(α) -1 s ) (|σ 2 I p || diag(α)|) 1 2 (4.47) ∝ exp -1 2 ( s -µ α , V -1 α (s -µ α ) -µ α , V -1 α µ α ) (|σ 2 I p || diag(α)|) 1 2 , ( 4.48)
where

V α = ( 1 σ 2 A A + diag(α) -1 ) -1 , µ α = 1 σ 2 V α A
x and the proportionality constant contains terms that do not depend on s or α.

From p(x, s, α), we get:

p(s|x, α) = N(s, µ α , Σ α ) (4.49)
Then, we have that:

p(α|x) = s p(α, s|x)ds (4.50) ∝ s p(x, s, α)ds (4.51) ∝ exp 1 2 ( µ α , V -1 α µ α ) (| diag(α)||V -1 α |) 1 2 (4.52) (4.53)
where we leave out terms that do not depend on α. The normalizing constant can be computed by summing over possible values of α:

p(α|x) = exp 1 2 ( µ α ,V -1 α µ α ) (| diag(α)||V -1 α |) 1 2 α,α j ∈A exp 1 2 ( µ α ,V -1 α µ α ) (| diag(α)||V -1 α |) 1 2 (4.54)
Then, we can obtain a formula in closed form for p(s|x) using

p(s|x) = α,α j ∈A p(s|x, α)p(α|x) (4.55)
The problem here is that the size of the set {α, α j ∈ A} is q p where q = |A| is the cardinal of A. This quantity quickly gets large when p increases making p(s|x) difficult to compute.

In this section, we have studied a simplified version of the model in [START_REF] Guo | A unified framework for group independent component analysis for multi-subject fMRI data[END_REF] and shown that it is difficult to fit. This is often the case with maximum likelihood approaches. Despite, their advantages they are often intractable.

Independent vector analysis

Independent vector analysis [START_REF] Lee | Independent vector analysis (IVA): multivariate approach for fMRI group study[END_REF] (IVA) models the data as a linear mixture of independent components x i = A i s i , where each component s ij of a given view i can depend on the corresponding component in other views:

s [j] = (s ij ) m i=1 are not independent. Introducing x ∈ R m×p such that x = [x 1 , . . . , x m ] , s ∈ R m×p such that s = [s 1 , . . . , s m ] and y ∈ R m×p such that y = [y 1 , . . . , y m ] where y i = W i x i with W i = A -1
i , the expected negative loglikelihood is given by:

L = -E[log(p(x))] = m i=1 -log(|W i |) -E[log(p(y))] = m i=1 -log(|W i |) + p j=1 -E[log(p s [j] (y [j] ))]
where we used the notation y [j] = (y ij ) m i=1 . The optimization can be carried out using alternate minimization keeping the mixing matrices of all subjects fixed but one. We can rely on the relative gradient as in section 4.1.2.3 and use update of the form W i ← (I -α i G i )W i where α i is given by backtracking line search and G i is the relative gradient given by:

G i = -I p + E[φ i (y i )y i ] (4.56)
where component j of φ i (y i ) is given by

φ ij (y i ) = ∂ -log(p s [j] (y [j] )) ∂y ij (4.57)
Practical implementations of this general model assume a distribution for p s [j] . In IVA-L [START_REF] Lee | Independent vector analysis (IVA): multivariate approach for fMRI group study[END_REF],

p s [j] (y [j] ) ∝ exp(- i (y ij ) 2 ) (4.58)
and therefore

φ ij (y i ) = y ij i (y ij ) 2 (4.59) In IVA-G [7] [152], p s [j] (y [j] ) = N(y [j] ; 0, Σ j ) (4.60)
and therefore

φ ij (y i ) = l Σ -1 j [il]y lj (4.61)
where Σ -1 j [il] is the coordinate i, l of Σ -1 j and y lj is the jth coordinate of y i .

In IVA-G, an estimate of Σ j is needed at each iteration. This is computed using the sample covariance:

Σ j = 1 n Y [j] Y [j] (4.62)
Second order extensions and Hessian approximations can be used in IVA as well. This is described in [START_REF] Anderson | Joint blind source separation with multivariate Gaussian model: Algorithms and performance analysis[END_REF]. Also note that although IVA-G and IVA-L are the two most popular implementations of the IVA framework, others exist (see for instance the work in [START_REF] Anderson | Independent vector analysis, the Kotz distribution, and performance bounds[END_REF]).

Hyperalignment

Hyperalignment is a model initially designed for fMRI data to reduce inter-subject variability [START_REF] James V Haxby | A common, highdimensional model of the representational space in human ventral temporal cortex[END_REF].

Let us assume we have access to the data of two subjects: x 1 , x 2 . Assuming these subjects are exposed to a time-locked stimuli (such as a movie), a possible alignment is given by the Procrustes transform:

min P∈R p×p ,PP =I p E[ Px 1 -x 2 2 ] (4.63)
This can be solved efficiently by

P = P(E[x 2 x 1 ]) (4.64)
where P is the projection on the orthogonal manifold:

P(M) = M(M M) -1 2 . In practice P(M) is computed by performing an SVD of M, M = U M D M V M so that P(M) = U M V M .
Hyperalignment is the combination of the Procrustes transform and an iterative procedure to produce a template from multiple alignments. In an initialization step, a random subject i is chosen and the alignment between all subjects s = i and the target are computed. The initial template t is given by the averaged aligned data. Then, all subjects are aligned to the current template t and the template is recomputed using the averaged aligned data. This procedure is repeated for a given number of iterations until convergence.

The intuition behind the iterative procedure is that averaging the aligned data will tend to move the template away from the initial target. However, they are no theoretical guarantees associated with this procedure and even no associated loss. We describe formally the method in Algorithm 1. The shared response model [START_REF] Chen | A reduced-dimension fMRI shared response model[END_REF] is a multi-view latent factor model. The data x 1 . . . x m are modeled as random vectors following the model:

x i = A i s + n i (4.65) A i A i = I p (4.66)
where x i ∈ R v is the data of view i, A i ∈ R p×v is the mixing matrix of view i, n i is the noise of view i and s ∈ R p are the shared components referred to as the shared response in fMRI applications. The mixing matrices A i are assumed to be orthogonal so that A i A i = I p . However, in general the matrix A i A i is different from identity. The noise n i is assumed to be Gaussian with covariance Σ i and independent across views. We assume the number of features v to be much larger than the number of components p: v p. The conceptual figure 4.1 illustrates an application of the shared response model to fMRI data. The mixing matrices are spatial topogra- phies specific to each subjects while the shared components give the common timecourses.

In [START_REF] Michael J Anderson | Enabling factor analysis on thousand-subject neuroimaging datasets[END_REF][START_REF] Chen | A reduced-dimension fMRI shared response model[END_REF], two versions of the shared response model are introduced which we now present.

Deterministic shared response model

Let us consider n observations of x i and s that we stack into matrices X i ∈ R v,n and S ∈ R p,n . The deterministic shared response model sees both the mixing matrices A i and the n observations of the shared response S as parameters to be estimated. The noise variance is fixed to a multiple of identity: ∀i, Σ i = σ 2 I v where σ is an hyper-parameter to choose. The model is optimized by maximizing the log-likelihood. The likelihood is given by: p(x) = i N(x i ; A i s, σ 2 I) and therefore the empirical expected negative log-likelihood is given up to a constant independent of A i and S by:

L = 1 n i A i S -X i 2 = 1 n S 2 -2 A i S, X i + X i 2 (4.67)
The negative log-likelihood L is optimized by performing alternate minimization on (A 1 . . . A m ) and S. Note that the hyper-parameter σ does not have an influence on the loss and can therefore be safely ignored.

The gradient with respect to S is given by i

A i (A i S -X i ) = i (S -A i X i ) yielding the closed form updates: S ← 1 m i (A i X i ) (4.68)
From (4.67), minimizing L with respect to A i is equivalent to maximizing A i , X i S and therefore we have:

A i ← P( 1 n X i S ) (4. 69
)
where P is the projection on the Stiefel manifold:

P(M) = M(M M) -1 2
. The complexity of Deterministic SRM is in Õ(n iter mpvn) where n is the number of samples and n iter the number of iterations. We monitor the convergence by looking at the ∞ norm of the gradient. Note that we can monitor the gradient without any increase in complexity. Indeed, after the updates with respect to each mixing matrix have been carried out, only the gradient with respect to S remains: i (S -A i x i ). The algorithm is stopped when the gradient falls below a chosen tolerance.

Probabilistic SRM

In Probabilistic SRM , Σ i = σ 2 i I v and the shared components are assumed to be Gaussian s ∼ N(0, Σ s ). In [START_REF] Chen | A reduced-dimension fMRI shared response model[END_REF], Σ s is only assumed to be definite positive. However, as will be seen in the FastSRM chapter (chapter 5), enforcing a diagonal Σ s ensures identifiability (provided the diagonal values are different). So we assume that Σ s is diagonal.

The model is optimized via the expectation maximization algorithm.

Denoting V[s|x] = ( i 1 σ 2 i I + Σ -1 s ) -1 and E[s|x] = V[s|x] i 1 σ 2 i A i x i , we have p(x, s) = i exp(-x i -A i s 2 2σ 2 i ) (2πσ 2v i ) 1 2 exp(-1 2 s, Σ -1 s s ) (2π|Σ s |) 1 2 (4.70) = c 1 exp(- 1 2 i 1 σ 2 i x i 2 -2 i 1 σ 2 i A i x i , s (4.71)
+ i 1 σ 2 i s 2 + s, Σ -1 s s ) (4.72) = c 2 (x) exp(- 1 2 s -E[s, x], V[s|x] -1 (s -E[s|x]) ) (4.73)
where

c 1 = 1 (2πσ 2v i ) 1 2 1 (2π|Σ s |) 1 2 and c 2 (x) = c 1 exp(-1 2 ( i 1 σ 2 i x i 2 - E[s, x], V[s|x] -1 E[s|x] )) are independent of s. Therefore s|x ∼ N(E[s|x], V[s, x])
The negative expected completed log-likelihood is given by

L = i 1 2 v log(σ 2 i ) + 1 2σ 2 i E[ x i -A i s 2 ] (4.74)
updates are therefore given by:

σ 2 i ← 1 v (E[ x i -A i E[s|x] 2 ] + diag(V[s|x]) 2 ) (4.75) A i ← P(E[x i E[s|x] ]) (4.76) Σ s ← V[s|x] + E[E[s|x]E[s|x] ] (4.77)
It is useful to access the log-likelihood to check the implementation of the algorithm and monitor the convergence. From equation (4.73), the likelihood is given by:

p(x) = c 2 (x) s exp(- 1 2 s -E[s, x], V[s|x] -1 (s -E[s|x]) )ds (4.78) = c 2 (x)(2π|V[s|x]|) 1 2
(4.79) replacing c 2 (x) by its expression and taking the log, the expected negative log-likelihood is (up to constants) given by:

E[-log(p(x))] = i v 2 log(σ 2 i ) + 1 2 log(|Σ s |) - 1 2 log(|V[s|x]|) + i 1 2 1 σ 2 i E[ x i 2 ] - 1 2 E[ E[s, x], V[s|x] -1 E[s|x] ] (4.80)
The complexity of Probabilistic SRM is Õ(n iter mpvn), the same as in Deterministic SRM. We can monitor the convergence by looking at the log-likelihood decrease at each iteration and stop the algorithm when the magnitude of the decrease is below some tolerance. The storage requirements of Deterministic or Probabilistic SRM are in Õ(mvn) which simply means that the dataset needs to hold in memory.

conclusion

In this chapter, we have reviewed several methods to perform unsupervised analysis of neuroimaging data. We have introduced methods most suited to the analysis of the data of a single subject such as ICA and have explored some of the extensions to multiple subjects such as CanICA, ConcatICA, IVA or SRM. In the next chapter, we present our first contribution: an efficient implementation of SRM that we call FastSRM.

Part II FA S T S R M

: A N E F F I C I E N T I M P L E M E N TAT I O N O F T H E S H A R E D R E S P O N S E M O D E L FA S T S R M T H E O R Y
In the previous chapter, we have decribed multiple unsupervised methods to analyze multiview data. As described in section 4.2.5, the shared response model [START_REF] Chen | A reduced-dimension fMRI shared response model[END_REF] (SRM) is a multi-view latent factor model. It sees the data (x i ) m i=1 as:

x i = A i s + n i (5.1)
A i A i = I p (5.2)
with n i ∼ N(0, Σ i ) where

Σ i = σ 2 I v in deterministic SRM and Σ i = σ 2 i I v in probabilistic SRM.
In practice, we have access to n observations of x i stacked in a matrix X i ∈ R v×n . The corresponding observations of s called S ∈ R k×n are unknown.

When working with high dimensional data, SRM is particularly interesting as it provides a principled way to perform dimension reduction. Note that this contrasts with ICA-like methods that do not incorporate dimension reduction in their model.

However SRM has initially been designed to work within regions of interest using only few subjects. When using full brain data, computational costs become high. In addition, memory requirements are difficult to meet since the full dataset needs to hold in memory.

Fortunately, these high costs can be reduced. Intuitively, since the shared response lives in a reduced space, a compressed representation of the input is good enough to find a suitable estimate of the shared response. FastSRM implements this idea. It turns out that there exists an optimal compression of the input for which we can obtain the same solution as with the full data.

the fastsrm algorithm

SRM algorithms use different set of parameters θ to represent the data. In deterministic SRM θ = (A i ) m i=1 , S where (A i ) m i=1 are the mixing matrices and S are the n observations of the shared response S while in probabilistic SRM θ = (A i ) m i=1 , Σ s , (σ i ) m i=1 where (A i ) m i=1 are the mixing matrices, (σ i ) m i=1 the noise levels and Σ s the components variance.

In fMRI, the classical approach used to reduce the data is to apply an atlas. A deterministic atlas such as [START_REF] Bellec | Multi-level bootstrap analysis of stable clusters in resting-state fMRI[END_REF] is a parcellation of the brain into r regions. Reducing an image using a deterministic atlas corresponds to averaging the signal within each region of the atlas. A probabilistic atlases such as [START_REF] Dadi | Fine-grain atlases of functional modes for fMRI analysis[END_REF] describes each region as a set of weights across the full brain. Therefore, the image reduction can be done with a matrix product.

In FastSRM we consider a set of view specific atlases U i ∈ R v×r such that U i U i = I r where r is the number of regions in the atlas. Data are reduced using z i = U i x i and an SRM algorithm is applied on data z i yielding parameters θ . The figure 5.1 illustrates this process.

Note that the parameters obtained with FastSRM θ are different from the parameters obtained with the corresponding SRM algorithm θ (the unmixing matrices in θ do not even have the same shape as the unmixing matrices in θ). However, as we will see in the next section, there exists atlases such as the correspondence between θ and θ can be made explicit.

From a computational stand point, the dimension reduction provides a large reduction in memory usage. Indeed as the original data are seen only once, it is no longer necessary to keep the full dataset in memory (we can load data X i one after the other and similarly for the atlases U i ). Therefore the memory consumption is only in Õ(vn) (where v is the number of voxels and n is the number of samples) which is lower than SRM by a factor of m, the number of subjects. The number of subjects is typically between 10 and 1000. This yields a practical benefit: on fMRI datasets with many subjects, one no longer needs a large cluster to run the shared response model but only a modern laptop. Additionally, low memory consumption reduces the risk of thrashing [START_REF] Peter | Thrashing: Its causes and prevention[END_REF], a phenomenon that causes large increase in computation time when the memory used is close to the total available memory in the hardware.

After preprocessing, the reduced representation z i is used instead of the original data x i yielding a time complexity of Õ(T preprocessing + n iter mpnr). Let us highlight that an experiment is often run multiple times such as when cross validated results are needed. In these cases, the pre-processing is performed only once and the apparent complexity becomes Õ(n iter mpnr) which is faster than SRM by a factor of v r . The number of regions in large atlases is about r = 1000 and in full brain data, the number of voxels is about 300 000 so that v r is typically about 1000.

It remains to show how to draw a correspondence between FastSRM and SRM which is addressed in the following section.

optimal atlases

In principle, FastSRM can be used with any atlas. However, in general, working with reduced data induces a loss of information and therefore there is little hope to recover the parameters that would have been obtained from SRM from the parameters of FastSRM. However, in this section, we show that there exists an optimal atlas in the sense that SRM and FastSRM yield the same results.

Let us consider x i = U x i z i a PCA of x i with the maximum number of components. As the number of samples n is lower than the number of features, U x i ∈ R v×n and z i ∈ R n . We also have U x i U x i = I. Therefore U x i constitutes a possible choice of subject specific atlas.

As the next property shows, U x i is an optimal atlas for deterministic FastSRM.

Proposition 9 (Optimal atlas for deterministic FastSRM). Let (A i ) i , S be the solution obtained by deterministic SRM on data (X i ) i and (A i ) i , S the solution obtained by deterministic FastSRM on data (X i ) i using atlases (U X i ) i where

X i = U x i Z i is a PCA of X i . Then A i = U X i A i and S = S .
Proof. Updates of the mixing matrices A i in deterministic SRM equation (4.69) can be written:

A i ← P( 1 n X i S ) = U X i P( 1 n Z i S ) (5.3)
where P is the projection on the Stiefel manifold:

P(M) = M(M M) -1 2 . Therefore we can look for A i as A i = U X i Ãi . Ãi is orthogonal. Indeed A i A i = I p (5.4) =⇒ Ã i U X i U X i Ãi = I p (5.5) =⇒ Ã i Ãi = I p (5.6)
Then, we use the fact that

X i -A i S 2 = U X i Z i -U X i Ãi S 2 = Z i -Ãi S 2 (5.7) so that A i = Ãi .
Therefore, the solution of deterministic SRM on data (z i ) m i=1 and (x i ) m i=1 are linked by the change of parameters A i = U x i A i and S = S . This concludes the proof.

In the case of probabilistic SRM we can obtain very similar results. However the algorithm applied on reduced data need to be slightly modified. We call probSRM(ψ) the probabilistic SRM algorithm modified such that updates

σ 2 i ← 1 v (E[ x i -A i E[s|x] 2 ] + diag(V[s|x]) 2 ) (5.8)
are replaced by updates

σ 2 i ← 1 ψ (E[ x i -A i E[s|x] 2 ] + diag(V[s|x]) 2 ) (5.9)
We have the following result:

Proposition 10 (Optimal atlas for probabilistic FastSRM). Let (A i ) i , σ i , Σ s be the solution obtained by probabilistic SRM on data x i and

(A i ) i , σ i , Σ s the solution obtained by ProbSRM(v) on data z i = U x i x i .Then A i = U x i A i , σ i = σ i and Σ s = Σ s .
Proof. Updates of the mixing matrices A i in probabilistic SRM equation (4.76) can be written:

A i ← U x i P(E[z i E[s|x i ] T ]) (5.10)
so we can look for A i as A i = U x i Ãi and, as in the deterministic case, Ãi is orthogonal. Therefore equality (5.7) holds.

Then we consider the expected negative log-likelihood of probabilistic srm:

L = i 1 2 v log(σ 2 i ) + 1 2 log(|Σ s |) + E[ s i 1 2σ 2 i x i -A i s 2 + 1 2 s, Σ -1 s s ds] (5.11) = i 1 2 v log(σ 2 i ) + 1 2 log(|Σ s |) + E[ s i 1 2σ 2 i z i -Ãi s 2 + 1 2 s, Σ -1 s s ds] (5.12)
where we use equality (5.7). Optimizing the log-likelihood via expectation maximization yields the exact same updates as probabilistic srm on data z i except that updates

σ 2 i ← 1 t (E[ z i -Ãi E[s|z] 2 ] + diag(V[s|z]) 2 ) (5.13)
are replaced by updates

σ 2 i ← 1 v (E[ z i -Ãi E[s|z] 2 ] + diag(V[s|z]) 2 ) (5.14) so that Ãi = A i .
Therefore, the updates in both algorithms are linked by

A i = U x i A i , σ i = σ i and Σ s = Σ s .
This concludes the proof.

The properties 9 and 10 show that no information is lost by replacing x i ∈ R v by its reduced representation z i ∈ R n . A key property of the optimal atlas U x i is that it is valid whether or not the model for deterministic (respectively probabilistic) SRM holds.

A complexity analysis shows that finding the optimal atlas becomes the limiting step of the pipeline. Even with fast implementations, the subject specific PCA is costly. However FastSRM only works on z i so we do not need to know the value of U x i . In practice, we observe data X i ∈ R v×n and we want to get Z i ∈ R n×n such that X i = U x i Z i . This can be done by performing an SVD of X i X i yielding

X i X i = V i D i V i and setting Z i = D 1 2 i V i . Although computing the product X i X i has time complexity Õ(vt 2
) there is exactly vt 2 multiplications and additions so it costs a lot less than the PCA on full data. When estimates of the mixing matrices are needed, they can be obtained by applying equation (5.3) in the deterministic SRM case and equation (5.10) in the probabilistic SRM case which only costs Õ(mvp 2 ). In practice the cost of the matrix products X i X i is often still the limiting step of the pipeline (this depends on the number of iterations) but as we show in the next chapter, it is much more efficient than performing SRM on the full data. Note than if memory allows it, these matrix products can be computed in parallel.

Up to now, we have assumed that the covariance of components is diagonal in probabilistic SRM. In the next section we justify this assumption.

identifiablity of the shared response model

We first show why deterministic SRM and probabilistic SRM without any restrictions on the covariance of the components can only recover unmixing matrices up to an arbitrary rotation.

Let us consider data X i generated from deterministic SRM (meaning deterministic SRM holds exactly for these data) with mixing matrices A i and shared response S. Then deterministic SRM with parameters A i = A i R and S = R ,where R ∈ R k×k is an orthogonal matrix, also generates X i . This shows that deterministic SRM is not identifiable.

Similarly, in the probabilistic SRM case, if data x i are generated from the probabilistic SRM model with parameters A i , Σ s , σ 2 i (where Σ s can be any symmetric positive definite matrix) then the probabilistic SRM model with parameters A i R, R Σ s R, σ 2 i where R ∈ R k×k is an orthogonal matrix also generates x i . This shows that if no constraints are imposed on Σ s , then probabilistic SRM is not identifiable.

In Proposition 11, we show that if Σ s is assumed to be diagonal, then, under weak assumptions, probabilistic SRM is identifiable.

Proposition 11 (Identifiability of probabilistic SRM). Let x i be generated from a probabilistic SRM with parameters A i , Σ s , σ 2 i (where Σ s is diagonal positive with distinct values on the diagonal) and assume there exists another set of parameters A i , Σ s , σ i 2 (where Σ s is diagonal positive with distinct values on the diagonal) that also generate x i . Then if m 3, A i = A i P , Σ s = PΣ s P and σ 2 i = σ i 2 where P is a sign and permutation matrix (meaning P is the product of a diagonal matrix with values in {-1, 1} and a permutation matrix).

Proof. Let us consider E[x i x j ] for i = j. We have:

A i Σ s A j = A i Σ s A j (5.15)
up to re-ordering, equation (5.15) gives two singular value decompositions of the same matrix. Therefore, by unicity of the singular value decomposition we have: Σ s = P i Σ s P j and A i = A i P i and A j = A j P j where P i and P j are sign and permutation matrices. Since there are more than three subjects, there exists subject z such that Σ s = P i Σ s P z and therefore P i Σ s P z = P i Σ s P j so that P j = P z . So all sign and permutations are the same and we call P their common value. Then we consider

E[x i x j ] = A i ΣA i + σ 2 I v = A i Σ s A i + σ 2 I v so we get σ 2 = σ 2
Proposition (11) justifies that Σ s should be assumed to be diagonal in probabilistic SRM. In addition, working with a diagonal covariance matrix allows to speed up slightly the computations (although this does not change the time complexity of the algorithm).

conclusion

In this chapter, we have presented FastSRM, an efficient implementation of SRM that uses optimal atlases to speed up computations and reduce memory requirements without loss of performance. We have also discussed the identifiability of SRM. In the next chapter, we measure the performance of FastSRM in practice and compare it to available implementations.

In the previous chapter we introduced the FastSRM algorithm, a framework to efficiently compute shared responses from the full brain data of multiple subjects. In this chapter, we investigate the practical benefits of FastSRM on synthetic and real fMRI data. FastSRM reduces the fitting time and memory requirements considerably, making it possible to compute shared responses using a laptop in a reasonable amount of time even on datasets that do not hold in RAM. The efficiency of FastSRM allows us to apply SRM on a large number of subjects. As an example, we apply FastSRM on 647 subjects from the CamCAN dataset and study how the mixing matrices of the SRM model are predictive of age.

comparing fitting time and performance of fastsrm and srm on synthetic data

We generate synthetic data x i according to the model of Probabilistic SRM. The parameters σ i , A i and Σ s are generated randomly. We sample the value of the subject specific noise level from a normal distribution: σ i ∼ N(0, 0.1). The mixing matrices A i are obtained by sampling their coefficient from a standardized normal distribution. Lastly, the covariance of the shared response Σ s is diagonal and the diagonal values are obtained by sampling from a Dirichlet distribution with parameter (1 . . . 1). We set the number of voxels to v = 125 000, the number of subjects to m = 10 and the number of components to p = 50. We generate n = 1000 samples. We benchmark deterministic SRM, probabilistic SRM and their FastSRM counterparts in terms of fitting time and performance. Algorithms are designated by the atlas they use and therefore SRM algorithms described in section 4.2.5.2 and 4.2.5.1 are refered to as None because no atlas is used and FastSRM algorithms will have the label Optimal. Note that it would be possible to use FastSRM with sub-optimal atlases (there exists a wide variety of atlases available [START_REF] Bellec | Multi-level bootstrap analysis of stable clusters in resting-state fMRI[END_REF][START_REF] Mensch | Extracting Universal Representations of Cognition across Brain-Imaging Studies[END_REF][START_REF] Schaefer | Local-global parcellation of the human cerebral cortex from intrinsic functional connectivity MRI[END_REF]) but without any guarantees that the performance are the same as SRM.

We use a number of iterations between 1 and 100 and report the performance, fitting time and a measure of convergence. In FastSRM, we do not compute the unmixing matrices but only the shared response. We measure the performance of an algorithm by computing the error between the true component S ∈ R p×n and the predicted component Ŝ ∈ R p×n using the quantity: mse( Ŝ, S) = min A∈R p×p A Ŝ -S 2 F = 65 S Ŝ † Ŝ -S 2 F . This way of measuring errors is insensitive to the indeterminacies in DetSRM. We measure the fitting time in seconds. Lastly, we measure convergence by computing the gradient ∞ norm in case of DetSRM given by max(abs(G)) where G is the gradient and use the distance between consecutive values of the loss for ProbSRM.

Results are plotted in figure 6.1. We empirically see that the optimal approach is equivalent to using no atlas in terms of performance. This is predicted by the theory in the previous chapter where we demonstrate that these two algorithms yield exactly the same output from the same input. In general probabilistic methods give much better results than their deterministic counterpart. This shows the superiority of likelihood based methods. In terms of fitting time, FastSRM is about a thousand time faster than SRM after 100 iterations. When no atlas is used, the number of iterations has a very strong impact on performance while it has a small impact when the optimal atlas is used. Lastly, looking at the convergence curves, we see that even after 100 iterations, algorithms did not fully converge. This means that in practice a much larger number of iterations is needed which would yield an even higher difference in fitting time between methods using no atlas and methods using the optimal atlas.

experiments on fmri data

We evaluate the performance of FastSRM on three fMRI datasets of subjects exposed to naturalistic stimuli: Sherlock, Forrest and Cam-CAN (more details about these datasets are available in section 3.1.5).

Comparing fitting time, memory usage and performance on a timesegment matching experiment

The timesegment matching experiment is first introduced in [START_REF] Chen | A reduced-dimension fMRI shared response model[END_REF]. In a nutshell, the time-segment matching accuracy measures the similarity between two multivariate time-series by trying to localize a timesegment in one time-series by correlation with the other. In the context of movie watching, this measure has a lot of sense: if we split the movies in scenes and compute a representation per scene and per subject, it makes sense to assume that different subjects watching the movie would still have closer representation of the same scenes than of different scenes. This explains why timesegment matching is a standard evaluation of SRM-like methods also used in [START_REF] Swaroop | A computational model of shared fine-scale structure in the human connectome[END_REF], [START_REF] Samuel A Nastase | Leveraging shared connectivity to aggregate heterogeneous datasets into a common response space[END_REF] or [START_REF] Zhang | A searchlight factor model approach for locating shared information in multi-subject fMRI analysis[END_REF].

We now describe more precisely the experimental design. We split the runs into a train and test set. After fitting the model on the training set, we apply the unmixing matrices W i = A -1 i of each subject on the test set yielding individual components matrices. We estimate the shared responses by averaging the individual compo- As expected, when optimal atlases are used, the performance is the same as if no atlas is used but the fitting time is much lower. This is even more pronounced when the number of iterations is high (and looking at convergence curves, we see that more iterations could be performed to be closer to a stationary point).

nents of each subjects but one. We select a target time-segment (9 consecutive timeframes) in the shared responses and try to localize the corresponding time segment in the components of the left-out subject using a maximum-correlation classifier. The time-segment is said to be correctly classified if the correlation between the sample and target time-segment is higher than with any other time-segment (partially overlapping time windows are excluded). We use 5-Fold cross-validation across runs: the training set contains 80% of the runs and the test set 20%, and repeat the experiment using all possible choices for left-out subjects. The mean accuracy is reported in Figure 6.2 (bottom). When the optimal atlas is used, the accuracy is the same as when no atlas is used but the fitting time is reduced by a factor 10 to 100 and so is the memory usage.

We would like to highlight here that these experiments are not exactly the same as in [START_REF] Chen | A reduced-dimension fMRI shared response model[END_REF] as we use full brain data and they use regions of interest. However, the code used for this experiment is very similar to the tutorial in https://brainiak.org/tutorials/11-SRM/.

Predict age from spatial components

Because FastSRM is fast and memory efficient, it enables large-scale analysis of fMRI recordings of subjects exposed to the same naturalistic stimuli. We use all 647 subjects of the CamCAN dataset and demonstrate the usefulness of FastSRM by showing that the spatial components that it extracts from movie watching data are predictive of age. A key asset of FastSRM is that these spatial components can be visualized and therefore provide meaningful insights.

Note that while it is possible to do the same study using SRM, the table in figure 6.4 shows that the memory usage of SRM is 10x more important than FastSRM, almost reaching the memory limits of our cluster and is about 200x slower. FastSRM could be easily used with a much higher number of subjects while keeping the computation time and memory requirements reasonable with SRM would be impossible.

We now describe our age prediction pipeline. Functionally matched spatial components A i are obtained using FastSRM. They are divided into two groups (train and test data) where the train set contains 80% of the data and the test set 20%. Within the train set we split again our data into two groups: the first group is used to train one Ridge model per spatial components to predict age, the second group is used to train a Random Forest to predict age from Ridge predictions. This way of stacking models is similar to the pipeline used in [START_REF] Rahim | Joint prediction of multiple scores captures better individual traits from brain images[END_REF]. We use 5 fold cross validation to split the train set (so that the number of samples used to train the Random Forest is the number of elements in the train set). After the Random Forest is trained, we re-train one Ridge model per component using the full train set. On the test set each Ridge model makes a prediction and the predictions are aggregated When the optimal atlas is used, the accuracy is the same as when no atlas is used but the fitting time is reduced by a factor 10 to 100 and so is the memory usage. We first learn the spatial components from fMRI data using FastSRM. We learn one Ridge model per spatial components to predict age across subjects. Then, these models are aggregated using a Random Forest (like in [START_REF] Rahim | Joint prediction of multiple scores captures better individual traits from brain images[END_REF]).

using the Random Forest model. An illustration of the process is available in Figure 6. [START_REF] Ablin | Faster independent component analysis by preconditioning with Hessian approximations[END_REF].

In each Ridge model, the coefficient that determines the level of l2 penalization is set by generalized cross validation, an efficient form of leave-one-out cross validation.

The train and test sets are chosen randomly. In Figure 6.4, we report the average mean absolute error (MAE) on the test set averaged over the 5 splits. FastSRM predicts age with an accuracy much better than chance resulting in a mean absolute error (MAE) of 7.5 years. Note that this is far from being the most accurate method. Using a combination of different modalities, it is possible to obtain a MAE of approximately 4.5 years [START_REF] Alexander Engemann | Combining electrophysiology with MRI enhances learning of surrogate-biomarkers[END_REF]. However, as will be shown in the next paragraph, our method extracts a set of components specific to each individual making it more interpretable than other approaches.

In order to assess which spatial components are most predictive of age, we assess the feature importance via the Gini index defined in [START_REF] Breiman | Random forests[END_REF] or [START_REF] Louppe | Understanding variable importances in forests of randomized trees[END_REF] that measures the relative reduction in Gini impurity brought by each feature. Feature importance varies with different splits. We use the averaged feature importance over the 5 splits of our pipeline and plot the 3 most important spatial components according to this metric in Figure 6.4. These spatial components represent respectively 16%, 12% and 8% of the total feature importance and they highlight the the visual dorsal pathway, the precuneus and the visual ventral pathway respectively. The fact that averaged spatial components are interpretable and meaningful allows us to study the influence of age on brain networks involved in movie-watching. In Figure 6.5, we plot the most important spatial component averaged within groups of ages. We see that these spatial components evolve with age allowing us to visually identify which regions are meaningful. It turns out that aging is mostly reflected in brain activity as a fading of activity in the spatial correlates of movie watching, particularly in the dorsal visual cortex.

conclusion

As studies using naturalistic stimuli will tend to become more common and large within and across subjects, we need scalable models in terms of computation time and memory usage. This is what FastSRM provides. We show that while FastSRM is provably equivalent to its SRM counterpart in terms of performance, it has much lower fitting time and memory requirements.

FastSRM allows large scale analysis of fMRI data of subjects exposed to naturalistic stimuli. As one example of such analysis, we show that it can be used to predict age from movie-watching data. Interestingly, although FastSRM is an unsupervised model, it extracts meaningful networks and as such constitutes a practical way of studying subjects exposed to naturalistic stimuli.

We also show that individual information can be extracted from the fMRI activity when subjects are exposed to naturalistic stimuli. Our predictive model is reminiscent of that of [START_REF] Bijsterbosch | The relationship between spatial configuration and functional connectivity of brain regions[END_REF], that have shown that ICA components obtained from the decomposition of resting state data carry important information on individual characteristics. Our model inherits from all the weaknesses of SRM including the fact that mixing matrices are assumed to be orthogonal which is rather unrealistic. In later chapters, we see how this constraint can be released.

Algorithm Memory usage (in Go) Fitting time (in minutes)

FastSRM 18 9 ProbSRM 178 2232 L R L R L R L R L R L R

M U LT I V I E W I C A T H E O R Y

In chapter 5 and chapter 6, we have introduced a fast version of the shared response model (SRM). While SRM provides a useful dimension reduction framework, it assumes orthogonality of the mixing matrices, which is not biologically plausible.

In this chapter, we propose a novel group ICA method called Mul-tiView ICA. In contrast to most Group ICA methods, MultiViewICA is grounded in a probabilistic model of the problem and comes with statistical guarantees such as asymptotic efficiency.

MultiViewICA models each subject's dataset as a linear combination of a common components matrix with additive Gaussian noise. Importantly, we consider that the noise is on the components and not on the sensors. This greatly simplifies the likelihood of the model which can even be written in closed-form.

Despite its simplicity, MultiView ICA allows for an expressive representation of inter-subject variability through subject-specific functional topographies (mixing matrices) and variability in the individual response (with noise in the component domain). To the best of our knowledge, this is the first time that such a tractable likelihood is proposed for multi-subject ICA. The likelihood formulation shares similarities with the usual ICA likelihood, which allows us to develop a fast and robust alternate quasi-Newton method for its maximization.

We first introduce the MultiView ICA model, and show that it is identifiable. We then write its likelihood in closed form, and maximize it using an alternate quasi-Newton method. We also provide a sensitivity analysis for MultiView ICA, and show that the choice of the noise parameter in the algorithm has little influence on the output.

multiview ica for shared response modelling

Model, likelihood and approximation

Given m subjects, we model the data of subject i as a random vector x i ∈ R p such that:

x i = A i (s + n i ), i = 1, . . . , m (7.1)
where s ∈ R p are the shared independent components, n i ∈ R p is individual noise and A i ∈ R p×p are the individual mixing matrices, assumed to be full-rank. In practice we have access to n observations of x i assumed independent and identically distributed that we stack into a matrix X i ∈ R p,n . For simplicity, we assume that the components share the same density δ, so that the independence assumption is p(s) = p j=1 δ(s j ). Finally, we assume that the noise is Gaussian decorrelated with variance σ 2 , n i ∼ N(0, σ 2 I p ), and that the noise is independent across subjects and independent from the components.

The assumption of additive white noise on the components models individual deviations from the shared components s. It is equivalent to having noise on the sensors with covariance σ 2 A i (A i ) , i.e. a scaled version of the data covariance without noise.

Since the components are shared by the subjects, there are many more observed variables than components in the model: there are p components, while there are p × m observations. Therefore, model (7.1) can be seen as an instance of undercomplete ICA. The goal of multiview ICA is to recover the mixing matrices A i from observations of the x i . The following proposition extends the standard idenfitiability theory of ICA [START_REF] Comon | Independent component analysis, a new concept?[END_REF] to multiview ICA, and shows that recovering the components/mixing matrices is a well-posed problem up to scale and permutation.

Proposition 12 (Identifiability of MultiView ICA). Consider x i , i = 1 . . . m, generated from (7.1). Assume that x i = A i (s + n i ) for some invertible matrices A i ∈ R p×p , independent non-Gaussian components s ∈ R p and Gaussian noise n i . Then, there exists a scale and permutation matrix P ∈ R p×p such that for all i, A i = A i P.

The proof is available in appendix A.1.1. We propose a maximum-likelihood approach to estimate the mixing matrices. We denote by W i = (A i ) -1 the unmixing matrices, and view the likelihood as a function of W i rather than A i .

To derive the likelihood, we start by conditioning on s. Then, we make a variable transformation from x i to n i = W i x i -s, as opposed to the transformation to s as is usual in ICA. Using the probability transformation formula, we obtain

p x i |s (x i |s) = |W i |p n i (W i x i -s) (7.2)
where p n i is the density of n i . Note that the x i are conditionally independent given s, so we have:

p x|s (x|s) = m i=1 |W i |p n i (W i x i -s) (7.3)
and we next get the joint density as:

p x,s (x, s) = p s (s) m i=1 |W i |p n i (W i x i -s) (7.4)
Integrating out s and taking the log and expectation gives the expected negative likelihood:

L(W 1 , . . . , W m ) = - m i=1 log |W i | -E[log s exp - 1 2σ 2 m i=1 W i x i -s 2 p(s)ds ] (7.5)
up to additive constants.

The integral in 7.5 after factorization, is given by

s p j=1 exp - 1 2σ 2 m i=1 (W ij x i -s j ) 2 δ(s j )ds (7.6)
where W ij is the j-th line of W i . Denote y ij = W ij x i and sj = 1 m m i=1 y ij . Fix j, and drop it to simplify notation. Then we need to solve the integral

s exp - 1 2σ 2 m i=1 (y i -s) 2 δ(s)ds = s exp - 1 2σ 2 [m(s -s) 2 + m i=1 (y i -s) 2 ] δ(s)ds = exp - 1 2σ 2 m i=1 (y i -s) 2 z exp - m 2σ 2 z 2 δ(s -z)dz
where we have made the change of variable z = s -s. The remaining integral simply means that δ is smoothed by a Gaussian kernel, which can be computed exactly if δ is a Gaussian mixture. We therefore define f(s) =log z exp -m 2σ 2 z 2 δ(s -z)dz . The expected negative log-likelihood becomes

L(W 1 , . . . , W m ) = - m i=1 log |W i | + 1 2σ 2 m i=1 E[ W i x i -s 2 ] + E[f( s)] (7.7)
Multiview ICA is then performed by minimizing L, and the estimated shared components are S = i W i X i m . The negative log-likelihood L is quite simple, and importantly, can be computed easily given the parameters of the model and the data; it does not involve any intractable integral.

For one subject (m = 1), L(W 1 ) simplifies to the negative loglikelihood of ICA and we recover Infomax [START_REF] Anthony | An informationmaximization approach to blind separation and blind deconvolution[END_REF][START_REF] Cardoso | Infomax and maximum likelihood for blind source separation[END_REF], where the component log-pdf is replaced with the smoothed f.

Alternate quasi-Newton method for MultiView ICA

The parameters of the model are estimated by minimizing L. We propose a combination of quasi-Newton method and alternate minimization for this task. First, L is non-convex: it is only defined when the W i are invertible, which is a non-convex set. Therefore, we only look for local minima as usual in ICA. We propose an alternate minimization scheme, where L is alternatively diminished with respect to each W i . When all matrices W 1 , . . . , W m are fixed but one, W i , L can be rewritten, up to an additive constant

L i (W i ) = -log |W i | + 1 -1/m 2σ 2 E[ W i x i - m m -1 s-i 2 + f( 1 m W i x i + s-i )] (7.8)
with s-i = 1 m j =i W j x j . This function has the same structure as the usual maximum-likelihood ICA cost function: it is written

L i (W i ) = -log |W i | + E[g(W i x i )]
, where g(y) = p j=1 f(

y j m + s-i,j ) + 1-1/m 2σ 2 (y j -m m-1 s-i,j ) 2
where s-i,j is the j-th component of s-i . Fast quasi-Newton algorithms [START_REF] Ablin | Faster independent component analysis by preconditioning with Hessian approximations[END_REF][START_REF] Zibulevsky | Blind source separation with relative newton method[END_REF] have been proposed for minimizing such functions. We employ a similar technique as [START_REF] Zibulevsky | Blind source separation with relative newton method[END_REF], which we now describe.

Quasi-Newton methods are based on approximations of the Hessian of L i . The relative gradient (resp. Hessian) [START_REF] Amari | A new learning algorithm for blind signal separation[END_REF][START_REF] Cardoso | Equivariant adaptive source separation[END_REF] of L i is defined as the matrix G i ∈ R p×p (resp. tensor H i ∈ R p×p×p×p ) such that as the matrix E ∈ R p×p goes to 0, we have L i ((

I p + E)W i ) L i (W i ) + G i , W i + 1
2 E, H i E . Standard manipulations yield:

G i = E[ 1 m f ( s)(y i ) + 1 -1/m σ 2 (y i - m m -1 s-i )(y i ) ] -I p (7.9)
where

y i = W i x i . (H i ) abcd = δ ad δ bc + δ ac E[ 1 m 2 f (s a ) + 1 -1/m σ 2 y ib y id ] (7.10) for a, b, c, d = 1 . . . p Newton's direction is then -(H i ) -1 G i .
However, this Hessian is costly to compute (it has p 3 non-zero coefficients) and invert (it can be seen as a big p 2 × p 2 matrix). Furthermore, to enforce that Newton's direction is a descent direction, the Hessian matrix should be regularized in order to eliminate its negative eigenvalues [START_REF] Nocedal | Numerical optimization[END_REF], and H i is not guaranteed to be positive definite. These obstacles render the computation of Newton's direction impractical. Luckily, if we assume that the signals in y i are independent, several coefficients cancel, and the Hessian simplifies to the approximation

(H i ) abcd = δ ad δ bc + δ ac δ bd Γ i ab (7.11) with (Γ i ) ab = E[ 1 m 2 f (s a ) + 1 -1/m σ 2 (y ib ) 2 ]
This approximation is sparse: it only has p(2p -1) non-zero coefficients. In order to better understand the structure of the approximation, we can compute the matrix (H i M) for M ∈ R p×p . We find (H i M) ab = (Γ i ) ab M ab + M ba : H i M ab only depends on M ab and M ba , indicating a simple block diagonal structure of H i . The operator H i is therefore easily regularized and inverted:

(H i ) -1 M ab = Γ i ba M ab -M ba Γ i ab Γ i ba -1 .
Finally, since this approximation is obtained by assuming that the y i are independent, the direction -H -1 i G i is close to Newton's direction when the y i are close to independence, leading to fast convergence. Algorithm 2 alternates one step of the quasi-Newton method for each subject until convergence. A backtracking line-search is used to ensure that each iteration leads to a decrease of L i . The algorithm is stopped when the maximum norm of the gradients over one pass on each subject is below some tolerance level, indicating that the algorithm is close to a stationary point.

= -H -1 i G i Find a step size ρ such that L i ((I p + ρD)W i ) < L i (W i ) with line search Update s = s + ρ m DW i x i , W i = (I p + ρD)W i , tol= max(tol, G i ) end end return Estimated unmixing matrices W i , estimated shared components s 7.1.

Robustness to model misspecification

Algorithm 2 has two hyperparameters: σ and the function f. The latter is usual for an ICA algorithm, but the former is not. We study the impact of these parameters on the separation capacity of the algorithm, when these parameters do not correspond to those of the generative model (7.1).

Proposition 13. We consider the cost function L in eq. (7.7) with noise parameters σ and function f. Assume sub-linear growth on f : |f (x)| c|x| α + d for some c, d > 0 and 0 < α < 1. Assume that x i is generated following model (7.1), with noise parameter σ and density of the component d which need not be related to σ and f. Then, there exists a diagonal matrix

Λ such that (Λ(A 1 ) -1 , . . . , Λ(A m ) -1 ) is a stationary point of L, that is G 1 , . . . , G m = 0 at this point.
The proof is available in appendix A.1.2.

The sub-linear growth of f is a customary hypothesis in ICA which implies that d has heavier-tails than a Gaussian, and in appendix A.1.2 we provide other conditions for the result to hold. In this setting, the shared components estimated by the algorithm are S = Λ(S + 1 m m i=1 N i ), which is a scaled version of the best estimate of the shared components under the Gaussian noise hypothesis.

This proposition shows that, up to scale, the true unmixing matrices are a stationary point for Algorithm 2: if the algorithm starts at this point it will not move. The question of stability is also interesting: if the algorithm is initialized close to the true unmixing matrices, will it converge to the true unmixing matrix? In the appendix A.1.3, we provide an analysis similar to [START_REF] Cardoso | Blind signal separation: statistical principles[END_REF], and derive sufficient numerical conditions for the unmixing matrices to be local minima of L.

related work

In contrast to ConcatICA or CanICA (see section 4.2.2.1), MultiView ICA maximizes a likelihood, which brings statistical guarantees like consistency or asymptotic efficiency. Furthermore MultiViewICA finds individual and shared independent components in a single step. This differs from ConcatICA or GroupICA that require additional steps when individual components are needed such as back-projection [START_REF] Vince D Calhoun | A method for making group inferences from functional MRI data using independent component analysis[END_REF] or dual-regression [START_REF] Christian F Beckmann | Group comparison of resting-state fMRI data using multi-subject ICA and dual regression[END_REF].

The approach of [START_REF] Guo | A unified framework for group independent component analysis for multi-subject fMRI data[END_REF] (see section 4.2.2.2) optimizes the more general model x i = A i s + n i . The likelihood for this model involves an intractable high dimensional integral that is cumbersome to evaluate, and is then optimized with an EM algorithm using an inexact E-step. Having the simpler model x i = A i (s + n i ) leads to a closed-form likelihood, that can then be optimized by more efficient means. Note that in MultiView ICA, the noise can be interpreted as individual variability rather than sensor noise.

The SR-ICA approach of [START_REF] Zhang | A searchlight factor model approach for locating shared information in multi-subject fMRI analysis[END_REF] performs dimension reduction, merging of individual data and independent component estimation. It is therefore similar to our method. However, they propose to modify the FastICA algorithm [START_REF] Hyvarinen | Fast and robust fixed-point algorithms for independent component analysis[END_REF] in a rather heuristic way, without specifying an optimization problem, let alone maximizing a likelihood. In the experiments on fMRI data in appendix A.1.4, we obtain better performance with MultiView ICA than the reported performance of SR-ICA.

One strength of our model is that we only assume that the mixing matrices are invertible and still enjoy identifiability whereas some other approaches impose additional constraints. For instance tensorial methods [START_REF] Christian | Tensorial extensions of independent component analysis for multisubject fMRI analysis[END_REF] assume that the mixing matrices are the same up to diagonal scaling. Other methods impose a common mixing matrix [START_REF] Vince D Calhoun | fMRI activation in a visual-perception task: network of areas detected using the general linear model and independent components analysis[END_REF][START_REF] Cong | Validating rationale of group-level component analysis based on estimating number of sources in EEG through model order selection[END_REF][START_REF] Vera A Grin-Yatsenko | Independent component approach to the analysis of EEG recordings at early stages of depressive disorders[END_REF][START_REF] Monti | A unified probabilistic model for learning latent factors and their connectivities from highdimensional data[END_REF]. Like PCA, the Shared Response Model [START_REF] Chen | A reduced-dimension fMRI shared response model[END_REF] (SRM) assumes orthogonality of the mixing matrices. While the model defines a simple likelihood and provides an efficient way to reduce dimension, the orthogonal constraint may not be plausible.

Deep Learning methods, such as convolutional auto-encoders (CAE), can also be used to find the subject specific unmixing [START_REF] Po-Hsuan Chen | A convolutional autoencoder for multi-subject fMRI data aggregation[END_REF]. While these nonlinear extensions of the aforementioned methods are interesting, these models are hard to train and interpret. In the experiments on fMRI data in appendix A.1.4, we obtain better accuracy with MultiView ICA than that of CAE reported in [START_REF] Po-Hsuan Chen | A convolutional autoencoder for multi-subject fMRI data aggregation[END_REF].

A different path to multi-subject ICA is to extract independent components with individual ICA in each subject and align them. We propose a simple baseline approach to do so called PermICA. Inspired by the heuristic of the hyperalignment method [START_REF] James V Haxby | A common, highdimensional model of the representational space in human ventral temporal cortex[END_REF] we choose a reference subject and first match the components of all other subjects to the components of the reference subject. The process is then repeated multiple times, using the average of previously aligned components as a reference. Finally, group components are given by the average of all aligned components. We use the Hungarian algorithm to align pairs of mixing matrices [START_REF] Tichavsky | Optimal pairing of signal components separated by blind techniques[END_REF]. Alternative approaches involving clustering have also been developed [START_REF] Bigdely-Shamlo | Measure projection analysis: a probabilistic approach to EEG source comparison and multi-subject inference[END_REF][START_REF] Esposito | Independent component analysis of fMRI group studies by selforganizing clustering[END_REF].

Lastly, IVA based methods (see section 4.2.3) estimate view-specific components but shared components are not modeled explicitly.

conclusion

In this chapter, we have proposed a novel unsupervised algorithm that reveals latent components observed through different views. Using an independence assumption, we have demonstrated that the model is identifiable. In contrast to previous approaches, the proposed model leads to a closed-form likelihood, which we then optimize efficiently using a dedicated alternate quasi-Newton approach. Our approach enjoys the statistical guarantees of maximum-likelihood theory, while still being tractable.

In the next chapter, we evaluate the performance of MultiView ICA on synthetic data, on EEG and fMRI data and compare its performance to other GroupICA methods.

In chapter 7, we have introduced MultiViewICA and presented its theoretical properties. In this chapter, we empirically verify through extensive experiments on fMRI and MEG data that MultiView ICA improves component identification with respect to competing methods, suggesting that the expressiveness and robustness of this model make it a useful tool for multivariate neural signal analysis.

experimental setting

In the following, the noise parameter in MultiviewICA is always fixed to σ = 1. We use the function f(•) = log cosh(•), giving f (•) = tanh(•) (f is called the learning function or non-linearity [START_REF] Hyvärinen | Independent component analysis in the presence of Gaussian noise by maximizing joint likelihood[END_REF]). We use the Infomax cost function [START_REF] Anthony | An informationmaximization approach to blind separation and blind deconvolution[END_REF] with the same non-linearity to run standard ICA, with the Picard algorithm [START_REF] Ablin | Faster independent component analysis by preconditioning with Hessian approximations[END_REF] for fast and robust minimization of the cost function. Picard is applied with the default hyper-parameters. The code for MultiViewICA is available online at https://github.com/hugorichard/multiviewica.

We compare the following methods to obtain p components: PermICA is described in section 7.2. SRM is the FastSRM algorithm described in part II. ConcatICA is described in section 4.2.2.1. PCA+ConcatICA corresponds to ConcatICA applied on subject data that have been first individually reduced by PCA with p components. CanICA is described in section 4.2.2.1. We define the chance level as the performance of an algorithm that computes unmixing matrices and projections to lower dimensional space by sampling random numbers from a standard normal distribution.

The subject specific dimension reduction in MultiView ICA, Per-mICA, ConcatICA and CanICA is performed with SRM in fMRI experiments and subject-specific PCA in MEG experiments.

The shape of the input data (v, t, m) depend on the dataset used and the experiment performed. A description of the datasets used is available in section 3.1.5 for fMRI data and in section 3.2.5 for MEG data. The values for p (the number of components) is given in the descriptions of the experiments. Note that in the special case of the experiments on the sherlock fMRI dataset, the experimental pipeline is available on github https://github.com/hugorichard/multiviewica/ tree/master/real _ data _ experiments (including downloading and preprocessing of the data).

Since the cost function L is non-convex, having a good initialization can make a difference in the final result. We propose a two stage approach. We begin by applying PermICA on the datasets, which gives us a first set of unimixing matrices W 1 , . . . , W m . Note that we could also use ConcatICA for this task. Next, we perform a diagonal scaling of the mixing matrices, i.e. we find the diagonal matrices Λ 1 , . . . , Λ m such that L(Λ 1 W 1 , . . . , Λ m W m ) is minimized. To do so, we employ Algorithm 2 but only take into account the diagonal of the descent direction at each step: the update rule becomes W i ← (I p + ρDiag(D))W i . The initial unmixing matrices for Algorithm 2 are then taken as Λ 1 W 1 , . . . , Λ m W m . Empirically, we find that this two stage procedure allows the algorithm to start close to a satisfactory solution. A summary of our quantitative results on real data is available in appendix A.6.

synthetic experiment

We validate our method on synthetic data generated according to the model in equation (7.1). The components are generated i.i.d. from a Laplace density d(x) = 1 2 exp(-|x|). The mixing matrices A 1 , • • • , A m are generated with i.i.d. entries following a normal law. Each compared algorithm returns a sequence of estimated unmixing matrices W 1 , . . . , W m . The performance of an algorithm is measured by the reconstruction error between the estimated components and the true components.

We use m = 10 datasets, p = 15 components and n = 1000 samples. Each experiment is repeated with 100 random seeds. We vary the noise level in the data generation from 10 -2 to 10 and display the performance of algorithms in Figure 8.1.

Multiview ICA has uniformly better performance than the other algorithms, which illustrates the strength of maximum-likelihood based methods. In accordance with results of section 7.1, it is able to separate the components even with misspecified noise parameter and component density.

fmri experiments

We use the Forrest, Sherlock, Raiders and CLIPS datasets described in section 3.1.5.

Reconstructing the BOLD signal of missing subjects

We want to show that once unmixing matrices have been learned, they can be used to predict evoked responses across subjects. This can be used to perform missing data imputation when the sessions of some subjects are missing. In [START_REF] Zhang | Transfer learning on fMRI datasets[END_REF], the authors consider multiple fMRI datasets that share a subpart of their subjects and use the fact that some subjects are shared to transfer information across datasets. Our reconstruction experiment measures a related quantity: the methods ability to predict data of left-out subjects from other subjects.

In this experiment we apply a 6 mm spatial smoothing to all datasets. We split the data into three groups. First, we randomly choose 80% of all runs from all subjects to form the training set. Then, we randomly choose 80% of subjects and take the remaining 20% runs as testing set. The left-out runs of the remaining subjects form the validation set. The compared algorithms are run on the training set and evaluated using the testing and validation sets. After an algorithm is run on training data, it defines for each subject a forward operator that maps individual data to the space spanned by components and a backward operator that maps the component space to individual data. For instance in ICA the forward operator is the product of the dimensionality reduction projection and unmixing matrix. We estimate the shared responses on the testing set by applying the forward operators on the testing data and averaging. Finally, we reconstruct the individual data from subjects in the validation set by applying the backward operators to the shared responses. We measure the difference between the true signal and the reconstructed one using voxel-wise R 2 score. The R 2 score between two series x ∈ R n and y ∈ R n is defined as R 2 (x, y) = 1 -

1 n Var(y) n t=1 (x t -y t ) 2 , where Var(y) = 1 n n t=1 (y t -1 n n t =1 y t ) 2
is the empirical variance of y. The R 2 score is always smaller than 1, and equals 1 when x = y. The experiment is repeated 25 times with random splits to obtain error bars. Figure 8.2 provides a visual summary of our experimental procedure.

The R 2 score per voxel depends heavily on which voxels are considered. For example voxels in the visual cortex are better reconstructed in the sherlock dataset than in the forrest dataset. Performances are therefore given in terms of mean R 2 score inside a region of interest (ROI) in order to leave out regions where there is no useful information. In order to determine the ROIs, we focus on the R2 score per voxel between the BOLD signal reconstructed by ConcatICA and the actual bold signal. We run ConcatICA with 10, 20 and 50 components and select the voxels that obtained a positive R2 score for all sets of components. We discard voxels with an R2 score above 80% as they visually correspond to artefacts and apply a binary opening using a unit cube as the structuring element. The chosen regions are plotted in figure 8.3.

In Figure 8.4 (top) we report the mean R 2 score within regions of interest. MultiView ICA has similar or better performance than the other methods on all datasets. This demonstrates its ability to capture inter-subject variability, making it a candidate of choice to handle missing data or perform transfer learning.

For completeness, we plot in Figure 8.5, for ConcatICA, SRM and MultiViewICA, the R2 score per voxel using 50 components for datasets sherlock, forrest, raiders and clips. As could be anticipated from the task definition, forrest obtains high reconstruction accuracy in the auditory cortices, while clips shows good reconstruction in the visual cortex (occipital lobe mostly); the richer sherlock and raiders datasets yield good reconstructions in both domains, but also in other systems (language, motor). We can also see that data reconstructed by Multi-ViewICA are a better approximation of the original data than other methods. This is particularly obvious for the clips datasets where it is

Between subjects time-segment matching

We reproduce the time-segment matching experiment described in section 6.2.1 and display the results in Figure 8.4 (bottom). MultiView ICA yields a consistent and substantial improvement in accuracy compared to other methods on the four datasets. We see a marked improvement on the sherlock and forrest datasets. A possible explanation lies in the preprocessing pipeline. Sherlock data undergo a 6 mm spatial smoothing and Forrest data are acquired at a higher resolution (7T vs 3T for other data). This affects the signal to noise ratio.

In order to investigate the practical impact of the choice of hyperparameter σ, we compute the accuracy of the MultiView ICA algorithm with different choice of σ. Results are reported in Figure 8.6. MultiviewICA performs consistently well for a wide range of noise parameter values, and only breaks at very high values. It supports the theoretical claim issued in Proposition 13 that the noise parameter is of little importance.

In appendix A.4, we plot the average forward operator across subjects of MultiView ICA and ConcatICA with 5 components on the forrest, sherlock, raiders and clips datasets.

Between-runs time-segment matching

We measure the ability of each algorithm to extract meaningful shared components that correlate more when they correspond to the same stimulus than when they correspond to distinct stimuli. We use the raiders-full dataset, which allows this kind of analysis because subjects watch some selected scenes from the movie twice, during the first two runs (1 and 2) and the last two (11 and 12). First, the forward operators are learned by fitting each algorithm with 20 components on the data of all 11 subjects using all 12 runs. We then select a subset of 8 subjects and the shared components are computed by applying the forward operators and averaging. We select a large target time-segment (50 timeframes) taken at random from run 1 and 2, and we try to localize the corresponding sample time-segment from the 10 last runs using a single component of the shared components.

The time-segment is said to be correctly classified if the correlation between the target and corresponding sample time-segment is higher than with any other time-segment (partially overlapping windows are excluded). In contrast to the between subject time-segment matching experiment, we obtain one accuracy score per component. We repeat the experiment 10 times with different subsets of subjects randomly chosen and report the mean accuracy of the three best performing components in Figure 8.7. Error bars correspond to a 95 % confidence interval. MultiView ICA achieves the highest accuracy.

We then focus on the 3 best performing components of MultiView ICA. For each component, we plot in Figure 8.8 (left) the shared components during two sets of runs where subjects were exposed to the same scenes of the movie. We then study the localisation of these components. We average the forward operators across subjects and plot the columns corresponding to the components of interest in Figure 8.8 (right). As each column is seen as a set of weights over all voxels, it represents a spatial map. We extract 20 components and report the mean accuracy of the three best performing components

The component 1 of the shared responses follows almost the same pattern in the two set of runs corresponding to the same scenes of the movie. The spatial map corresponding to component 1 highlights the language network. In component 2, the temporal patterns during the viewing of identical scenes are also very similar. The corresponding spatial map highlights the visual network especially the visual dorsal pathway. In component 3, there exists a similarity however less striking than with the two previous components. The corresponding spatial map highlights a contrast between the spatial attention network and the auditory network.

phantom meg data

We demonstrate the usefulness of our approach on MEG data using the Sinusoidal Phantom MEG dataset where m = 8 dipoles at different locations produce a known sinusoidal oscillation (more details in section 3.2.5). 100 epochs are available. For each dipole, we chose N e = 2, . . . , 16 epochs at random among our set of 100 epochs and concatenate them in the temporal dimension. We then apply algorithms on these data to extract p = 20 shared components. As we know the true component (the timecourse of the dipole), we can compute the reconstruction error of each component as the squared norm of the difference between the estimated component and the true component, after normalization to unit variance and fixing the sign. We only retain the component yielding minimal error. We also estimate for each forward operator the localization of the component by performing dipole fitting using its column corresponding to the component of minimal error. We then compute the distance of the estimated dipole to the true dipole. These metrics are reported in figure 8.9 when the number of epochs considered N e varies. MultiView ICA requires fewer epochs to correctly reconstruct and localize the true component.

experiment on camcan dataset

Finally, we apply MultiView ICA on the CamCAN dataset [START_REF] Taylor | The Cambridge Centre for Ageing and Neuroscience (Cam-CAN) data repository: structural and functional MRI, MEG, and cognitive data from a cross-sectional adult lifespan sample[END_REF]. A detailed description of the CamCAN dataset is available in section 3.2.5.

We use the magnetometer data from the MEG of 200 subjects chosen randomly. Each subject is repeatedly presented an audio-visual stimulus. The MEG signal corresponding to these trials are then timeaveraged to isolate the evoked response, yielding individual data. MultiView ICA is then applied to extract 20 shared components. 9 components were found to correspond to noise by visual inspection, and the 11 remaining are displayed in Figure 8.9. We observe that MultiView ICA recovers a very clean sequence of evoked potentials with sharp peaks for early components and slower responses for late components. In order to visualize their localization, we perform component localization for each subject by solving the inverse problem using sLORETA [START_REF] Domingo | Standardized lowresolution brain electromagnetic tomography (sLORETA): technical details[END_REF], providing a component estimate for each component. Then, we register each component estimate to a common reference brain. Finally, the component estimates are averaged, and thresholded maps are displayed in Figure 8.9. Individual maps corresponding to each component are displayed in Appendix A.3. The figure highlights both early auditory and visual cortices, also sug- gesting a propagation of the activity towards the ventral regions and higher level visual areas.

conclusion

In this chapter, we have demonstrated the usefulness of MultiView ICA for neuroimaging group studies both on fMRI and MEG data, where it outperforms other methods. A limiting aspect of MultiView ICA is the assumption that the noise variance is the same across subjects. This is limiting because it does not properly model between-subjects variability. In the next chapter, we propose an extension of MultiView ICA with a more realistic noise model.

Part IV S H A R E D I C A S H A R E D I C A T H E O R Y

In chapter 7 and chapter 8, we have introduced MultiView ICA, a well principled method to perform shared response modeling. While MultiView ICA yields good practical results, it does not model subject specific deviations from the shared response. Yet, the magnitude of the response may differ across subjects [START_REF] Penny | Random effects analysis[END_REF], as does any noise due to heart beats, respiratory artefacts or head movements [START_REF] Thomas | Noise contributions to the fMRI signal: An overview[END_REF]. This drawback is shared by most GroupICA methods that often rely on single subject ICA to recover the shared response. In addition, such methods are typically unable to separate Gaussian components.

In contrast, the framework of Independent vector analysis (IVA) [START_REF] Anderson | Joint blind source separation with multivariate Gaussian model: Algorithms and performance analysis[END_REF][START_REF] Lee | Independent vector analysis (IVA): multivariate approach for fMRI group study[END_REF] allows subject specific variability in the unmixed data. Current implementations such as IVA-L [START_REF] Lee | Independent vector analysis (IVA): multivariate approach for fMRI group study[END_REF], IVA-G [START_REF] Anderson | Joint blind source separation with multivariate Gaussian model: Algorithms and performance analysis[END_REF], IVA-L-SOS [START_REF] Bhinge | Extraction of time-varying spatiotemporal networks using parameter-tuned constrained IVA[END_REF], IVA-GGD [START_REF] Anderson | Independent vector analysis: Identification conditions and performance bounds[END_REF] or IVA with Kotz distribution [START_REF] Anderson | Independent vector analysis, the Kotz distribution, and performance bounds[END_REF] estimate view-specific components but shared components are not modeled explicitly. Studying the components post IVA can give great insights about which components are shared and by which subjects (see [START_REF] Long | Independent vector analysis for common subspace analysis: Application to multi-subject fMRI data yields meaningful subgroups of schizophrenia[END_REF]) and individual responses that are closed enough could be merged to produce a shared response. However, the flexibility in modeling different components from each subject comes at the cost of not having a well principled way of aggregating the components into a common shared response.

In this chapter, we introduce Shared ICA (ShICA), where each dataset is modeled as a linear transform of shared independent components contaminated by additive Gaussian noise. ShICA allows for principled extraction of the shared components (or responses) in addition to view-specific components. Since it incorporates a statistically sound noise model, it enables optimal inference of the shared responses (minimizing the mean squared error).

We first analyse the theoretical properties of the ShICA model, before providing powerful inference algorithms. First, we exhibit necessary and sufficient conditions for ShICA to be identifiable (previous work only shows local identifiability [START_REF] Anderson | Independent vector analysis: Identification conditions and performance bounds[END_REF]), in the presence of Gaussian or non-Gaussian components. We then introduce an algorithm called ShICA-J that uses Multiset CCA to fit the model when all the components are assumed to be Gaussian. We exhibit necessary and sufficient conditions for Multiset CCA to be able to fit the model (previous work only gives sufficient conditions [START_REF] Li | Joint blind source separation by multiset canonical correlation analysis[END_REF]) and provide examples on which ShICA-J can recover the mixing matrices, while Multiset CCA cannot. We next point out a practical problem, namely that even a small sampling noise can lead to large rotations of unmixing matrices when Multiset CCA is used. To address this issue and recover the correct unmixing matrices, we propose to apply joint diagonalization to the result of Multiset CCA. We further introduce ShICA-ML, a maximum likelihood estimator of ShICA that models non-Gaussian components using a Gaussian mixture model. While ShICA-ML yields more accurate components, ShICA-J is significantly faster and offers a great initialization to ShICA-ML. 9.1 shared ica (shica): an identifiable multi-view model

We assume a similar generative model as in MultiViewICA:

x i = A i (s + n i ) (9.1)
Like in MultiView ICA we assume that the shared components are statistically independent p(s) = p j=1 p(s j ), and that the individual noises are Gaussian and independent from the shared components. We assume n i ∼ N(0, Σ i ), where the matrices Σ i are assumed diagonal and positive. This contrasts with MultiView ICA: the individual noise variances are learned and are not assumed to be the same across components or subjects. We further assume that there are at least 3 views: m 3.

In contrast to almost all existing works, we assume that some components (possibly all of them) may be Gaussian, and denote G the set of Gaussian components: s j ∼ N(0, 1) for j ∈ G. The other components are non-Gaussian: for j / ∈ G, s j is non-Gaussian.

identifiability The parameters of the model are Θ = (A 1 , . . . , A m , Σ 1 , . . . , Σ m ). We are interested in the identifiability of this model: given observations x 1 , . . . , x m generated with parameters Θ, are there some other Θ that can generate the same observations? Let us consider the following assumption that requires that the individual noises for Gaussian components are sufficiently diverse:

Assumption 14 (Noise diversity in Gaussian components). For all j, j ∈ G, j = j , the sequences (Σ ij ) i=1...m and (Σ ij ) i=1...m are different where Σ ij is the j, j entry of Σ i

It is readily seen that there is one trivial set of indeterminacies in the problem: if P ∈ R p×p is a sign and permutation matrix the parameters (A 1 P, . . . , A m P, P Σ 1 P, . . . , P Σ m P) also generate x 1 , . . . , x m . The following theorem shows that under the above assumption, these are the only indeterminacies of the problem.

Theorem 15 (Identifiability). We suppose Assumption 14. We let Θ = (A 1 , . . . , A m , Σ 1 , . . . , Σ m ) another set of parameters, and assume that they also generate x 1 , . . . , x m . Then, there exists a sign and permutation matrix P such that for all i, A i = A i P, and Σ i = P Σ i P.

Proof. By hypothesis, the covariances verify

C ij = E[x i x j ] = A i (I p + δ ij Σ i )A j = A i (I p + δ ij Σ i )A j
for all i, j. We let P i = A -1 i A i . The previous relationship for j = i gives P i P j = I p . Because there are more than 3 views, there is another integer k / ∈ {i, j}, and we have P i P k = P j P k = I p . This shows that P i = P j : all these matrices are equal, and we call P their common value. The previous equation also gives PP = I p , so P is orthogonal. We have that s + n i and s + n i have independent components and s + n i = P(s + n i ). Lemma 21 in appendix B.1 (a direct consequence of classical ICA results [START_REF] Comon | Independent component analysis, a new concept?[END_REF], Theorem 10) gives P = Π -1 ΩΠ where Π and Π are sign and permutation matrices such that the first g components of Π(s + n i ) and Π (s + n i ) are Gaussian, and Ω is a block diagonal matrix given by

Ω = Ω g 0 0 I p-g
where Ω g is orthogonal. We call A (g) the first g × g block of a matrix A so that Ω (g) = Ω g . Then, considering only the Gaussian components, we can write for i = j: (ΠΣ i ) (g) = Ω g (Π Σ i ) (g) Ω g for all i. This, combined with Assumption 14, implies that Ω g is a sign and permutation matrix (see Lemma 22 in appendix B.1) and therefore P is a sign and permutation matrix. Then it follows that I + Σ i = P(I + Σ i )P and therefore Σ i = PΣ i P so Σ i = P Σ i P.

Identifiability in the Gaussian case is a consequence of the identifiability results in [START_REF] Vía | Joint blind source separation from second-order statistics: Necessary and sufficient identifiability conditions[END_REF] and in the general case, local identifiability results can be derived from the work of [START_REF] Anderson | Independent vector analysis: Identification conditions and performance bounds[END_REF]. However local identifiability only shows that for a given set of parameters there exists a neighborhood in which no other set of parameters can generate the same observations [START_REF] Thomas | Identification in parametric models[END_REF]. In contrast, the proof of Theorem 15 shows global identifiability.

Theorem 15 shows that the task of recovering the parameters from the observations is a well-posed problem, under the sufficient condition of Assumption 14. We also note that Assumption 14 is necessary for identifiability. For instance, if j and j are two Gaussian components such that Σ ij = Σ ij for all i, then a global rotation of the components j, j yields the same covariance matrices. The current work assumes m 3. In appendix B.2 we give an identifiability result for m = 2.

estimation of components with noise diversity via joint-diagonalization

We now consider the computational problem of efficient parameter inference. This section considers components with noise diversity, while the next section deals with non-Gaussian components.

Fitting ShICA via Multiset CCA

If we assume that the components are all Gaussian, the covariance of the observations given by C ij = E[x i x j ] = A i (I p + δ ij Σ i )A j are sufficient statistics and methods using only second order information, like Multiset CCA, are candidates to estimate the parameters of the model. Consider the matrix C ∈ R pm×pm containing m × m blocks of size p × p such that the block i, j is given by C ij . Consider the matrix D identical to C excepts that the non-diagonal blocks are filled with zeros. Multiset CCA (using the SUMCORR cost function under some constraint as described in section 4.2.1) consists in the following generalized eigenvalue problem:

Cu = λDu, λ > 0, u ∈ R pm . ( 9.2)
Consider the matrix U = [u 1 , . . . , u p ] ∈ R mp×p formed by concatenating the p leading eigenvectors of the previous problem ranked in decreasing eigenvalue order. Then, consider U to be formed of m blocks of size p × p stacked vertically and define (W i ) to be the i-th block. These m matrices are the output of Multiset CCA. We also denote λ 1 • • • λ p the p leading eigenvalues of the problem.

The next theorem shows that we only need λ 1 . . . λ p to be distinct for Multiset CCA to solve ShICA: Assumption 16 (Unique eigenvectors). λ 1 . . . λ p are distinct. Theorem 17. We suppose Assumption 16 (only). Then, there exists a permutation matrix P and scale matrices Γ i such that W i = PΓ i A -1 i for all i.

Proof. Let us denote W ∈ R mp×mp the block diagonal matrix with block i given by (A i ) -1 . We have Cu = λDu ⇐⇒ WCW z = λWDW z where u = W z. We call z a reduced eigenvector. Each block in WCW and in WDW is diagonal so any reduced eigen-

vector z =     z 1 . . . z m     is such that the matrix Z = [z 1 . . . z m ] has ex-
actly one non-zero line. Following Lemma 23 in appendix B.1, the first p leading reduced eigenvectors z 1 , . . . , z p all have different first non-zero coordinates. Therefore the concatenation of the first p leading reduced eigenvectors is given by [z 1 , . . .

z p ] =     Γ 1 . . . Γ m     P
where P ∈ R p×p is a permutation matrix and Γ i ∈ R p×p is a diagonal matrix. Therefore, the first p eigenvectors are given by

[u 1 . . . u p ] =     W 1 . . . W m     =     (A -1 1 ) Γ 1 P . . . (A -1 m ) Γ m P     and so W i = PΓ i A -1 i
This theorem means that solving the generalized eigenvalue problem (9.2) allows to recover the mixing matrices up to a scaling and permutation: this form of generalized CCA recovers the parameters of the statistical model. Note that Assumption 16 is also a necessary condition. Indeed, if two eigenvalues are identical, the eigenvalue problem is not uniquely determined.

Note that the link between Multiset CCA and probabilistic models has been studied in other contexts [START_REF] Li | Joint blind source separation by multiset canonical correlation analysis[END_REF] or [START_REF] Adali | Diversity in independent component and vector analyses: Identifiability, algorithms, and applications in medical imaging[END_REF] although the authors do not use the same formulation of MultisetCCA as ours. We highlight here, that the different versions of Multiset CCA are not equivalent and may very well produce different results.

We have two different Assumptions, 14 and 16, the first of which guarantees theoretical identifiability as per Theorem 15 and the second guarantees consistent estimation by Multiset CCA as per Theorem 17. Next we will discuss their connections, and show some limitations of the Multiset CCA approach. To begin with, we have the following result about the eigenvalues of the problem (9.2) and the Σ ij . Proposition 18. For j p, let λ j the largest solution of

m i=1 1 λ j (1+Σ ij )-Σ ij = 1.
Then, λ 1 , . . . , λ p are the p largest eigenvalues of problem (9.2).

It is easy to see that we then have λ 1 , . . . , λ p greater than 1, while the remaining eigenvalues are lower than 1. From this proposition, two things appear clearly. First, Assumption 16 implies Assumption 14. Indeed, if the λ j 's are distinct, then the sequences (Σ ij ) i must also be different from the previous proposition. This is expected as from Theorem 17, Assumption 16 implies identifiability, which in turn implies Assumption 14.

Proposition 18 also allows us to derive cases where Assumption 14 holds but not Assumption 16. The following proposition shows that we can chose parameters of the model so that the model is identifiable but it cannot be solved using Multiset CCA: Proposition 19. Assume that for two integers j, j , the sequence (Σ ij ) i is a permutation of (Σ ij ) i , i.e. that there exists a permutation of {1, . . . , p}, π, such that for all i, Σ ij = Σ π(i)j . Then, λ j = λ j .

In this setting, Assumption 14 holds so ShICA is identifiable, while Assumption 16 does not hold, so Multiset CCA cannot recover the unmixing matrices.

Sampling noise and improved estimation by joint diagonalization

The consistency theory for Multiset CCA developed above is conducted under the assumption that the covariances C ij are the true covariances of the model, and not approximations obtained from observed samples. In practice, however, a serious limitation of Multiset CCA is that even a slight error of estimation on the covariances, due to "sampling noise", can yield a large error in the estimation of the unmixing matrices, as will be shown next. We begin with an empirical illustration. We take m = 3, p = 2, and Σ i such that λ 1 = 2 + ε and λ 2 = 2 for ε > 0.

In this way, we can control the eigengap of the problem, ε. We take W i the outputs of Multiset CCA applied to the true covariances C ij . Then, we generate a perturbation ∆ = δ • S, where S is a random positive symmetric pm × pm matrix of norm 1, and δ > 0 controls the scale of the perturbation. We take ∆ ij the p × p block of ∆ in position (i, j), and Wi the output of Multiset CCA applied to the covariances C ij + ∆ ij . We finally compute the sum of the Amari distance between the W i and Wi : the Amari distance measures how close the two matrices are, up to scale and permutation [START_REF] Amari | A new learning algorithm for blind signal separation[END_REF]. Fig 9 .1 displays the median Amari distance over 100 random repetitions, as the perturbation scale δ increases. The different curves correspond to different values of the eigen-gap ε. We see clearly that the robustness of Multiset CCA critically depends on the eigen-gap, and when it is small, even a small perturbation of the input (due, for instance, to sampling noise) can lead to large estimation errors.

This problem is very general and well studied [START_REF] Stewart | Error and perturbation bounds for subspaces associated with certain eigenvalue problems[END_REF]: the mapping from matrices to (generalized) eigenvectors is highly non-smooth. However, the gist of our method is that the span of the leading p eigenvectors is smooth, as long as there is a large enough gap between λ p and λ p+1 . For our specific problem we have the following bounds, derived from Prop. 18.

Proposition 20. We let σ max = max ij Σ ij and σ min = min ij Σ ij . Then, λ p 1 + m-1 1+σ max , while λ p+1 1 - 1 1+σ min .
As a consequence, we have λ p -λ p+1 m-1 1+σ max + 1 1+σ min : the gap between these eigenvalues increases with m, and decreases with the noise power.

Algorithm 3: ShICA-J Input: Covariances Cij = E[x i x j ] ( Wi ) i ← MultisetCCA(( Cij ) ij ) Q ← JointDiag(( Wi Cii W i ) i ) Γ ij ← Q Wi Cij W j Q (Φ i ) i ← Scaling((Γ ij ) ij ) return Unmixing matrices (Φ i Q Wi ) i
In this setting, when the magnitude of the perturbation ∆ is smaller than λ p -λ p+1 , [START_REF] Stewart | Error and perturbation bounds for subspaces associated with certain eigenvalue problems[END_REF] indicates that Span([W 1 , . . . , W m ] ) Span([ W1 , . . . , Wm ] ), where [W 1 , . . . , W m ] ∈ R pm×p is the vertical concatenation of the W i 's. In turn, this shows that there exists a matrix

Q ∈ R p×p such that W i Q Wi for all i. (9.3)
We propose to use joint-diagonalization to recover the matrix Q. Given the Wi 's, we consider the set of symmetric matrices Ki = Wi Cii W i , where Cii is the contaminated covariance of x i . Following Eq. ( 9.3), we have Q Ki Q = W i Cii W i , and using Theorem 17, we have

Q Ki Q = PΓ i A -1 i Cii A - i Γ i P . Since Cii is close to C ii = A i (I p + Σ i )A i , the matrix PΓ i A -1 i Cii A - i Γ i P is almost diag- onal.
In other words, the matrix Q is an approximate diagonalizer of the Ki 's, and we approximate Q by joint-diagonalization of the Ki 's. In Fig 9 .1, we see that this procedure mitigates the problems of multiset-CCA, and gets uniformly better performance regardless of the eigen-gap. In practice, we use a fast joint-diagonalization algorithm [START_REF] Ablin | Beyond Pham's algorithm for joint diagonalization[END_REF] to minimize a joint-diagonalization criterion for positive symmetric matrices [START_REF] Dinh | Joint approximate diagonalization of positive definite Hermitian matrices[END_REF]. The estimated unmixing matrices

U i = Q Wi (9.4)
correspond to the true unmixing matrices only up to some scaling: the information that the components are of unit variance is lost.

Scale estimation

We form the matrices Γ ij = U i Cij U j . In order to estimate the scalings, we solve

L Φ = min (Φ i ) i =j Φ i diag(Γ ij )Φ j -I p 2 F (9.5)
where the Φ i are diagonal matrices. The gradient is given by

∂L Φ ∂Φ i = 2 j =i (Φ i diag(Y ij )Φ j -I p )Φ j (9.6)
Therefore we get

∂L Φ ∂Φ i = 0 (9.7) ⇐⇒ 2 j =i (Φ i diag(Y ij )Φ j -I p )Φ j = 0 (9.8) ⇐⇒ Φ i j =i diag(Y ij )Φ 2 j - j =i
Φ j = 0 (9.9)

⇐⇒ Φ i = j =i Φ j j =i diag(Y ij )Φ 2 j (9.10)
We then iterate formula (9.10) over i until convergence. The final estimates of the unmixing matrices are given by

(Φ i U i ) m i=1 (9.11)
The full procedure, called ShICA-J, is summarized in Algorithm 3.

Estimation of noise covariance and inference of shared components

In practice, it is important to estimate noise co-variances Σ i in order to take advantage of the fact that some views are noisier than others. As it is well known in classical factor analysis, modelling noise variances allows the model to virtually discard variables, or subjects, that are particularly noisy. Using the ShICA model with Gaussian components, we derive noise covariances estimate directly from maximum likelihood. We use an expectation-maximization (EM) algorithm, which is especially fast because noise updates are in closed-form. Following derivations given in appendix B.3.1, the sufficient statistics in the E-step are given by

E[s|x] = m i=1 Σ -1 i + I -1 m i=1 Σ -1 i y i (9.12) V[s|x] = ( m i=1 Σ -1 i + I) -1 (9.13)
Incorporating the M-step we get the following updates that only depend on the covariance matrices:

Σ i ← diag( Ĉii -2V[s|x] m j=1 Σ -1 j Ĉji + V[s|x] m j=1 m l=1 Σ -1 j Ĉjl Σ -1 l V[s|x] + V[s|x]) (9.14)
9.3 shica-ml: maximum likelihood for non-gaussian components ShICA-J only uses second order statistics. However, the ShICA model (9.1) allows for non-Gaussian components. We now propose an algorithm for fitting the ShICA model that assumes non-Gaussian components so that it can separate Gaussian and non-Gaussian components. We estimate the parameters by maximum likelihood. Since most non-Gaussian components in real data are super-Gaussian, we assume that the non-Gaussian components s have the super-Gaussian density

p(s j ) = 1 2 N(s j ; 0, 1 2 ) + N(s j ; 0, 3 2 ) (9.15)
We propose to maximize the expected log-likelihood using a generalized EM [START_REF] Arthur P Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF][START_REF] Radford | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF]. Derivations are available in Appendix B.4. Like in the previous section, the E-step is in closed-form yielding the following sufficient statistics:

E[s j |x] = α∈{ 1 2 , 3 2 } θ α α ȳj α+ Σj α∈{0.5,1.5} θ α (9.16) V[s j |x] = α∈{ 1 2 , 3 2 } θ α Σj α α+ Σj α∈{0.5,1.5} θ α (9.17)
where θ α = N( ȳj ; 0, Σj + α), ȳj

= i Σ -1 ij y ij i Σ -1 ij
and Σj = ( i Σ -1 ij ) -1 with y i = W i x i . Noise updates are in closed-form and given by:

Σ i ← diag(E[(y i -E[s|x])(y i -E[s|x]) ]) + V[s|x]) (9.18)
However, no closed-form is available for the updates of unmixing matrices. We therefore perform quasi-Newton updates given by

W i ← (I -ρ( H W i ) -1 G W i )W i (9.19)
where ρ ∈ R is chosen by backtracking line-search

H W i a,b,c,d = δ ad δ bc + δ ac δ bd E[(y ib ) 2 ] Σ ia (9.20)
is an approximation of the Hessian of the expected negative complete log-likelihood and

G W i = -I + (Σ i ) -1 E[(y i -E[s|x])(y i ) ] (9.21)
is the gradient.

We alternate between computing the statistics E[s|x], V[s|x] (E-step) and updates of parameters Σ i and W i for i = 1 . . . m (M-step). Let us highlight that our EM algorithm and in particular the E-step resembles the one used in [START_REF] Moulines | Maximum likelihood for blind separation and deconvolution of noisy signals using mixture models[END_REF]. However because they assume noise on the sensors and not on the components, their formula for E[s|x] involves a sum with 2 p terms whereas we have only 2 terms. The resulting method is called ShICA-ML. minimum mean squared error estimates in shica In ShICA-J as well as in ShICA-ML, we have a closed-form for the expected components given the data E[s|x], shown in equation (9.12) and (9.16) respectively. This provides minimum mean squared error estimates of the shared components, and is an important benefit of explicitly modelling shared components in a probabilistic framework.

related work

ShICA combines theory and methods coming from different branches of "component analysis". It can be viewed as a GroupICA method, as an extension of Multiset CCA, as an Independent Vector Analysis method or, crucially, as an extension of SRM. In the setting studied here, ShICA improves upon all existing methods.

ShICA inherits from all the advantages of MultiView ICA. Unlike CanICA or ConcatICA (see section 4.2.2.1), it optimizes a proper likelihood and unlike tensorial methods [START_REF] Christian | Tensorial extensions of independent component analysis for multisubject fMRI analysis[END_REF] or SRM, it does not assume any structure on the mixing matrices.

Compared to the likelihood based method of Guo presented in section 4.2.2.2 or to MultiView ICA, ShICA allows different subjects to have different noise variances. This last point is crucial as it allows to separate Gaussian components. In addition, the estimation in ShICA-ML relies on a very efficient closed form E-step. This differs from the likelihood based method of Guo that does not have a closed form E-step and has to rely on a first order approximation.

ShICA-J uses the Multiset CCA presented in section 4.2.1 which is one of the fastest as it reduces to solving a generalized eigenvalue problem. The fact that CCA solves a well defined probabilistic model has first been studied in [START_REF] Francis | A probabilistic interpretation of canonical correlation analysis[END_REF] where it is shown that CCA is identical to multiple battery factor analysis [START_REF] Michael | Factor analysis of multiple batteries by maximum likelihood[END_REF] (restricted to 2 views). This latter formulation differs from our model in that the noise is added on the sensors and not on the components which makes the model unidentifiable. Identifiable variants and generalizations can be obtained by imposing sparsity on the mixing matrices such as in [START_REF] Archambeau | Sparse probabilistic projections[END_REF][START_REF] Klami | Group factor analysis[END_REF][START_REF] Daniela | Extensions of sparse canonical correlation analysis with applications to genomic data[END_REF] or non-negativity [START_REF] Deleus | Functional connectivity analysis of fMRI data based on regularized multiset canonical correlation analysis[END_REF]. Lastly, the work in [START_REF] Li | Joint blind source separation by multiset canonical correlation analysis[END_REF] exhibits a set of sufficient (but not necessary) conditions under which a well defined model can be learnt by the formulation of Multiset CCA used in ShICA-J. The set of conditions we exhibit in this work are necessary and sufficient. We further emphasize that basic Multiset CCA provides a poor estimator, as explained in section 9.2.2.

Let us highlight that ShICA can be seen as a particular instance of IVA where subject specific components s ij are such that s ij = s j + n ij . However, current implementations such as IVA-L [START_REF] Lee | Independent vector analysis (IVA): multivariate approach for fMRI group study[END_REF], IVA-G [START_REF] Anderson | Joint blind source separation with multivariate Gaussian model: Algorithms and performance analysis[END_REF], IVA-L-SOS [START_REF] Bhinge | Extraction of time-varying spatiotemporal networks using parameter-tuned constrained IVA[END_REF], IVA-GGD [START_REF] Anderson | Independent vector analysis: Identification conditions and performance bounds[END_REF] or IVA with Kotz distribution [START_REF] Anderson | Independent vector analysis, the Kotz distribution, and performance bounds[END_REF] estimate view-specific components but shared components are not modeled explicitly. Studying the components post IVA can give great insights about which components are shared and by which subjects (see [START_REF] Long | Independent vector analysis for common subspace analysis: Application to multi-subject fMRI data yields meaningful subgroups of schizophrenia[END_REF]) and individual responses that are closed enough could be merged to produce a shared response. However, the flexibility in modeling different components from each subject comes at the cost of not having a well principled way of aggregating the components into a common shared response. In contrast, ShICA specifically enables extraction of shared components from the subject specific components via its minimum mean squared error estimate.

The IVA theory provides global identifiability conditions in the Gaussian case (IVA-G) [START_REF] Vía | Joint blind source separation from second-order statistics: Necessary and sufficient identifiability conditions[END_REF] and local identifiability conditions in the general case [START_REF] Anderson | Independent vector analysis: Identification conditions and performance bounds[END_REF] from which local identifiability conditions of ShICA could be derived. However, in this work, we provide global identifiability conditions for ShICA. Lastly, IVA can be performed using joint diagonalization of cross covariances [START_REF] Congedo | Orthogonal and non-orthogonal joint blind source separation in the least-squares sense[END_REF][START_REF] Li | Joint blind source separation by generalized joint diagonalization of cumulant matrices[END_REF] although multiple matrices have to be learnt and cross-covariances are not necessarily symmetric positive definite, which makes the algorithm slower and less principled.

Let us point out that extracting a shared response from multiple dataset is also the goal of SRM. Some deep variants [START_REF] Po-Hsuan Chen | A convolutional autoencoder for multi-subject fMRI data aggregation[END_REF] release the orthogonality constrain but they are much more computationally demanding. ShICA leverages ICA theory to provide a much more powerful model of shared responses.

limitations The main limitation of this work is that the model does not reduce the dimension inside each view. In line with other methods, such view-specific dimension reduction has to be done by some external method, typically view-specific PCA. Using specialized methods for the estimation of covariances should also be of interest for ShICA-J, where it only relies on sample covariances. Finally, ShICA-ML uses a simple model of a super-Gaussian distribution, while modeling the non-gaussianities in more detail in ShICA-ML should improve the performance.

conclusion

In this chapter, we introduced the ShICA model as a principled unifying solution to the problems of shared response modelling and GroupICA. ShICA is able to use both the diversity of Gaussian variances and non-Gaussianity for optimal estimation. We presented two algorithms to fit the model: ShICA-J, a fast algorithm that uses noise diversity, and ShICA-ML, a maximum likelihood approach that can separate Gaussian and non-Gaussian components. ShICA algorithms come with principled procedures for shared components estimation, as well as adaptation and estimation of noise levels in each view (subject) and component. In the next chapter, we evaluate ShICA using both real and synthetic data and compare with competitive approaches.

S H A R E D I C A I N P R A C T I C E

In the previous chapter, we have introduced the ShICA model, a principled unifying solution to the problems of shared response modeling and GroupICA. In this chapter we evaluate its practical utility on synthetic data and on brain imaging data. i are obtained by sampling from a uniform density between 0 and 1. Non-Gaussian components are generated from a Laplace distribution and their noise standard deviations are equal. We study 3 cases where either all components are Gaussian, all components are non-Gaussian or half of the components are Gaussian and half are non-Gaussian. We vary the number of samples n between 10 2 and 10 5 and display in Fig 10 .1 the mean Amari distance across subjects as a function of n. The experiment is repeated 100 times using different seeds. We report the median result and error bars represent the first and last deciles. When all components are Gaussian (Fig. 10.1 (a)), CanICA cannot separate the components at all. In contrast ShICA-J, ShICA-ML and Multiset CCA are able to separate them, but Multiset CCA needs many more samples

109

to reach the same Amari distance as ShICA-J or ShICA-ML, which shows that correcting for the rotation due to sampling noise improves the results. Looking at error bars, we also see that the performance of Multiset CCA varies quite a lot with the random seeds: this shows that depending on the sampling noise, the rotation can be very different from identity. MultiViewICA does achieve separation but obtains relatively poor separation performance compared to ShICA or Multiset CCA. When none of the components are Gaussian (Fig. 10.1 (b)), only CanICA, ShICA-ML and MultiView ICA are able to separate the components, as other methods do not make use of non-Gaussianity. Finally, in the hybrid case (Fig. 10.1 (c)), ShICA-ML is able to separate the components very well as it can make use of both non-Gaussianity and noise diversity. As we can see, MultiView ICA yields decent performance though uniformly worse than ShICA-ML. Also note that error bars are very large showing that for some seeds it gives poor results. Overall, MultiView ICA is a lot less reliable than ShICA-ML.

Separation performance in function of non-Gaussianity

We generate data according to model (9.1). Components s are generated using s j = d(x) with d(x) = x|x| α-1 and x ∼ N(0, 1). We impose noise diversity: the noise of view i has standard deviation Σ

1 2
i (obtained by sampling from a uniform density between 0 and 1). Mixing matrices A i are generated by sampling their coefficients from a standardized Gaussian law. The number of samples is fixed to n = 10 5 and we vary α between 0.8 and 2. Each experiment is repeated 40 times using different seeds in the random number generator. We use p = 4 components and m = 5 views. We display in Fig 10 .2 the mean Amari distance across subjects. The experiment is repeated 100 times using different seeds. We report the median result and error bars represent the first and last deciles. When α is close to 1 (components are almost Gaussian), ShICA-J, ShICA-ML and multiset CCA can separate components well (but multiset CCA reaches higher Amari distance than ShICA). In this regime, MultiViewICA yields much higher Amari distance than ShICA-J, ShICA-ML or Multiset CCA but is still better than CanICA which cannot separate components at all. As non-Gaussianity (α) increases, ICA based methods yield better results but ShICA-ML yields uniformly lower Amari distance.

Computation time

We generate components using a slightly super Gaussian density: s j = d(x) with d(x) = x|x| 0.2 and x ∼ N(0, 1). We vary the number of samples n between 10 2 and 10 4 . We compute the mean Amari distance across subjects and record the computation time. The experiment is repeated 40 times. We plot the Amari distance as a function of the . Each point corresponds to the Amari distance/computation time for a given number of samples and a given seed. We then consider for a given number of samples, the median Amari distance and computation time across seeds and plot them in the form of a thick line. From Fig 10 .3, we see that ShICA-J is the method of choice when speed is a concern while ShICA-ML yields the best performance in terms of Amari distance at the cost of an increased computation time. The thick lines for ShICA-J and Multiset CCA are quasi-flat, indicating that the number of samples does not have a strong impact on the fitting time as these methods only work with covariances. On the other hand CanICA or MultiviewICA computation time is more sensitive to the number of samples.

experiments on brain imaging data

The shape of the input data (v, t, m) depend on the dataset used and the experiment performed. A description of the datasets used is available in section 3.1.5 for fMRI data and in section 3.2.5 for MEG data. The values for p (the number of components) is given in the descriptions of the experiments. Note that in the special case of the experiments on the sherlock fMRI dataset, the experimental pipeline is available on github https://github.com/hugorichard/multiviewica/ tree/master/real _ data _ experiments (including downloading and preprocessing of the data).

Robustness w.r.t intra-subject variability in MEG

In the following experiments we consider the Cam-CAN dataset [START_REF] Taylor | The Cambridge Centre for Ageing and Neuroscience (Cam-CAN) data repository: structural and functional MRI, MEG, and cognitive data from a cross-sectional adult lifespan sample[END_REF]. We use the magnetometer data from the MEG of m = 100 subjects chosen randomly among 496. Each subject is repeatedly presented three audio-visual stimuli. For each stimulus, we divide the trials into two sets and within each set, the MEG signal is averaged across trials to isolate the evoked response. This procedure yields 6 chunks of individual data (2 per stimulus). We study the similarity between shared components corresponding to repetitions of the same stimulus. This gives a measure of robustness of each ICA algorithm with respect to intra-subject variability. Data are first reduced using a subjectspecific PCA with p = 10 components. Algorithms are run 10 times with different seeds on the 6 chunks of data, and shared components are extracted. When two chunks of data correspond to repetitions of the same stimulus they should yield similar components. For each component and for each stimulus, we therefore measure the 2 distance between the two repetitions of the stimulus. This yields 300 distances per algorithm that are plotted on Fig 10.4.

The components recovered by ShICA-ML have a much lower variability than other approaches. The performance of ShICA-J is competitive with Multiview ICA while being much faster to fit. Multiset CCA yields satisfying results compared with ShICA-J. However we see that the number of components that do not match at all across trials is greater in Multiset CCA.

The mixing operators in ShICA define spatial maps. In appendix B.5, we plot the average spatial maps across subjects. We use the Elektra Phantom MEG dataset described in 3.2.5 where dipoles at m = 32 different locations emit the same signal. We reduce the data by applying view specific PCA with k = 20 components and algorithms are applied on the reduced data. We select the component that is closer to the true one and compute the L2 norm between the predicted component and the true one after normalization. Then we attempt to recover the position of each dipole by performing dipole fitting on the mixing operator of each view (using only the column corresponding to the true component). The localization error is defined as the mean l2 distance between the true localization and the predicted localization where the mean is computed across dipoles. Each epoch corresponds to 301 samples and 20 epochs are available in total. We vary the number of epochs between 2 and 18 and display in Fig 10 .5 the reconstruction error and the localization error as a function of the number of epochs used. ShICA-ML outperforms other methods. ShICA-J gives satisfying results while being much faster.

MEG Phantom Sinusoidal components

We reproduce the phantom experiment presented in section 8.4. ShICA-ML outperforms other methods. ShICA-J gives satisfying results while being much faster.

Reconstructing the BOLD signal of missing subjects

We reproduce the experiental pipeline described in section 8.3.1. ShICA-ML yields the best R 2 score in all datasets and for any number of components. ShICA-J yields competitive results with respect to Multiview ICA while being much faster to fit.

conclusion

In this chapter, we have shown the practical benefits of ShICA. On simulated data, ShICA clearly outperforms all competing methods in terms of the trade-off between statistical accuracy and computation time. On brain imaging data, ShICA gives more stable decompositions for comparable computation times, and more accurately predicts the data of one subject from the data of other subjects, making it a good candidate to perform transfer learning. As ShiCA only involves lin-ear transforms, decisions based on its output are easier to interpret, making it accessible to practitioners.

In the next section, we show that ICA can be used to perform data augmentation in fMRI.

Part

V C O N D I C A 11 C O N D I C A T H E O R Y
In this chapter, we present an ICA-based method to achieve data augmentation for fMRI data.

Advances in computational cognitive neuroimaging research are related to the availability of large amounts of labeled brain imaging data, since classifiers used to decode brain maps have large samplecomplexity. However, such data are scarce.

To tackle this problem, data generation is an attractive approach, as it could potentially compensate for the shortage of data. Conditional Generative Adversarial Networks (CGANs) are promising generative models [START_REF] Goodfellow | Generative adversarial nets[END_REF] designed for computer vision. However, such improvements have not yet carried over to brain imaging. A likely reason is that CGANs are ill-suited to the noisy, high-dimensional and small-sample data available in functional neuroimaging. Furthermore the training of CGANs is notoriously unstable and there are many hyper-parameters to tune.

In this work, we introduce Conditional ICA: a novel data augmentation technique using ICA together with conditioning mechanisms to generate surrogate brain imaging data and improve image classification performance. Conditional ICA benefits from the abundant resting state data and can be trained with only few labeled samples.

Figure 11.1: Conditional ICA approach. Our method aims to generate surrogate data from Task and Rest fMRI data by synthesizing statistical maps that qualitatively fit the distribution of the original maps. These can be used to improve the accuracy of machine learning models that identify contrasts from the corresponding brain activity patterns.

methods

Spatial Dimension reduction

The outline of the proposed approach is presented in Fig. 11.1. While brain maps are high-dimensional, they span a smaller space than that of the voxel grid. For the sake of tractability, we reduce the dimension of the data by projecting the voxel values on the highresolution version of the Dictionaries of Functional Modes DiFuMo atlas [START_REF] Dadi | Fine-grain atlases of functional modes for fMRI analysis[END_REF], i.e. with p = 1024 components. The choice of dimension reduction technique generally has an impact on the results. However, we consider this question to be out of the scope of the current study and leave this to future work.

A generative model for task data

Consider a task dataset X task in R p,n where n is the number of observations (samples) and p = 1024 the number of components in the atlas. X task can be seen as n observations of a random vector x task ∈ R p . Let us consider how to learn the distribution of x task . Assuming a Gaussian distribution is standard in this setting, yet, as shown later, it misses key distributional features. Moreover, we consider a model that subsumes the distribution of any type of fMRI data (task or rest): a linear mixture of k p independent temporal signals. We therefore use temporal ICA to learn a dimension reduction and unmixing matrix W task ∈ R k,p such that the k components of W task x task are as independent as possible.

A straightforward method to generate new task data would be to independently sample them from the distribution of the components. This is easy because such distribution has supposedly independent marginals. We apply an invertible quantile transform q task to the components of W task x task so that the distribution of q task (W task x task ) has standardized Gaussian marginals. Since it also has independent marginals, it is given by N(0 k , I k ) from which we can easily sample.

However task datasets have a small number of samples (10 ∼ 10 2 ). As a result, there are too few samples to learn a high quality unmixing matrix. In contrast, resting state datasets have a large number of samples (10 4 ∼ 10 5 ). Therefore, we replace the unmixing matrix learned on task data W task by the unmixing matrix learned on resting state data W rest .

We form z task = W rest x task and learn its quantile transform q. The encoding model is thus given by:

z task = q(W rest x task ) (11.1)
However, the independence assumption no longer holds and thus a latent structure among the marginals of z task has to be taken into Figure 11.2: Conditional ICA approach in depth. The approach proceeds by learning a temporal ICA of rest data X rest ∈ R p,n , resulting in independent components and unmixing matrix W rest ∈ R k,p . Applying the unmixing matrix to the task data, we obtain samples in the component space that we map to a normal distribution, yielding Z task ∈ R k,n . Then, we estimate the covariance Λ ∈ R k,k (all classes are assumed to have the same covariance) and the class-specific means µ 1 , . . . , µ C ∈ R k according to Ledoit-Wolf's method. For each class c, we can draw random samples Ztask c ∈ R k,n fakes from the resulting multivariate Gaussian distribution N(µ c , Λ) and obtain fake data Xc ∈ R p,n fakes by applying the inverse quantile transform and re-mixing the data using the pseudo inverse of the unmixing matrix. We append these synthetic data to the actual data to create our new augmented dataset on which we train classifiers.

account. In addition the generative model needs to be conditioned to each class. We therefore assume that the samples in class c, x task c are such that:

q(W rest x task c ) ∼ N(µ c , Λ) (11.2)
In order to maximize the number of samples used to learn the parameters of the model, we assume that the quantile transform q and the latent covariance Λ do not depend on the class c. However, the mean µ c , that can be learned efficiently using just a few tens of samples, depends on class c. Λ is learned using all task samples from a standard shrunk covariance estimator Λ = Σ(1 -α) + α k tr(Σ)I k where α is given by the Ledoit-Wolf formula [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF] and Σ is the sample covariance of z task .

The generative model of data for brain maps in a certain class c is given by the pseudo inverse of the encoding model:

x c = (W rest ) † q -1 ( ) (11.3)
with ∼ N(µ c , Λ) and (W rest ) † is the Moore Penrose inverse of W rest An overview of our generative method is shown in Fig. 11.2.

related work

In image processing, data augmentation is part of standard toolboxes and typically includes operations like cropping, rotation, translation. On fMRI data these methods do not make much sense as brain data are not invariant to such transformations. More advanced techniques [START_REF] Zhuang | Fmri data augmentation via synthesis[END_REF] are based on generative models such as CGANs or variatonal autoencoders [START_REF] Diederik | Auto-encoding variational bayes[END_REF]. Although CGAN-based method are powerful they are slow and difficult to train [START_REF] Arjovsky | Wasserstein GAN[END_REF]. Our method is not an adversarial procedure, however it relates to other powerful generative models such as variatonal autoencoders [START_REF] Diederik | Auto-encoding variational bayes[END_REF] with which it shares strong similarities. Indeed the analog of the encoding function in the variational auto-encoder is given by e(x) = Λ -1 2 q(W rest x) in our model and the analog to the decoding function in the variational auto-encoder is given by d(z) = (W rest ) † q -1 (Λ 1 2 z) in our model. As in the variational autoencoder, e approximately maps the empirical data distribution to a standardized Gaussian distribution, while the reconstruction error defined by the difference in l2 norm d(e(x)) -x 2 2 must remain small. Lastly, another classical generative model related to ours is normalizing flows. We note that when W rest is square (no dimension reduction in ICA), the decoding operator d is invertible (its inverse is e) making our model an instance of normalizing flows [START_REF] Jimenez | Variational inference with normalizing flows[END_REF]. A great property is thus the simplicity and reduced cost of data generation.

conclusion

In this chapter, we introduced Conditional ICA, a fast generative model for task data. Conditional ICA is essentially a linear generative model with pointwise non-linearity, which makes it cheap, easy to instantiate on new data, and to introspect. In the next chapter, we look at the performance of Conditional ICA on fMRI data.

In the previous chapter, we have introduced Conditional ICA, an efficient generative model for task fMRI data. In this chapter, we benchmark Conditional ICA as a generative model of task data against various augmentation methods by assessing their ability to improve classification accuracy on a large task fMRI dataset. We find that Conditional ICA yields highest accuracy improvements. In particular, Conditional ICA outperforms conditional GANs [START_REF] Mirza | Conditional generative adversarial nets[END_REF] while being much easier to optimize and interpret. Lastly, we show on 8 different datasets that the use of Conditional ICA results in systematic improvements in classification accuracy ranging from 1% to 5%.

dataset, data augmentation baselines and classifiers used

The unmixing matrices are learned on the rest HCP dataset [START_REF] David C Van Essen | The WU-Minn human connectome project: an overview[END_REF] using 200 subjects. These data were used after standard preprocessing, including linear detrending, band-pass filtering ([0.01, 0.1]Hz) and standardization of the time courses. The other 8 datasets [111, 121-123, 127, 137, 149] are obtained from the Neurovault repository [START_REF] Krzysztof | NeuroVault. org: a web-based repository for collecting and sharing unthresholded statistical maps of the human brain[END_REF]. The classes used in each of these datasets correspond to the activation maps related to contrasts (such as "face vs tools") present in the set of tasks of each dataset. In table 12.1, we give references to the datasets used as well as the total number of samples (subjects), the size of train and test sets in each of the cross validation splits and the number of classes in each dataset. We consider 4 alternative augmentation methods: ICA, Covariance, ICA + Covariance and CGANs. When no augmentation method is applied, we use the Original label.

The ICA method applies ICA to X task to generate unmixing matrices W task and components S task = W task X task . To generate a sample xc from class c, we sample independently from each component restricted to the samples of class c yielding stask The ICA + Covariance method combines the augmentation methods ICA and Covariance: samples are drawn following the ICA approach, CGANs can generate fake data from a given class. In the CGAN method, the generator and discriminator have a mirrored architecture with 2 fully connected hidden layer of size (256 and 512). The number of epochs, batch size, momentum and learning rate are set to 20k, 16, 0.9, 0.01 and we use the Leaky RELU activation function.

We evaluate the performance of augmentation methods through the use of classifiers: logistic regression (LogReg), linear discriminant analysis with Ledoit-Wold estimate of covariance (LDA), perceptron with two hidden layers (MLP) and random forrests (RF). The hyperparameters in each classifier are optimized through an internal 5-Fold cross validation. We set the number of iterations in each classifier so that convergence is reached. The exact specifications are given in table 12.2.

comparing classification accuracy gains on task hcp dataset

In order to compare the different augmentation methods, we measure their relative benefit in the context of multi-class classification. We use 787 subjects from the HCP task dataset that contains 23 classes and randomly split the dataset into a train set that contains 100 subjects and a test set that contains 687 subjects. In each split, we run augmentation methods on the train set to generate fake samples corresponding to 200 subjects. These samples are then appended to the train set, resulting in an augmented train set on which the classifiers are trained. The results displayed in table 12.3 show that Conditional ICA always yields a higher accuracy than when no augmentation method is applied. The gains are over 1% on all classifiers tested excepts with the random forest classifier which yields much lower accuracy than other methods. By contrast, ICA+Covariance and ICA lead to a decrease in accuracy while the Covariance approach leads to non-significant gains.

In In practice, the computational over-head induced by Conditional ICA is negligible. methods, the computational over-head induced by CondICA is very low.

12.3 gains in accuracy brought by conditional ica on eight datasets.

In this experiment, we assess the gains brought by Conditional ICA data augmentation on the eight different task fMRI datasets refered to in section 12.1. The experimental pipeline is exactly the same as with the HCP task dataset. We report in Fig. 12.1 the cross-validated accuracy of classifiers with and without augmentation. We notice that the effect of data augmentation is consistent across datasets, classifiers and splits, with 1% to 5% net gains. Lastly, we provide a sensitivity analysis on the number of components used in CondICA in figure 12.2. CondICA gives good performance for numbers of components between 800 and 1000 components. In all experiments we used k = 900 components.

conclusion

When Conditional ICA is used as a data augmentation method, it yields consistent improvement in classification accuracy: on 8 tasks fMRI datasets, we observe an increase in accuracy between 1% and 5% depending on the dataset and the classifier used. Importantly, this performance was obtained without any fine-tuning of the method, showing its reliability. One can also notice that our experiments cover datasets with different cardinalities, from tens to thousand, and different baseline prediction accuracy.

The systematic performance improvement CondICA yields makes it a promising candidate for data augmentation in a wide range of contexts. Future work may focus on its applicability to other decoding tasks such as the diagnosis of Autism Spectrum Disorder (ASD) [START_REF] Nicha C Dvornek | Identifying autism from resting-state fMRI using long short-term memory networks[END_REF][START_REF] Eslami | ASD-DiagNet: a hybrid learning approach for detection of autism spectrum disorder using fMRI data[END_REF][START_REF] Eslami | Auto-ASD-network: a technique based on deep learning and support vector machines for diagnosing autism spectrum disorder using fMRI data[END_REF] or Attention-Deficit/Hyperactivity Disorder detection (ADHD) [START_REF] Mao | Spatio-temporal deep learning method for adhd fmri classification[END_REF]. Other extensions of the present work concern the adaptation to individual characteristics (e.g. age) prediction where fMRI has shown some potential.

An optimal transport approach to hyperalignment

In this work, we benchmark optimal transport, ridge regression and scaled Procrustes to align the data of two subjects. Optimal Transport and Ridge regression outperformed alternatives in that task. The implementation of the methods developed in this thesis is freely available on Github https://github.com/hugorichard. Some of the code we wrote has made its way to bigger packages.

Mvlearn

Mvlearn [START_REF] Perry | mvlearn: Multiview Machine Learning in Python[END_REF] is a Python package for multiview learning tools. It offers reference implementations for algorithms and methods related to multiview learning. Its API is close to the scikit-learn [START_REF] Abraham | Machine learning for neuroimaging with scikit-learn[END_REF] one, making it easy to learn. We have implemented the GroupICA, GroupPCA and MultiViewICA modules of mvlearn.

Brainiak

Brainiak [START_REF] Kumar | BrainIAK: The brain imaging analysis kit[END_REF][START_REF] Kumar | BrainIAK tutorials: Userfriendly learning materials for advanced fMRI analysis[END_REF], is a Python package that applies machine learning methods to neuroimaging data. Its API is the same as in scikit-learn and it includes modules such as Representational Similarity Analysis or Shared response modeling. We have implemented the FastSRM module of Brainiak.

conclusion

In this thesis, we have presented three methods to perform component analysis of multi-subject neuroimaging data and a data augmentation method for fMRI data.

First, in chapter 5 and chapter 6, we have developed an atlas based procedure that is shown to accelerate significantly the existing procedure for performing dimension reduction of fMRI data in a multi-subject context with provably no loss of performance. It is now possible to apply these algorithms in big datasets where the number of subjects is of the order of several hundreds, with several thousand samples and several hundred thousand features.

Then, we have proposed in chapter 7 and chapter 8 a novel unsupervised algorithm, MultiViewICA, that reveals latent sources observed through different views. Using an independence assumption, we have demonstrated that the model is identifiable, provided that the latent sources are not Gaussian. In contrast to previous approaches, the proposed model leads to a closed-form likelihood, which we then optimize efficiently using a dedicated alternate quasi-Newton approach. Therefore, MultiViewICA enjoys the statistical guarantees of maximum-likelihood theory, while still being tractable. MultiViewICA outperforms other unsupervised methods used to process fMRI and MEG data in the context of shared response modeling. However, it assumes the same level of noise in all subjects, which does not model properly between-subjects variability.

In chapter 9 and chapter 10, we have extended MultiViewICA in orter to deal with source noise heteroscedasticity. In practice ShICA outperforms MultiViewICA and other unsupervised methods used in the context of shared response modeling.

Lastly, we have introduced in chapter 11 and chapter 12 a data augmentation method based on ICA that outperforms deep learning algorithms in terms of decoding accuracy while being much faster.

future work and perspectives

Combined with the FastSRM algorithm, MultiViewICA and ShICA yield a novel way to make use of multi-view data. They yield a set of operators per view that map the data of each view to a shared response, reducing the variability between views. In principle, reducing the variability between views should facilitate the understanding of the data and therefore increase the performance of classification based tasks such as contrast maps labeling, automatic diagnosis or age prediction. Investigating to which extend these benefits are observed could be the topic of future research.

A second practical direction would be to apply our methods to different neuroimaging settings in which assumptions differ from ours. This thesis is geared towards naturalistic imaging, where the temporal response is assumed to be shared across subjects. We see at least two other neuroimaging settings in which our methods can be useful. The first one is the analysis of resting state data assuming that spatial topographies are shared across subjects. In this setting, the spatial topographies become the common components while the mixing operators correspond to a set of time-courses. In practice, transposing the data is enough to enforce such assumptions. A second one is the analysis of MEG / EEG data assuming both a spatial and temporal mixing. This can be done in practice by stacking the features of consecutive samples.

In terms of methods, we have treated separately dimension reduction (chapter 5) and source identification (chapter 7 and chapter 9). Future work might focus on understanding how these two steps can be performed jointly. Another possible extension in the case of naturalistic stimuli, is to assume that mixing matrices are close to each other. Indeed, as such matrices represent spatial topographies, such prior makes sense. Lastly, our data augmentation method does not make use of our understanding of multiview datasets. This constitutes an exciting direction of research. We fix a subject i. Since s has independent components, so does s + n i . Following [START_REF] Comon | Independent component analysis, a new concept?[END_REF], Theorem 11, there exists a scale-permutation matrix P i such that A i = A i P i . As a consequence, we have s + n i = P i (s + n i ) for all i.

Then, we focus on subject 1 and subject i = 1:

s + n 1 -(s + n i ) = P 1 (s + n 1 ) -P i (s + n i ) (A.1) n 1 -n i = P 1 (s + n 1 ) -P i (s + n i ) (A.2) ⇐⇒ P 1 s -P i s = P i n in i + n 1 -P 1 n 1 (A.3)

Since the right hand side of equation (A.3) is a linear combination of Gaussian random variables, this would imply that P 1 s -P i s is also Gaussian. However, given that s is assumed to be non-Gaussian, the equality can only hold if P 1 = P i and both the right and the left hand side vanish. Therefore, the matrices P i are all equal, and there exists a scale and permutation matrix P such that A i = A i P.

a.1.2 Proof of Prop. 13

We consider W i = Λ(A i ) -1 , where Λ is a diagonal matrix. We recall x i = A i (s + n i ), so that y i = W i x i = Λ(s + n i ). The gradient of L is given by equation (7.9): ]. One the one hand, we have G(0) = -1. On the other hand, if we assume for instance that f has sub linear growth (i.e. |f (x)| c|x| α + d for some α < 1) or that φ is positive, we find that G(+∞) = +∞. Therefore, G cancels, which concludes the proof.

G i = 1 m E[f ( s)(s + n i ) ]Λ + 1 -1/m σ 2 ΛE[   n i - 1 m -1 j =i n j   (s + n i ) ]Λ -I p = 1 m E[f (Λ(s + 1 m j n j ))(s + n i ) ]Λ + σ 2 (1 -1/m) σ 2 Λ 2 -I p (A.

a.1.3 Stability conditions

We consider W i = Λ(A i ) -1 where Λ is such that the gradients G i all cancel. We consider a small relative perturbation of W i of the form W i ← (I p + E i )W i , and consider the effect on the gradient. We define ∆ i = G i (I p + E 1 )W 1 , . . . , (I p + E m )W m . Denoting C = 1-1/m σ 2 and ñ = 1 m m i=1 n i , we find:

∆ i = ∆ i 1 + C∆ i 2 -I p where ∆ i 1 = E[ 1 m f   Λ(s + ñ) + 1 m m j=1 E j Λ(s + n j )   (s + n i ) Λ(I p + E i ) ] (A.5)
and

∆ i 2 = E[   Λn i - 1 m -1 j =i Λn j + E i Λ(s + n i ) - 1 m -1 j =i E j Λ(s + n j )   (s + n i ) Λ(I p + E i ) ] (A.6)
The first term is expanded at the first order, denoting S = m j=1 E j :

∆ i 1 = E[ 1 m f (Λ(s + ñ)) 1 m m j=1 E j Λ(s + n j ) + f (Λ(s + ñ)) (s + n i ) Λ(I p + E i ) ] = E[ 1 m f (Λ(s + ñ))(s + n i ) Λ(I p + E i ) + 1 m 2 S f (Λ(s + ñ))(s 2 ) Λ 2 + 1 m 2 E i f (Λ(s + ñ))((n i ) 2 ) Λ 2 ] (A.7)
The symbol denotes the element-wise multiplication, f (s) =     f (s 1 )

. . .

∆ i 2 = E[σ 2 Λ 2 (I p + E i ) + (1 + σ 2 )E i Λ 2 - 1 m -1 (S -E i )Λ 2 ] (A.8)
Combining this, we find:

∆ i = (E i ) + E i Γ E + S Γ S (A.9)
where

Γ E = 1 m 2 E[f (Λ(s + ñ))((n i ) 2 ) ] + (1 - 1 m ) σ 2 σ 2 + 1 σ 2 Λ 2 Γ S = 1 m 2 E[f (Λ(s + ñ))(s 2 ) ] - 1 mσ 2 Λ 2
are p × p matrices, independent of the subject. This linear operator is the Hessian block corresponding to the i-th subject: Denoting H the Hessian, it is the mapping H(E 1 , . . . , E m ) = (∆ 1 , . . . , ∆ m ).

The coefficient ∆ i ab only depends on (E i ab , E i ba , E 1 ab , . . . , E m ab ). Therefore, the Hessian is block diagonal with respect to the blocks of coordinates (E

              
The positivity of H is equivalent to the positivity of this operator for all pairs a, b. We now assume ββ > 0.

First, we should note that K(ε, ε , β, β ) is congruent to Since I m ⊗ M α and γ1 ⊗ I 2 commute, the minimum value of Sp(K) is min(I m ⊗ M α ) + min(γSp(1)) = 1 2 (α + α -(α -α ) 2 + 4) + m min(0, γ). Since we assumed ββ > 0 we have γ > 0. This is similar to the usual ICA case, we find that the condition is αα > 1.

K(ε β β , ε
If the following conditions hold for all pair of components a, b, the components are a local minimum of the cost function:

• Γ S ab Γ S ba 0 • Γ E ab Γ E ba > 1 a.1.

Reproducing time-segment matching experiment

We reproduce the time-segment matching experiments described in [START_REF] Po-Hsuan Chen | A convolutional autoencoder for multi-subject fMRI data aggregation[END_REF] and [START_REF] Zhang | A searchlight factor model approach for locating shared information in multi-subject fMRI analysis[END_REF] and use two fold classification over runs instead of 5fold as we have done in chapter 8. We used the sherlock data available at http://arks.princeton.edu/ark:/88435/dsp01nz8062179 and the full brain mask provided in the Python package associated with the paper. We applied high-pass filtering (140 s cutoff) and the time series of each voxel were normalized to zero mean and unit variance.

The results are available in

a.2 related work

The following table describes some usual method for extracting shared components from multiple subjects datasets. The column "Modality/-Components" describes the type of data for which each algorithm was initially proposed, even though each algorithm could be applied on any type of data. The components type can be either temporal if extracted components are time courses or spatial if they are spatial patterns. Proof. The proof is in two parts. First, we show that there exist some coefficients α 1 , . . . , α m such that the matrix i α i Σ i has distinct coefficients on the diagonal. Then, since we have i α i Σ i = P ( i α i Σ i ) P , and the diagonal i α i Σ i has distinct entries, we can invoke the unicity of the eigenvalue decomposition for symmetric matrices, which shows that P is necessarily a permutation matrix. Now, the only thing left is to prove is that Assumption 2 implies the existence of this linear combination.

Method

We assume by contradiction that any linear combination of the Σ i has two equal entries.

For α = [α 1 , . . . , α m ], we let S(α) = diag( i α i Σ i ) ∈ R p , where diag(•) extracts the diagonal entries. The operator S is linear. We now define for j, j p the linear form jj (α) = S(α) j -S(α) j ∈ R. The assumption on the linear combinations of Σ i simply rewrites: For all α ∈ R m , there exists j, j p such that jj (α) = 0. From a set point of view, this relationship writes j,j Ker( jj ) = R m .

Since the jj are all linear forms, the Ker( jj ) are subspaces of dimensions m or m -1, and since their union is of dimension m, there exists j, j such that Ker( jj ) = R m , i.e. such that jj = 0.

As a consequence, we have for all α, S(α) j = S(α) j . This implies that the sequences (Σ ij ) i and (Σ ij ) i are equal, which contradicts Assumption 2.

We have therefore shown that Assumption 2 implies the existence of a linear combination of the Σ i that has distinct entries, which concludes the proof. Proof. We sort the eigenvectors in p groups of m vectors so that all vectors in group l have their l-th coordinate different from 0. Let z (l) be an eigenvector in group l and let us call w l ∈ R m the non-zero coordinates of this eigenvector: ∀i ∈ {1 . . . m}, w li = z (l) l+(i-1)p . We have:

A l w l = B l w l λ l (B.5) We now show that the biggest eigenvalue of (B.5) is strictly above 1 while all others are strictly below 1. The core of the proof comes from the study of the eigenvalues of a matrix modified by a rank 1 matrix. The reasoning we use here follows [START_REF] Golub | Some modified matrix eigenvalue problems[END_REF] (end of section 5). The characteristic polynomial is given by: P(λ) = det((I + K) -1 K -λI + (I + K) -1 uu ) (B.8)

A l =        1 + Σ 1l 1 . . . 1

Let us introduce

∝ det(I + ((I + K) -1 K -λI) -1 (I + K) -1 uu ) (B.9)

where we implicitly focus here on eigenvalues λ such that det((I + K) -1 K -λI) = 0 ⇐⇒ ∀i, λ = k i 1+k i . We then use the following property: Let A ∈ R a,b and B ∈ R b,a we have det(I a + AB) = det(I b + BA).

Let us call χ(λ) = det(I + ((I + K) -1 K -λI) -1 (I + K) -1 uu ) we have: [START_REF] Archambeau | Sparse probabilistic projections[END_REF] where k i = Σ il > 0. Taking the derivative we get

χ(λ) = 1 + u ((I + K) -1 K -λI) -1 (I + K) -1 u (B.10) = 1 + m i=1 1 1 + k i 1 k i 1+k i -λ (B.
χ (λ) = m i=1 1 1 + k i 1 ( k i 1+k i -λ) 2 > 0 (B.12)
Trivially, ∀i, k i 1+k i < 1. We also have

χ(1) = 1 + m i=1 1 1 + k i 1 k i 1+k i -1 = 1 -m < 0 (B.13)
and lim λ→+∞ χ(λ) = 1 so as χ is continuous and strictly increasing on [1, +∞[. Therefore, it reaches 0 only once on this interval (excluding 1 since we know χ(1) = 0). Therefore the greatest eigenvalue λ * is strictly above 1 while all other eigenvalues are strictly below 1.

Note that because χ > 0, λ * is of multiplicity 1. In the analysis above we ignored those eigenvalues λ such that λ = k i 1+k i for some i. However since k i 1+k i < 1, none of these eigenvalues can be the largest one.

Finally, the p first eigenvectors belong to different groups (the corresponding eigenvalues are all strictly above 1). This shows that these eigenvectors have different first non-zero coordinates.

b.2 identifiability results for m < 3

We have a slightly weaker identifiability result when m = 2.

Proposition 24. Let m = 2, and suppose that the scalars (1 + Σ 1j )(1 + Σ 2j ) for j = 1 . . . p are all different. We let Θ = (A 1 , A 2 , Σ 1 , Σ 2 ) that also generates x 1 , x 2 . Then, there exists a permutation and scale matrix P such that A 1 = A 1 P and A 2 = A 2 P -.

Proof. We let P = A -1 1 A 1 . Since C 12 = I p , it holds A -1 2 A 2 = P -. Then, we have I p + Σ 1 = P(I p + Σ 1 )P . This means that there exists U ∈ O p such that P = (I p + Σ 1 ) 1 2 U(I p + Σ 1 ) -1 2 . Since P -(I p + Σ 2 )P -1 = I p + Σ 2 , we find U(I p + Σ 1 )(I p + Σ 2 )U = (I p + Σ 1 )(I p + Σ 2 ). By identification, U is a permutation matrix, and P is a scale and permutation matrix.

As a consequence, when there are only two subjects, it is possible to recover the components and noise levels up to a scaling factor. When there is only one view, m = 1, there is a global rotation indeterminacy: A 1 (I p + Σ 1 )A 1 = A 1 (I p + Σ 1 )A 1 for A 1 = A 1 (I p + Σ 1 ) 1 2 U(I p + Σ 1 ) -1 2 where U is any orthogonal matrix. In this case, we lose identifiability. where c does not depend on Σ i or W i Therefore we get closed-form updates for Σ i :

Σ i ← diag(E[(y i -E[s|x])(y i -E[s|x]) ] + V[s|x]) (B.56)
We update W i by performing a quasi-Newton step. We use the relative gradient G W i and H W i defined by

J(W i + εW i ) = J(W i ) + ε, G W i + 1 2 ε, H W i ε (B.57)
We get: where the Hessian approximation is exact when the unmixed data have truly independent components.

J(W i + εW i ) =
Updates for W i are then given by W i ← (I -ρ( H W i ) -1 G W i )W i , where ρ is chosen by backtracking line-search. We alternate between computing the statistics E[s|x] and Var[s|x] (E-step) and updates of parameters Σ i and W i for i = 1 . . . m (M-step).

160 a. 4 A

 4 Average forward operators on fMRI datasets 161 a.5 Synthetic benchmark using additive noise on the sensors 162 a.6 Summary of our quantitative results 162 b shica 167 b.1 Lemmas 167 b.2 Identifiability results for m < 3 170 b.3 EM E-step and M-step for ShICA with Gaussian components 171 b.3.1 E-step 171 b.3.2 M-step 171 b.4 EM E-step and M-step for ShICA with non-Gaussian components 171 contents xi b.4.1 E-step 171 b.4.2 M-step 174 b.5 CamCAN spatial maps 176

1

 1 Controlled experiments in cognitive brain imaging

  published work (Spotlight) H. Richard et al. "Modeling Shared responses in Neuroimaging Studies through MultiView ICA." In: Advances in Neural Information Processing Systems 33. Dec. 2020 1.2.3 Shared ICA for neuroimaging data

  published work H. Richard et al. "Shared Independent Component Analysis for Multi-Subject Neuroimaging." In: Advances in Neural Information Processing Systems 33. Dec. 2021 1.2.4 Conditional ICA

  published work Badr Tajini, Hugo Richard, and Bertrand Thirion. "Functional Magnetic Resonance Imaging data augmentation through conditional ICA." In: MICCAI (2021) (Oral, co-first authorship with equal contribution)

θ

  * and from our assumption about J n we get J n (θ 0 ) P ----→ n→∞ J(θ * ) = -I(θ * ).

Figure 3 . 1 :

 31 Figure 3.1: A model of the hemodynamic response function

Figure 3 . 2 :

 32 Figure 3.2:Contrast "sentence listening versus sentence reading" computed from the fMRI data of one subject in the IBC dataset[START_REF] Luísa Pinho | Individual Brain Charting dataset extension, second release of high-resolution fMRI data for cognitive mapping[END_REF][START_REF] Luísa Pinho | Individual Brain Charting, a high-resolution fMRI dataset for cognitive mapping[END_REF]. This contrast was downloaded from neurovault[START_REF] Krzysztof | NeuroVault. org: a web-based repository for collecting and sharing unthresholded statistical maps of the human brain[END_REF] (collection 2138 subject 1 session 0) .

Algorithm 1 : 4 . 2 . 5

 1425 Hyperalignment Input: Data X 1 , . . . , X m ∈ R p×n , number of iterations n iter Select a random subject i ∼ U(1, m) Initialize the alignment operators for s = 1 . . . m do if s = i then P st = I p end else P st = P(X t X s ) end end T = m s=1 P st X s m Main loop for it = 1 . . . n iter do Align data and the current template for s = 1 . . . m do P st = P(T X s ) end Compute the template as the mean of aligned data T = m s=1 P st X s m end return Estimated template T and operators P st The shared response model (SRM)

Figure 4 . 1 :

 41 Figure 4.1: Shared response model: The raw fMRI data are modeled as a weighted combination of subject-specific spatial components with additive noise. The weights are shared between subjects and constitute the shared response to the stimuli.

59

 59

Figure 5 . 1 :

 51 Figure 5.1: FastSRM algorithm In step 1, data x i are projected onto an atlas U i that may depend on the subject i (top). In step 2 a SRM algorithm is applied on reduced data to compute the shared response.

Figure 6 . 1 :

 61 Figure 6.1: Benchmark of SRM algorithms on synthetic data: Performance, fitting time and convergence of SRM algorithms in the deterministic (left) or probabilistic (right) case. As expected, when optimal atlases are used, the performance is the same as if no atlas is used but the fitting time is much lower. This is even more pronounced when the number of iterations is high (and looking at convergence curves, we see that more iterations could be performed to be closer to a stationary point).

Figure 6 . 2 :

 62 Figure 6.2: Benchmark of SRM algorithms on fMRI data (top) Timesegment matching accuracy (middle) Fitting time (bottom) Memory usage.When the optimal atlas is used, the accuracy is the same as when no atlas is used but the fitting time is reduced by a factor 10 to 100 and so is the memory usage.

Figure 6 . 3 :

 63 Figure 6.3: Experiment -Predict age from spatial components extracted using FastSRM:We first learn the spatial components from fMRI data using FastSRM. We learn one Ridge model per spatial components to predict age across subjects. Then, these models are aggregated using a Random Forest (like in[START_REF] Rahim | Joint prediction of multiple scores captures better individual traits from brain images[END_REF]).

Figure 6 . 4 :

 64 Figure 6.4: Age prediction from spatial components: (top) FastSRM predicts age with a better accuracy than chance resulting in a mean absolute error (MAE) of 7.5 years. (middle) FastSRM is more than 200x faster than ProbSRM and uses 10x less memory, hence it scales better than ProbSRM. (bottom)The three most important spatial components in terms of the reduction in Gini impurity they bring (see Gini importance or Feature importance in[START_REF] Breiman | Random forests[END_REF],[START_REF] Louppe | Understanding variable importances in forests of randomized trees[END_REF]). From top to bottom, the most important spatial component (feature importance: 16%) highlights the visual dorsal pathway, the second most important spatial component (feature importance: 12%) highlights the precuneus and the third most important spatial component (feature importance: 8%) highlights the visual ventral pathway.

Figure 6 . 5 :

 65 Figure 6.5: Evolution of the most predictive spatial component with age: (Top) Spatial component most predictive of age averaged within groups of different age (18-35, 36-48, 48-61, 61-74, 74-88). (Bottom) Mean activation in the region highlighted by the mask on the left. We see that the activity in the dorsal pathway decreases with age, which explains why this spatial component is a good predictor of age.

Algorithm 2 :

 2 Alternate quasi-Newton method for MultiView ICA Input: Dataset (x i ) m i=1 , initial unmixing matrices W i , noise parameter σ, function f, tolerance ε Set tol= +∞, s = 1 m p i=1 W i x i while tol> ε do tol = 0 for i = 1 . .. m doCompute y i = W i x i , s-i = s -1 m y i , gradient G i (eq. (7.9)) and Hessian H i (eq. (7.11)) Compute the search direction D

Figure 8 . 1 :

 81 Figure 8.1: Synthetic experiment: reconstruction error of the algorithms on data following model x i = A i (s + n i ).

Figure 8 . 2 :Figure 8 . 3 :

 8283 Figure 8.2: Reconstructing the BOLD signal of missing subjects -Experimental procedure 80% of the runs are used to compute forward operators (spatial components) for every subject (left). Then the forward operators and data from the left-out runs of 80% of the subjects one are used to compute the shared response during the left-out runs. Lastly, the shared response during the left-out runs and the forward operators of the test subjects are used to predict the data of the test subjects during the left-out runs. The performance of the model is measured by comparing the prediction and true data using the R 2 score.

Figure 8 . 4 :

 84 Figure 8.4: Top: Reconstructing the BOLD signal of missing subjects. Mean R 2 score between reconstructed data and true data (higher is better). Bottom: Between subjects time-segment matching. Mean classification accuracy. Error bars represent a 95 % confidence interval over cross validation splits.

Figure 8 . 5 :

 85 Figure 8.5: Reconstructing the BOLD signal of missing subjects: Reconstruction R2 score per voxel We plot for ConcatICA, SRM and Multi-ViewICA, the R2 score per voxel using 50 components for datasets sherlock, forrest, raiders and clips. We can see that data reconstructed by MultiViewICA are more faithful reproduction of the original data than other methods.

Figure 8 . 6 :

 86 Figure 8.6: Effect of the parameter σ: We compute the accuracy of the MultiViewICA pipeline on the time-segment matching experiment for various values of the σ hyperparameter over a grid. The accuracy varies only marginally with σ.

Figure 8 . 7 :

 87 Figure 8.7: Between runs time-segment matching.Interesting components correlates more when they correspond to the same stimulus (same scenes of the movie) than when they correspond to distinct stimuli (different scenes). We extract 20 components and report the mean accuracy of the three best performing components

Figure 8 . 8 :

 88 Figure 8.8: Between-runs time segment matching: spatial maps and timecourses Left: Timecourses of the 3 shared components yielding the highest accuracy. The two displayed set of runs correspond to the same scenes in the movie. Right: Localisation of the same shared components in the brain

AmplitudeFigure 8 . 9 :

 89 Figure 8.9: Top: Experiment on MEG Phantom data. Reconstruction error is the norm of the difference between the estimated and true component. Localization error is the distance between the estimated and true dipole. Middle and Bottom: Experiment on 200 subjects from the CAM-can dataset Middle: Time course of 11 shared components (one color per component). We recover clean evoked potentials. Bottom: Associated brain maps, obtained by averaging component estimates registered to a common reference.

Figure 9 . 1 :

 91 Figure 9.1: Amari distance between true mixing matrices and estimates of Multiset CCA when covariances are perturbed. Different curves correspond to different eigengaps. When the gap is small, a small perturbation can lead to complete mixing. Jointdiagonalization (colored dotted lines) fixes the problem.

10. 1

 1 synthetic experiment In the following synthetic experiments, data are generated according to model (9.1) with p = 4 components and m = 5 views and mixing matrices are generated by sampling coefficients from a standardized Gaussian. 10.1.1 Separation performance: different use cases Gaussian components are generated from a standardized Gaussian and their noise has standard deviation Σ

Figure 10 . 1 :

 101 Figure 10.1: Separation performance: Algorithms are fit on data following model 9.1 (a) Gaussian components with noise diversity (b) Non-Gaussian components without noise diversity (c) Half of the components are Gaussian with noise diversity, the other half is non-Gaussian without noise diversity.

Figure 10 . 2 :Figure 10 . 3 :

 102103 Figure 10.2: Separation performance in function of non-Gaussianity Separation performance of algorithms for sub-Gaussian α < 1 and super-Gaussian α > 1 components

Figure 10 . 4 :

 104 Figure 10.4: Robustness w.r.t intra-subject variability in MEG: (top) 2 distance between shared components corresponding to the same stimuli in different trials. (bottom) Fitting time.

3

 3

10. 2 . 2

 22 MEG Phantom experiment 10.2.2.1 Elekta Phantom

Figure 10 . 5 : 5 Figure 10 . 6 :

 1055106 Figure 10.5: MEG Phantom (Elekta): (left) L2 distance between the predicted and actual component (middle) Mean error (in mm) between predicted and actual dipoles localization (right) Fitting time (in seconds)

Figure 10 . 7 :Figure 10 . 8 :

 107108 Figure 10.7: Reconstructing the BOLD signal of missing subjects. (left) Mean R 2 score between reconstructed data and true data. (right) Fitting time.

c

  and mix the data: xc = (W task ) † stask c . The Covariance method generates a new sample of synthetic data in class c by sampling from a Multivariate Gaussian with mean µ c and covariance Σ, where µ c is the class mean and Σ is the covariance of centered task data estimated using Ledoit-Wolf method. In brief, it assumes normality of the data per class.

Figure 12 . 1 :

 121 Figure 12.1: Accuracy of models for eight multi-contrast datasets. Cross validated accuracy of two linear (LDA and LogReg) and one non-linear classifier (MLP) with or without using data augmentation. The improvement yielded by data augmentation is displayed in red. Black error bars indicate standard deviation across splits while white error bars indicate standard deviation across splits with no augmentation.

Figure 12 . 2 :

 122 Figure 12.2: Accuracy of augmented discriminative models when varying k. We use 100 train subjects from the HCP task dataset to train Conditional ICA with k components and generate 200 fake subjects. Classifiers are trained on the train and fake subjects and tested on the left-out 687 subjects. We repeat the procedure for various values of k using 5 random splits per value and report the mean accuracy across splits as a function of k. The dotted line represents the number of components that has been used in our experiments (k = 900).

  published work Thomas Bazeille et al. "Local optimal transport for functional brain template estimation." In: International Conference on Information Processing in Medical Imaging. Springer. 2019, pp. 237-248 13.2.4 Software

  published work Ronan Perry et al. "mvlearn: Multiview Machine Learning in Python." In: Journal of Machine Learning Research 22.109 (2021), pp. 1-7. url: http://jmlr.org/papers/v22/20-1370.html

  published work Manoj Kumar et al. "BrainIAK: The brain imaging analysis kit." In: Aperture (2020). url: https://osf.io/db2ev/

1 a. 1 . 1

 111 Part VII A P P E N D I C E S A M U LT I V I E W I C A a.1 proofs of section 7.Proof of Prop. 12

153 E 2 λ 2 - 1 ,

 153221 , G i is diagonal and constant across subjects (becauseE[f (Λ(s + 1 m j n j ))(n i ) ] = [f (Λ(s + 1 m j n j ))(n i ) ]). Let us therefore consider only its coefficient (a, a), and let λ = Λ aa :G i aa = G(λ) = φ(λ)λ + σ 2 (1 -1/m) σ where φ(λ) = 1 m E[f (λ(s a + 1 m j n j a ))(s a + n i a )

  , the second term gives at the first order:

  1 ab , E 1 ba , . . . , E m ab , E m ba ). Denote ε = Γ E ab , ε = Γ E ba , β = Γ Sab and β = Γ S ba . The linear operator for the block is:K(ε, ε , β, β ) =

Figure A. 1 :

 1 Figure A.1: Reproducing the time-segment matching experiment of [37] [157] Mean classification accuracy -error bars represent 95% confidence interval

  Figure A.1.

Figure A. 2 :

 2 Figure A.2: Synthetic experiment with model x i = A i s + n i

Lemma 23 .

 23 Let us consider the following eigenvalue problem: m ∈ R p,p are positive diagonal matrices and I is the identity matrix. If the first p eigenvalues are distincts, the first p eigenvectors z 1 , . . . , z p , z i ∈ R mp have different first non-zero coordinates.

  K l = diag(Σ 1l . . . Σ ml ) and u = us drop the index l in the notations for simplicity. The problem can be rewritten (uu + K)w = (I + K)wλ (B.6) ⇐⇒ (I + K) -1 (uu + K)w = wλ (B.7)

b. 3

 3 em e-step and m-step for shica with gaussian components b.[START_REF] Ablin | Faster independent component analysis by preconditioning with Hessian approximations[END_REF].1 E-stepThe derivations are the same as in section B.4.1 but the sum over α ∈ 1 2 , 3 2 is replaced by just α = 1.b.3.2 M-stepThe function to minimize in the M-step is then given by:J = E[log p(x, s)] y i -E[s|x])(y i -E[s|x]) ] + V[s|x] ) (B.16) + c (B.17)where c does not depend on Σ i Therefore we get closed-form updates for Σ i :Σ i ← diag(E[(y i -E[s|x])(y i -E[s|x]) ] + V[s|x]) (B.18) Plugging in the closed-form formula for E[s|x] and V[s|x] we get updates that only depends on the covariances Ĉij = E[x i x j ]. Σ i ← diag( Ĉii -2V[s|x] e-step and m-step for shica with non-gaussian components b.4.1 E-stepThe complete likelihood is given byp(x, s) = i p(x i |s)p(s) = ( i Σ -1 ij ) -1. Hiding variable that do not depend on s we obtain i p(x i |s) ∝ j N( ȳj ; s j , Σj ) j ; ȳj , Σj )N(s j ; 0, α) (B.36) Let us now prove a second Lemma: Lemma 26. N(x; y, ν)N(x, 0, α) = N(y; 0, ν + α)N(x; αy α + ν , να α + ν ) Proof. We have N(x; y, ν)N(x, 0, α) = exp -(x-y)

) α∈{0. 5 , 1 .-

 51 5} θ α (B.[START_REF] Barry Erhardt | Comparison of multi-subject ICA methods for analysis of fMRI data[END_REF] where θ α = N( ȳj ; 0, Σj + α).So we obtain the desired result:E[s j |x] = α∈{0.5,1.5} θ α α ȳj α+ Σj α∈{0.5,1.5} θ α (B.51) V[s j |x] = α∈{0.5,1.5} θ α Σj α α+ Σj α∈{0.5,1.5} θ α (B.52) b.4.2 M-stepThe function to minimize in the M-step is then given by:J = E[log p(x, s)] log(|W i |) + log(|Σ i |)+ (B.54) y i -E[s|x])(y i -E[s|x]) ] + V[s|x] ) + c (B.55)

-

  log(|W i |)log(|I k + ε|) -E[log(N(y i + εy i ; s; Σ i ))] + const (B.58) = J(W i )tr(ε) εy i , (Σ i ) -1 (y i -s) ]+ E[ (y i -s), (Σ i ) -1 εy i ] + E[ εy i , (Σ i ) -1 εy i ] + o( ε 2 ) ε a,b E[(Σ i ) -1 (y i -s)(y i ) ] ε a,b E[(Σ i ) -1 εy i (y i ) ] b E[(Σ i ) -1 (y i -s)(y i ) ] b (Σ i ) -1 a,a ε a,d E[y i (y i ) ] a,b = -δ a,b + E[(Σ i ) -1 (y i -s)(y i ) ] ,c,d = δ a,d δ b,c + δ a,c E[y ib y id ] Σ ia (B.63)We approximate the Hessian byH W i a,b,c,d = δ ad δ bc + δ ac δ bd E[(y ib ) 2 ] Σ ia (B.64)

  approaches a Gaussian density as n gets large. This property is called asymptotic normality.

68) < ε (2.69) Where equation (2.68) follows from Cauchy-Schwarz, for the scalar product a, b → E[ab ]. While Proposition 2 states that θn P ----→ n→∞ θ * , Proposition 4 goes further in the analysis. It states that √ n θnθ *

  5 mm. The raiders-full

	Dataset	Subjects Runs	Average run	Voxels
				length	(per subject)
				(in timeframes)
		m			v
	CLIPS	10	17	325	212445
	SHERLOCK	16	5	395	212445
	RAIDERS	10	9	330	212445
	FORREST	19	7	465	212445
	CamCAN	647	1	193	212445

Table 3 .

 3

1: Datasets description

Table 12 . 1 :

 121 Datasets used in the experiments. The table provides references to the datasets that were used for our experiments, with the number of subjects, the number of classes, the number of subjects in train and test set in each cross validation split and the collection number in Neurovault but with some additive non-isotropic Gaussian noise. As in ICA, we estimate W task and S task from X task via ICA. Then we consider R task = X task -W task S task and estimate the covariance Σ R of R

	123

task via LedoitWolf's method. We then generate a data sample xc from class c as with ICA and add Gaussian noise ñ ∼ N(0, Σ R ). Samples are thus generated as xc + ñ.

Table 12 . 2 :

 122 Optimizers and hyper-parameters of classifiersFor each classifier, we give the optimization method used as well as the value of hyperparameters.

	Models	LDA LogReg MLP RF
	Original	0.893 0.874 0.779 0.782
	ICA	0.814 0.840 0.803 0.778
	Covariance	0.895 0.876 0.819 0.780
	ICA + Covariance 0.816 0.840 0.815 0.780
	CGANs	0.874 0.874 0.726 0.779
	Conditional ICA 0.901 0.890 0.832 0.783

Table 12 .3: Comparing augmentation methods based on classification ac- curacy on task HCP dataset We

 12 compare augmentation methods based on the classification accuracy (Acc) obtained by 2 linear classifiers (LDA and LogReg) and two non-linear classifier (MLP and RF) trained on augmented datasets on HCP Task fMRI data. We report the mean accuracy across 5 splits.

Table 12 . 4 :

 124 table 12.4, we give the running-time of the CGAN method and Conditional ICA. This shows that in contrast to deep learning based Running time. We display the running time of conditional ICA and conditional GAN (CGANs) methods used to generate synthetic task fMRI data. Conditional ICA is several orders of magnitude faster than CGANs.

	Methods	Running-time (secs)
	CGANs	11015.1 (≈ 3,05 hr)
	Conditional ICA	62 s

  Table A.2, Table A.3 and Table A.4. All methods are compared upon extraction of components with the same dimensionality (20 components). Table A.3: Reconstructing the BOLD signal of missing subjects: Summary of our quantitative results. We report the mean R2 score across cross-validation splits.

	shica
	Dataset	Method	R2 score Confidence interval
	clips	Chance	0.000	[0.000 ,0.000]
		CanICA	0.110	[ 0.097 , 0.123]
		PCA + ConcatICA	0.075	[ 0.058 , 0.092]
		ConcatICA	0.077	[ 0.059 , 0.094]
		PermICA	0.099	[ 0.087 , 0.111]
		SRM	0.081	[ 0.069 , 0.094]
		MultiViewICA	0.114	[ 0.099 , 0.128]
	forrest	Chance	0.000	[0.000 ,0.000]
		CanICA	0.181	[ 0.169 , 0.193]
		PCA + ConcatICA	0.072	[ 0.054 , 0.090]
		ConcatICA	0.081	[ 0.062 , 0.099]
		PermICA	0.098	[ 0.090 , 0.106]
		SRM	0.180	[ 0.168 , 0.193]
		MultiViewICA	0.191	[ 0.177 , 0.204]
	raiders	Chance	0.000	[0.000 ,0.000]
		CanICA	0.136	[ 0.122 , 0.149]
		PCA + ConcatICA	0.063	[ 0.045 , 0.080]
		ConcatICA	0.062	[ 0.043 , 0.081]
		PermICA	0.107	[ 0.091 , 0.124]
		SRM	0.138	[ 0.121 , 0.154]
		MultiViewICA	0.144	[ 0.124 , 0.164]
	sherlock	Chance	0.000	[0.000 ,0.000]
		CanICA	0.156	[ 0.141 , 0.172]
		PCA + ConcatICA	0.087	[ 0.065 , 0.108]
		ConcatICA	0.091	[ 0.070 , 0.112]
		PermICA	0.067	[ 0.055 , 0.078]
		SRM	0.164	[ 0.147 , 0.181]
		MultiViewICA	0.161	[ 0.142 , 0.180]

C O N C L U S I O N

13.1 a note about resources used All the code is written in Python. We use Matplotlib for plotting [START_REF] Hunter | Matplotlib: A 2D graphics environment[END_REF] , scikit-learn for machine-learning pipelines [START_REF] Pedregosa | Scikitlearn: Machine learning in Python[END_REF], MNE for MEG processing [START_REF] Gramfort | MEG and EEG data analysis with MNE-Python[END_REF], Nilearn for fMRI processing and for its CanICA implementation [START_REF] Abraham | Machine learning for neuroimaging with scikit-learn[END_REF], Brainiak [START_REF] Kumar | BrainIAK tutorials: Userfriendly learning materials for advanced fMRI analysis[END_REF] for its SRM implementation.

13.2 contributions outside of the scope of the thesis Some of the contributions we have made during the thesis extended beyond the scope of multivariate decompositions. We now present these contributions succinctly.

A deep approach to model complex stimuli

In this work, we learn a model to predict fMRI data of subjects watching a movie from the activities of a deep neural network exposed to the same movie. The neural network is previously trained to perform action recognition on a large corpus of movies. The association of activity in visual areas with the different layers of the deep architecture displays complexity-related contrasts across visual areas and reveals a striking foveal/peripheral dichotomy.

The mixing is constrained to be the same across all subjects.

An additional related model is described in [START_REF] Gresele | The Incomplete Rosetta Stone problem: Identifiability results for Multi-view Nonlinear ICA[END_REF]. Similarly to our work, the ICA model has noise on the components side. However, the model involves nonlinear mixings, which are computationally unfeasible to optimize via maximum likelihood; a contrastive learning scheme is therefore adopted, and the likelihood is not derived in closed form.

No evaluation on neuroimaging datasets is presented. Proof. From [START_REF] Comon | Independent component analysis, a new concept?[END_REF], Theorem 10: Assume s = Ps , if the column j of P has more than one non-zero element then s j is Gaussian. Let us define permutations Π 1 , Π 1 such that the first g components of Π 1 s and Π 1 s are Gaussian and P 1 = Π 1 P(Π 1 ) -1 . We can see that P 1 is orthogonal.

We have Π 1 s = P 1 Π 1 s . So the last p -g columns of P 1 contain at most one non-zero element. Using orthogonality of P 1 this non-zero element has value 1 or -1 and is also the only one in its line. Let us focus on column l > g. Assume column l has its non-zero element at index k g. Then line k in P 1 is only non-zero at index l and therefore (Π 1 s) k (which is Gaussian) is equal to (Π 1 s ) l (which is not). Therefore column l can only have its non-zero element at an index greater than g. This shows that P 1 is block diagonal P 1 = O g 0 0 P 2 where O g is orthogonal and P 2 is a sign and permutation matrix. Lemma 22. Assume that Assumption 2 holds for Σ i , and that there is an orthogonal matrix P and diagonal matrices Σ i such that for all i, Σ i = PΣ i P . Then, P is a permutation matrix.

where Σ ij is the coefficient j, j of Σ i and y i = Wx i .

Let us introduce a first lemma: We use m = 496 subjects and fit ShICA-ML with p = 10 components. We localize the components of each subject using sLORETA [START_REF] Domingo | Standardized lowresolution brain electromagnetic tomography (sLORETA): technical details[END_REF].

Then components are registered to a common brain and averaged. Thresholded maps are displayed below along with the time courses of each component. Spatial maps obtained with ShICA-ML highlight the ventral visual cortex and auditive cortex. The results suggest that the response of the auditive cortex is faster and lasts a shorter time than the response of the ventral visual cortex.