Multi-scale Mining of fMRI data with Hierarchical Structured Sparsity

Introduction

Functional magnetic resonance imaging (or fMRI) is a widely used functional neuroimaging modality. Modeling and statistical analysis of fMRI data are commonly done through a linear model, called general linear model (GLM) in the community, that incorporates information about the different experimental conditions and the dynamics of the hemodynamic response in the design matrix. The experimental paradigm consists of a sequence of stimuli, e.g., visual and auditory stimuli, which are included as regressors in the design matrix after convolution with a suitable hemodynamic filter. The resulting model parameters—one coefficient per voxel and regressor—are known as activation maps. They represent the local influence of the different experimental conditions on fMRI signals at the level of individual voxels. The most commonly used approach to analyze these activation maps is called classical inference. It relies on mass-univariate statistical tests (one for each voxel), and yields so-called statistical parametric maps (SPMs) . Such maps are useful for functional brain mapping, but classical inference has some limitations: it suffers from multiple comparisons issues and it is oblivious of the multivariate structure of fMRI data. Such data exhibit natural correlations between neighboring voxels forming clusters with different sizes and shapes, and also between distant but functionally connected brain regions.

To address these limitations, an approach called reverse inference (or “brain-reading”) was recently proposed. Reverse inference relies on pattern recognition tools and statistical learning methods to explore fMRI data. Based on a set of activation maps, reverse inference estimates a function that can then be used to predict a target (typically, a variable representing a perceptual, cognitive or behavioral parameter) for a new set of images. The challenge is to capture the correlation structure present in the data in order to improve the accuracy of the fit, which is measured through the resulting prediction accuracy. Many standard statistical learning approaches have been used to construct prediction functions, among them kernel machines (SVM, RVM) or discriminant analysis (LDA, QDA) . For the application considered in this paper, earlier performance results indicate that we can restrict ourselves to mappings that are linear functions of the data.

Learning the parameters $({\mathbf{w}},b)$ remains challenging since the number of features ($10^{4}$ to $10^{5}$ voxels) exceeds by far the number of samples (a few hundreds of volumes). The prediction function is therefore prone to overfitting in which the learning set is predicted precisely whereas the algorithm provides very inaccurate predictions on new samples (the test set). To address this issue, dimensionality reduction attempts to find a low dimensional subspace that concentrates as much of the predictive power of the original set as possible for the problem at hand.

Feature selection is a natural approach to perform dimensionality reduction in fMRI, since reducing the number of voxels makes it easier to identify a predictive region of the brain. This corresponds to discarding some columns of ${\mathbf{X}}$. This feature selection can be univariate, e.g., analysis of variance (ANOVA) , or multivariate. While univariate methods ignore joint information between features, multivariate approaches are more adapted to reverse inference since they extract predictive patterns from the data as a whole. However, due to the huge number of possible patterns, these approaches suffer from combinatorial explosion, and some costly suboptimal heuristics (e.g., recursive feature elimination ) can be used. That is why ANOVA is usually preferred in fMRI. Alternatively, two more adapted solutions have been proposed: regularization and feature agglomeration.

Regularization is a way to encode a priori knowledge about the weight vector ${\mathbf{w}}$. Possible regularizers can promote for example spatial smoothness or sparsity which is a natural assumption for fMRI data. Indeed, only a few brain regions are assumed to be significantly activated during a cognitive task. Previous contributions on fMRI-based reverse inference include . They can be presented through the following minimization problem:

For the squared loss, when setting $\lambda_{1}$ to 0, the model is called ridge regression, while when $\lambda_{2}=0$ it is known as Lasso or basis pursuit denoising (BPDN) . The essential shortcoming of the Elastic net is that it does not take into account the spatial structure of the data, which is crucial in this context . Indeed, due to the intrinsic smoothing of the complex metabolic pathway underlying the difference of blood oxygenation measured with fMRI , statistical learning approaches should be informed by the 3D grid structure of the data.

In order to achieve dimensionality reduction, while taking into account the spatial structure of the data, one can resort to feature agglomeration. Precisely, new features called parcels are naturally generated via averaging of groups of neighboring voxels exhibiting similar activations. The advantage of agglomeration is that no information is discarded a priori and that it is reasonable to hope that averaging might reduce noise. Although, this approach has been successfully used in previous work for brain mapping , existing work does typically not consider the supervised information (i.e., the target $y$) while exploring the parcels. A recent approach has been proposed to address this issue, based on a supervised greedy top-down exploration of a tree obtained by hierarchical clustering . This greedy approach has proven to be effective especially for inter-subject analyzes, i.e., when training and evaluation sets are related to different subjects. In this context, methods need to be robust to intrinsic spatial variations that exist across subjects: although a preliminary co-registration to a common space has been performed, some variability remains between subjects, which implies that there is no perfect voxel-to-voxel correspondence between volumes. As a result, the performances of traditional voxel-based methods are strongly affected. Therefore, averaging in the form of parcels is a good way to cope with inter-subject variability. This greedy approach is nonetheless suboptimal, as it explores only a subpart of the whole tree.

Based on these considerations, we propose to integrate the multi-scale spatial structure of the data within the regularization term $\Omega$, while preserving convexity in the optimization. This notably guarantees global optimality and stability of the obtained solutions. To this end, we design a sparsity-inducing penalty that is directly built from the hierarchical structure of the spatial model obtained by Ward’s algorithm using a contiguity-constraint . This kind of penalty has already been successfully applied in several contexts, e.g., in bioinformatics, to exploit the tree structure of gene networks for multi-task regression , in log-linear models for the selection of potential orders , in image processing for wavelet-based denoising , and also for topic models . Other applications have emerged in natural language and audio processing .

We summarize here the contributions of our paper:

We explain how the multi-scale spatial structure of fMRI data can be taken into account in the context of reverse inference through the combination of a spatially constrained hierarchical clustering procedure and a sparse hierarchical regularization.

We provide a convex formulation of the problem and propose an efficient optimization procedure.

We conduct an experimental comparison of several algorithms and formulations on fMRI data and illustrate the ability of the proposed method to localize in space and in scale some brain regions involved in the processing of visual stimuli.

The rest of the paper is organized as follows: we first present the concept of structured sparsity-inducing regularization and then describe the different regression/classification formulations we are interested in. After exposing how we handle the resulting large-scale convex optimization problems thanks to a particular instance of proximal methods—the forward-backward splitting algorithm, we validate our approach both in a synthetic setting and on a real dataset.

Combining agglomerative clustering with sparsity-inducing regularizers

As suggested in the introduction, it is possible to construct a tree-structured hierarchy of new features on top of the original voxels using hierarchical clustering. Moreover, spatial constraints can be enforced in the clustering algorithm so that the underlying voxels corresponding to each of these features form localized spatial patterns on the brain similar to the ones we hope to retrieve . Once these features constructed, instead of selecting features in the tree greedily, we propose to cast the feature selection problem as supervised learning problem of the form (1). One of the qualities of the greedy approach however is that it is only allowed to select potentially more noisy features, corresponding to smaller clusters, after the a priori more stable features associated with ancestral clusters in the hierarchy have been selected. As we will show, it is possible to construct a convex regularizer $\Omega$ that has the same property, i.e. that respects the hierarchy, and prioritizes the selection of features in the same way. Naturally, the regularizer has to be constructed directly from the hierarchical clustering of the voxels.

To this end, we consider hierarchical agglomerative clustering procedures . These begin with every voxel ${\mathbf{x}}^{j}$ representing a singleton cluster $\{j\}$, and at each iteration, a selected pair of clusters—according to a criterion discussed below—is merged into a single cluster. This procedure yields a hierarchy of clusters represented as a binary tree $\mathcal{T}$ (also often called a dendrogram) , where each nonterminal node is associated with the cluster obtained by merging its two children clusters. Moreover, the root of the tree $\mathcal{T}$ is the unique cluster that gathers all the voxels, while the leaves are the clusters consisting of a single voxel. From now on, we refer to each nonterminal node of $\mathcal{T}$ as a parcel, which is the union of its children’s voxels (see Figure 1).

Among different hierarchical agglomerative clustering procedures, we use the variance-minimizing approach of Ward’s algorithm . In short, two clusters are merged if the resulting cluster minimizes the sum of squared differences of the fMRI signal within all clusters (also known as inertia criterion). More formally, at each step of the procedure, we merge the clusters $c_{1}$ and $c_{2}$ that minimize the criterion

where we have introduced the average vector $\langle{\mathbf{X}}\rangle_{c}\triangleq\frac{1}{|c|}\sum_{j\in c}{\mathbf{x}}^{j}$. In order to take into account the spatial information, we also add connectivity constraints in the hierarchical clustering algorithm, so that only neighboring clusters can be merged together. In other words, we try to minimize the criterion $\Delta(c_{1},c_{2})$ only for pairs of clusters which share neighboring voxels (see Algorithm 1). This connectedness constraint is important since the resulting clustering is likely to differ from standard Ward’s hierarchical clustering.

where $A_{k}$ stands for the set of ancestors of a node $k$ in $\mathcal{T}$ (including itself).

2 Hierarchical sparsity-inducing norms

In the perspective of inter-subject validation, the augmented space of variables can be exploited in the following way: Since the information of single voxels may be unreliable, the deeper the node in $\mathcal{T}$, the more variable the corresponding parcel’s intensity is likely to be across subjects. This property suggests that, while looking for sparse solutions of (1), we should preferentially select the variables near the root of $\mathcal{T}$, before trying to access smaller parcels located further down in $\mathcal{T}$.

When these groups form a partition of the space of variables, the resulting penalty has been shown to improve the prediction performance and/or interpretability of the learned models, provided that the block structure is relevant (e.g., see and references therein).

If the groups overlap , richer structures can then be encoded. In particular, we follow who first introduced hierarchical sparsity-inducing penalties. Given a node $j$ of $\mathcal{T}$, we denote by $g_{j}\subseteq\{1,\dots,q\}$ the set of indices that record all the descendants of $j$ in $\mathcal{T}$, including itself. In other words, $g_{j}$ contains the indices of the subtree rooted at $j$; see Figure 1. If we now denote by ${\mathcal{G}}$ the set of all $g_{j},\ j\in\{1,\dots,q\}$, that is, ${\mathcal{G}}\triangleq\{g_{1},\dots,g_{q}\}$, we can define our hierarchical penalty as

The family of norms with the previous property is actually slightly larger and we consider throughout the paper norms $\Omega$ of the form :

Supervised learning framework

In this first setting, we naturally consider the squared loss function, so that problem (1) can be reduced to

Note that in this case, we have omitted the intercept $b$ since we can center the vector ${\mathbf{y}}$ and the columns of ${\mathbf{X}}$ instead. Prediction for a new fMRI volume ${\mathbf{x}}$ is then simply performed by computing the dot product ${\mathbf{x}}^{\top}\!{\mathbf{w}}^{*}$.

2 Classification

A standard way of addressing multi-class classification problems consists in using a multi-logit model, also known as multinomial logistic regression (see, e.g., and references therein). In this case, class-conditional probabilities are modeled for each class by a softmax function, namely, given a fMRI volume ${\mathbf{x}}$, the probability of having the $k$-th class label reads

The parameters $\{{\mathbf{W}},{\mathbf{b}}\}$ are then learned by maximizing the resulting (conditional) log-likelihood, which leads to the following optimization problem:

Whereas the regularization term is separable with respect to the different weight vectors ${\mathbf{w}}^{k}$, the loss function induces a coupling in the columns of ${\mathbf{W}}$. As a result, the optimization has to be carried out over the entire matrix ${\mathbf{W}}$. In this setting, and given a new fMRI volume ${\mathbf{x}}$, we make predictions by choosing the label that maximizes the class-conditional probabilities (6), that is, $\text{argmax}_{k\in\{1,\dots,c\}}\text{Prob}(y=k|{\mathbf{x}};{\mathbf{W}}^{*},{\mathbf{b}}^{*})$.

In Section 5, we consider another multi-class classification scheme. The “one-versus-all” strategy (OVA) consists in training $c$ different (real-valued) binary classifiers, each one being trained to distinguish the examples in a single class from the observations in all remaining classes. In order to classify a new example, among the $c$ classifiers, the one which outputs the largest (most positive) value is chosen. In this framework, we consider binary classifiers built from both the squared and the logistic loss functions. If we denote by ${\bar{\mathbf{Y}}}\in\{-1,1\}^{n\times c}$ the indicator response matrix defined as ${\bar{\mathbf{Y}}}_{i}^{k}\triangleq 1$ if ${\mathbf{y}}_{i}=k$ and $-1$ otherwise, we obtain

By invoking the same arguments as in Section 3.1, the vector of intercepts ${\mathbf{b}}$ is again omitted in the above problem with the squared loss. Moreover, given a new fMRI volume ${\mathbf{x}}$, we predict the label $k$ that maximizes the response ${\mathbf{x}}^{\top}\![{\mathbf{w}}^{*}]^{k}$ among the $c$ different classifiers. The case of the logistic loss function parallels the setting of the multinomial logistic regression, where each of the $c$ “one-versus-all” classifiers leads to a class-conditional probability; the predicted label is the one corresponding to the highest probability.

The formulations we have reviewed in this section can be solved efficiently within the same optimization framework that we now introduce.

Optimization

The convex minimization problem (1) is challenging, since the penalty $\Omega$ as defined in (5) is non-smooth and the number of variables to consider is large (we have $q\approx 10^{5}$ variables in the following experiments). These difficulties are well addressed by forward-backward splitting methods, which belong to the broader class of proximal methods. Forward-backward splitting schemes date back (at least) to and have been further analyzed in various settings (e.g., see ); for a thorough review of proximal splitting techniques, we refer the interested readers to .

Our convex minimization problem (1) can be handled well by such techniques since it is the sum of two semi-lower continuous, proper, convex functions with non-empty domain, and where one element—the loss function $\mathcal{L}({\mathbf{y}},{\mathbf{X}},.)$—is assumed differentiable with Lipschitz-continuous gradient (which notably covers the cases of the squared and simple/multinomial logistic functions, as introduced in Section 3).

This operator was initially introduced by Moreau to generalize the projection operator onto a convex set; for a complete study of the properties of $\text{Prox}_{\lambda\Omega}$, see . Based on definition (7), and given the current iterate ${\mathbf{w}}^{(k)}$,For clarity of the presentation, we do not consider the optimization of the intercept that we leave unregularized in all our experiments. the typical update rule of forward-backward splitting methods has the formFor simplicity, we only present a constant-stepsize scheme; adaptive line search can also be used in this context and can lead to larger stepsizes .

where $L>0$ is a parameter which is a upper bound on the Lipschitz constant of the gradient of $\mathcal{L}$. In the light of the update rule (8), we can see that solving efficiently problem (7) is crucial to enjoy good performance. In addition, when the non-smooth term $\Omega$ is not present, the previous proximal problem (8), also known as the implicit or backward step, exactly leads to the standard gradient update rule.

In the subsequent experiments, we focus on accelerated multi-step versions of forward-backward splitting methods (see, e.g., ),More precisely, we use the accelerated proximal gradient scheme (FISTA) taken from . The Matlab/C++ implementation we use is available at http://www.di.ens.fr/willow/SPAMS/. where the proximal problem (8) is not solved for a current estimate, but for an auxiliary sequence of points that are linear combinations of past estimates. These accelerated versions have increasingly drawn the attention of a broad research community since they can deal with large non-smooth convex problems, and their convergence rates on the objective achieve the complexity bound of $O(1/k^{2})$, with $k$ denoting the iteration number. As a side comment, note that as opposed to standard one-step forward-backward splitting methods, nothing can be said about the convergence of the sequence of iterates themselves. In our case, the cost of each iteration is dominated by the computation of the gradient (e.g., $O(np)$ for the squared loss) and the proximal operator, whose time complexity is linear, or close to linear, in $p$ for the tree-structured regularization .

Experiments and results

We now present experimental results on simulated data and real fMRI data.

The choice of the weights and of the correlation introduced in images aim at illustrating how the hierarchical regularization estimates weights at different resolutions in the image. The targets were simulated by forming ${\mathbf{w}}^{\top}{\mathbf{x}}^{(i)}$ corrupted with an additive white noise (SNR=10dB). The loss used was the squared loss as detailed in Section 3.1. The regularization parameter was estimated with two-fold cross-validation (150 images per fold) on a logarithmic grid of 30 values between $10^{3}$ and $10^{-3}$.

The components of the estimated weight vector ${\mathbf{w}}^{*}$ at different scales are presented in the images of Fig. 2 , with each image corresponding to a different depth in the tree. For a given tree depth, an image is formed from the corresponding parcellation. All the voxels within a parcel are colored according to the associated scalar in ${\mathbf{w}}^{*}$. It can be observed that all three patterns are present in the weight vector but at different depths in the tree. The small activation in the top-right hand corner shows up mainly at scale 3 while the bigger patterns appear higher in the tree at scales 5 and 6. This simulation clearly illustrates the ability of the method to capture informative spatial patterns at different scales.

2 Description of the fMRI data

We apply the different methods to analyze the data of ten subjects from an fMRI study originally designed to investigate object coding in visual cortex (see for details). During the experiment, ten healthy volunteers viewed objects of two categories (each one of the two categories is used in half of the subjects) with four different exemplars in each category. Each exemplar was presented at three different sizes (yielding $12$ different experimental conditions per subject). Each stimulus was presented four times in each of the six sessions. We averaged data from the four repetitions, resulting in a total of $n=72$ volumes by subject (one volume of each stimulus by session). Functional volumes were acquired on a 3-T MR system with eight-channel head coil (Siemens Trio, Erlangen, Germany) as T2*-weighted echo-planar image (EPI) volumes. Twenty transverse slices were obtained with a repetition time of 2s (echo time, 30ms; flip angle, $70^{\circ}$; $2\times 2\times 2$-mm voxels; $0.5$-mm gap). Realignment, spatial normalization to MNI space, slice timing correction and GLM fit were performed with the SPM5 softwarehttp://www.fil.ion.ucl.ac.uk/spm/software/spm5.. In the GLM, the time course of each of the $12$ stimuli convolved with a standard hemodynamic response function was modeled separately, while accounting for serial auto-correlation with an AR(1) model and removing low-frequency drift terms with a high-pass filter with a cut-off of 128s ($7.8\times 10^{-3}$Hz). In the present work we used the resulting session-wise parameter estimate volumes. Contrary to a common practice in the field the data were not smoothed with an isotropic Gaussian filter. All the analysis are performed on the whole acquired volume.

The four different exemplars in each of the two categories were pooled, leading to volumes labeled according to the three possible sizes of the object. By doing so, we are interested in finding discriminative information to predict the size of the presented object.

This can be reduced to either a regression problem in which our goal is to predict the class label of the size of the presented object (i.e., $y\in\{0,1,2\}$),An interesting alternative would be to consider some real-valued dimension such as the field of view of the object. or a three-category classification problem, each size corresponding to a category. We perform an inter-subject analysis on the sizes both in regression and classification settings. This analysis relies on subject-specific fixed-effects activations, i.e., for each condition, the six activation maps corresponding to the six sessions are averaged together. This yields a total of 12 volumes per subject, one for each experimental condition. The dimensions of the real data set are $p\approx 7\times 10^{4}$ and $n=120$ (divided into three different sizes). We evaluate the performance of the method by cross-validation with a natural data splitting, leave-one-subject-out. Each fold consists of 12 volumes. The parameter $\lambda$ of all methods is optimized over a grid of 30 values of the form $2^{k}$, with a nested leave-one-subject-out cross-validation on the training set. The exact scaling of the grid varies for each model to account for different $\Omega$.

3 Methods involved in the comparisons

First of all, when the regularization $\Omega$ as defined in (5) is employed, we consider three settings of values for $(\eta_{g})_{g\in{\mathcal{G}}}$ which leverage the tree structure $\mathcal{T}$. More precisely, we set $\eta_{g}=\rho^{\textrm{depth}(g)}$ for $g$ in ${\mathcal{G}}$, with $\rho\in\{0.5,1,1.5\}$ and where $\textrm{depth}(g)$ denotes the depth of the root of the group $g$ in $\mathcal{T}$. In other words, the larger $\rho$, the more averse we are to selecting small (and variable) parcels located near the leaves of $\mathcal{T}$. As the results illustrate it, the choice of $\rho$ can have a significant impact on the performance. More generally, the problem of selecting $\rho$ properly is a difficult question which is still under investigation, both theoretically and practically, e.g., see .

The greedy approach from is included in the comparisons, for both the regression and classification tasks. It relies on a top-down exploration of the tree $\mathcal{T}$. In short, starting from the root parcel that contains all the voxels, we choose at each step the split of the parcel that yields the highest prediction score. The exploration step is performed until a given number of parcels is reached, and yields a set of nested parcellations with increasing complexity. Similarly to a model selection step, we chose the best parcellation among those found in the exploration step. The selected parcellation is thus used on the test set. In the regression setting, this approach is combined with Bayesian ridge regression, while it is associated with a linear support vector machine for the classification task (whose regularization parameter, often referred to as $C$ in the literature , is found by nested cross-validation in $\{0.01,0.1,1\}$).

3.2 Classification setting

4 Results

We present results comparing our approach based on the hierarchical sparsity-inducing norm (5) with the regularization listed in the previous section. For each method, we computed the cross-validated prediction accuracy and the percentage of non-zero coefficients, i.e., the level of sparsity of the model.

In terms of sparsity, we can see, as expected, that ridge regression does not yield any sparsity and that the Lasso solution is very sparse (in the feature space, with approximately $7\times 10^{4}$ voxels). Our method yields a median value of 9.36% of non-zero coefficients (in the augmented space of features, with about $1.4\times 10^{5}$ nodes in the tree). The maps of weights obtained with Lasso and the hierarchical regularization for one fold, are reported in Fig. 3. The Lasso yields a scattered and overly sparse pattern of voxels, that is not easily readable, while our approach extracts a pattern of voxels with a compact structure, that clearly outlines brain regions expected to activate differentially for stimuli with different low-level visual properties, e.g., sizes; the early visual cortex in the occipital lobe at the back of the brain. Interestingly, the patterns of voxels show some symmetry between left and right hemispheres, especially in the primary visual cortex which is located at the back and center of the brain. This observation is consistent with the current understanding in neuroscience that the symmetric parts of this brain region process respectively the visual contents of each of the visual hemifields. The weights obtained at different depth level in the tree, corresponding to different scales, show that the largest coefficients are concentrated at the higher scales (scale 6 in Fig. 3), which suggest that the object sizes cannot be well decoded at the voxel level but require features corresponding to larger clusters of voxels.

4.2 Classification results

Conclusion

In this article, we introduced a hierarchically structured regularization, which takes into account the spatial and multi-scale structure of fMRI data. This approach copes with inter-subject variability in a similar way as feature agglomeration, by averaging neighboring voxels. Although alternative agglomeration strategies do exist, we simply used the criterion which appears as the most natural, Ward’s clustering, and which builds parcels with little variance.

Results on a real dataset show that the proposed algorithm is a promising tool for mining fMRI data. It yields similar or higher prediction accuracy than reference methods, and the map of weights it obtains exhibit a cluster-like structure. It makes them easily readable compared to the overly sparse patterns found by classical sparsity-promoting approaches.

Finally, it should be mentioned that the performance achieved by this approach in inter-subject problems suggests that it could potentially be used successfully for medical diagnosis, in a context where brain images –not necessarily functional images– are used to classify individuals into diseased or control population. Indeed, for difficult problems of that sort, where the reliability of the diagnostic is essential, the stability of models obtained from convex formulations and the interpretability of sparse and localized solutions are quite relevant to provide a credible diagnostic.

Acknowledgments

The authors acknowledge support from the ANR grants ViMAGINE ANR-08-BLAN-0250-02 and ANR 2010-Blan-0126-01 “IRMGroup”. The project was also partially supported by a grant from the European Research Council (SIERRA Project).

References