G Varoquaux

B Thirion

Statistical learning for resting-state fMRI: successes and challenges

In the absence of external stimuli, fluctuations in cerebral activity can be used to reveal intrinsic structures. Well-conditioned probabilistic models of this so-called resting-state activity are needed to support neuroscientific hypotheses. Exploring two specific descriptions of resting-state fMRI, namely spatial analysis and connectivity graphs, we discuss the progress brought by statistical learning techniques, but also the neuroscientific picture that they paint, and possible modeling pitfalls.

Introduction

The study of the fluctuations in on-going brain activity, present at rest, opens a window on intrinsic functional brain structures. Indeed, correlation patterns observed in fMRI without external stimuli mirror, to some extent, task-related networks [START_REF] Biswal | Functional connectivity in the motor cortex of resting human brain using echo-planar MRI[END_REF]. Using the resting brain to accumulate knowledge on its function requires statistical models of the signal, though too often descriptive statistics, such as sample correlations, are used to draw inferences. The estimation of such models is limited, at the subject level, by the scarcity of the data available, and, at the population level, by intersubject variability.

Here we discuss what, from our point of view, statistical machine learning has brought to resting-state modeling, namely well-conditioned models and estimators in high-dimensional settings due to a mature approach to regularization. Learning higher-dimensional models brings richer neuroscientific pictures, but all the modeling issues cannot be addressed with machine learning. Due to its brevity, this paper presents only a small view of the resting-state modeling literature. First we revisit the spatial extraction of resting-state networks and show that adding penalization to models yields finer-grained descriptions. Second we consider inferring connectivity graphs, for which regularization is crucial but can also shape our understanding of brain wiring.

2 Spatial analysis: better models define finer networks 2.1 A handful of large-scale brain networks or a score of regions?

Early insight on resting-state came from seed-analysis studies, such as in the motor cortex [START_REF] Biswal | Functional connectivity in the motor cortex of resting human brain using echo-planar MRI[END_REF], that highlighted brain networks correlating with a given seed region during rest. However, a wide range of different networks can be obtained by varying the choice of the seed [START_REF] Cole | Advances and pitfalls in the analysis and interpretation of resting-state fmri data[END_REF]. For this reason, investigators have turned to multivariate methods that do not require such arbitrary choice, in particular Independent Component Analysis (ICA), to extract the spatial modes of restingstate activity. These approaches have identified a small number of large-scale brain networks consistently found across studies and subjects [START_REF] Damoiseaux | Consistent resting-state networks across healthy subjects[END_REF]. Recent studies argued that these networks could be divided in a set of smaller consistent brain regions [START_REF] Varoquaux | A group model for stable multi-subject ICA on fMRI datasets[END_REF][START_REF] Kiviniemi | Functional segmentation of the brain cortex using high model order group PICA[END_REF]. From a multivariate analysis point of view, the complexity of the models that can be learned is limited by the amount of data at hand. Is the small number of networks an over simplification due to data scarcity? At which scale can one extract brain networks?

The importance of penalizing: ICA to sparse decompositions

While data-driven methods such as ICA are often called model-free, they rely on the following simple multivariate decomposition model of the signal: Y = U V + N, where Y ∈ R n×p are the fMRI data: n observations of images with p voxels, U ∈ R n×k and V ∈ R k×p are respectively k time series and k spatial maps learned from the data, and N ∈ R n×p some additional noise. This model needs additional constraints to be well-posed. ICA imposes low mutual-information between the spatial maps. However, it has no intrinsic noise model: it fits all the available signal. Thus, in practice, ICA is most often combined with an preliminary PCA step to select a smaller signal subspace accounting for most of the original signal variance [START_REF] Beckmann | Probabilistic independent component analysis for functional magnetic resonance imaging[END_REF]. In addition, output maps are thresholded to remove background noise. The resulting maps that are considered as neuroscientifically meaningful are spatially structured and highlight a small number of activated brain regions. ICA also segments artifacts in the fMRI data, such as vascular or movement related signal. However, if the model order, i.e. the number of components, is increased, some output maps show little structure and are difficult to interpret. This limitation stems from the inability of ICA to reject noise: at high model order the PCA step used for this purpose includes signal that cannot be well decomposed according to the ICA criterion. The two step estimation PCA + ICA is necessary as the ICA criterion cannot be combined with a criterion controlling explained variance. In Bayesian terms, it corresponds to an improper prior. A solution to this limitation is to frame the model in terms of sparsity, rather than independence. Indeed ICA selects sparse maps as a byproduct of the independence criterion, and this sparsity is central to its success on fMRI data [START_REF] Daubechies | Independent component analysis for brain fMRI does not select for independence[END_REF][START_REF] Varoquaux | ICA-based sparse features recovery from fMRI datasets[END_REF].

An ℓ 1 penalization can impose a sparse prior on the spatial maps V in the decomposition model. Combined with PCA in a sparse PCA procedure, it regularizes model estimation by constraining the fitted signal [START_REF] Varoquaux | Multi-subject dictionary learning to segment an atlas of brain spontaneous activity[END_REF]. As ℓ 1 penalizations perform poorly with correlated features, it is useful to impose sparsity on groups of neighboring voxels, for instance using sparse structured penalties. A Sparse structured PCA (SSPCA) [START_REF] Jenatton | Structured sparse principal component analysis[END_REF] extracts the regions comprising the brain networks in separate maps (see Fig. 1). At high model order SSPCA remains stable and interpretable as it learns finer-grained segmentations of the brain. Some ICA-extracted maps show little structure and appear as noise. On the opposite, all the SSPCA maps segment regions, but their spatial regularity may be considered as forced. MSDL maps are more appealing: less noisy and less geometric.

Multi-subject analysis, goldmine or minefield?

To increase statistical power, data is often pooled across subjects. Yet intersubject modeling of resting-state is challenging as the spatial maps differ across subjects and the time-series are not related. Most often single-subject procedures are applied to the time-concatenated volumes of the different subjects [START_REF] Calhoun | A method for making group inferences from fMRI data using independent component analysis[END_REF]. Such concatenation models the data in the different subjects as drawn from the same distribution, i.e. sharing the same spatial maps V. However, for ICA-based approaches, separating signal-space selection at the subject and the group level stabilizes estimation in high-model-order settings and extracts more meaningful maps [START_REF] Varoquaux | A group model for stable multi-subject ICA on fMRI datasets[END_REF]. In a sparse decomposition framework, group-level maps can be formulated as a prior on subject-level models, as in the Multi-Subject Dictionary Learning (MSDL) approach [START_REF] Varoquaux | Multi-subject dictionary learning to segment an atlas of brain spontaneous activity[END_REF]. The corresponding joint estimation of subjectlevel and group-level maps yields a reference population atlas of brain networks and corresponding subject-specific counter-parts. While multi-subject analyses do stabilize model estimation (see Fig. 1), at very high model order the spatial maps may reflect modes of inter-subject variability, such as spatial gradients of some of the networks [START_REF] Varoquaux | Multi-subject dictionary learning to segment an atlas of brain spontaneous activity[END_REF].

3 Connectivity graphs: methods shape our understanding

Disentangling indirect effects to infer functional connections

Given brain regions, an important problem is to quantify their functional interactions from fRMI observations. Traditionally, an opposition is made between functional connectivity models, that describe the correlations in the time-series, and effective connectivity models, reflecting the influence neural systems exerts over each other. In our opinion, this dichotomy is over-emphasized and there is a continuum of models ranging from descriptive statistics on correlation matrices to inference of directed graphical models. Unless the number of regions is very small, graphical model estimation is ill-posed in the small-sample limit as the likelihood of a connection between a pair of regions is a function of the signal observed in these regions conditionally on all the other regions. Indeed, isolating the signal specific to that pair is necessary to disambiguate direct links from indirect effects. The simplest probabilistic model for fMRI is the multivariate normal, for which estimating graphical models boils to partial-correlation estimation. As with effective connectivity, these graphical models depict coupling between brain regions, yet as they involve only second-order statistics, they are close to correlations analysis usually classified as functional connectivity [START_REF] Mcintosh | Moving between functional and effective connectivity[END_REF].

Regularizing inverse covariance estimation

As partial correlations are given by the inverse covariance matrix, the sparsity pattern of this matrix captures the graph of conditional independence between brain regions and it is a good candidate for brain connectivity graphs [START_REF] Smith | Network modelling methods for fMRI[END_REF]. If the number of observations is not very large compared to the number of regions, the estimation of the inverse covariance must be regularized. Different regularization schemes yield different priors for connectivity structures [START_REF] Varoquaux | Markov models for fMRI correlation structure: is brain functional connectivity small world, or decomposable into networks[END_REF]. To fit to the data, ℓ 1 penalized estimators create overly dense graphs. Greedy methods based on the PC algorithm, that proceed by pruning the graph, fail in the presence of high-degree nodes [START_REF] Varoquaux | Markov models for fMRI correlation structure: is brain functional connectivity small world, or decomposable into networks[END_REF]. Finally top-down approaches decomposing the graph in sub-graphs -known as decomposable models in Gaussian graphs-are limited by the absence of small separable sub-systems in the brain [START_REF] Varoquaux | Markov models for fMRI correlation structure: is brain functional connectivity small world, or decomposable into networks[END_REF].

In our opinion, a promising approach for extracting brain connectivity relies on imposing a common graph structure across subjects to estimate jointly multiple subject-specific models [START_REF] Varoquaux | Brain covariance selection: better individual functional connectivity models using population prior[END_REF]. A penalty based on mixed-norms can set such population-level sparsity in a well-behaved -convex-estimation problem. The resulting population prior accumulates samples for better structure estimation.

Focusing on graph modularity rather than small-worldness

A hallmark of many current brain connectivity studies is their focus on smallword properties of brain graphs: efficient transport with very few connections [START_REF] Bullmore | Complex brain networks: graph theoretical analysis of structural and functional systems[END_REF]. However, in empirically-defined networks, these properties are very fragile to noise [START_REF] Bialonski | From brain to earth and climate systems: Small-world interaction networks or not?[END_REF][START_REF] Varoquaux | Markov models for fMRI correlation structure: is brain functional connectivity small world, or decomposable into networks[END_REF]. For instance a bad control on false detections of connections, as with non-regularized inverse covariance estimates, will create spurious random links in essence fabricating a small-world network. On the other hand, small-world connectivity is hard to recover [START_REF] Bento | Which graphical models are difficult to learn?[END_REF] and estimation procedures fail on high-degree nodes or graphs with many large loops. Given these pitfalls, what evidence on brain connectivity can we gather from fMRI?

The failure of decomposable models and greedy approaches suggests that a good representation of brain connectivity displays a large graph bandwidth and contains hub nodes. The brain networks identified by spatial analysis suggest that brain connectivity contains large structures. An appealing picture is that of a modular graph composed of communities with inter-communities connections large in comparison with extra-communities connections [START_REF] He | Uncovering intrinsic modular organization of spontaneous brain activity in humans[END_REF][START_REF] Bullmore | Complex brain networks: graph theoretical analysis of structural and functional systems[END_REF]. Identifying these communities is a clustering task. It is general a non-convex problem with little guarantees of success from the learning theory standpoint. It can be partly relaxed in a better-behaved graph-cut problem [START_REF] Varoquaux | Brain covariance selection: better individual functional connectivity models using population prior[END_REF]. The resulting functionalsegregation structure highlights from the connectivity graphs large clique-like sub-systems that match well the known brain networks (Fig. 2) [START_REF] Varoquaux | Brain covariance selection: better individual functional connectivity models using population prior[END_REF]. The fullbrain connectivity graph in which they are embedded also reveals their functional interactions (Fig. 2).

Conclusion: with great priors comes great responsibility

Probabilistic modeling of resting-state is a hard problem because it entails fitting very-high-dimensional models in unsupervised settings. Going beyond simple pictures requires moving away from descriptive statistics and applying state-of-theart statistical learning techniques such as adequate penalizations on the various estimation steps. However penalizations are nothing but well-posed formulations of priors, and they can shape the output of the methods. These priors must not only be consistent with current neuroscientific knowledge but also robust to the confounding structured noise of fMRI. Performing principled model selection to pick the best performing prior is still an open question, in the absence of clearly defined statistical task. Indeed, simply selecting for goodness of fit can inflate artificially model complexity to explain structured noise or intersubject variability [START_REF] Varoquaux | Multi-subject dictionary learning to segment an atlas of brain spontaneous activity[END_REF]. For instance, it fails to capture an optimal scale for the description. In addition, while model fit on test data may guide the choice of a penalty, it is intrinsic to the model postulated.

Over-interpretation of plausible results is a recurrent temptation. Given the shortage of observations compared to the complexity of the models, how much can we conclude from the data, and how much is embedded in our priors? We believe that existing large multi-subject and multi-centric databases [START_REF] Biswal | Toward discovery science of human brain function[END_REF] can help to overcome these questions as they call for less regularization and can establish and calibrate priors for smaller studies.

Fig. 1 .

 1 Fig.1. Output of different decomposition models: from left to right, ICA, SSPCA and MSDL. top: default-mode network; for ICA maps are thresholded but for SSPCA and MSDL the transparency indicates zeros resulting from sparsity. bottom: outline of all maps. Some ICA-extracted maps show little structure and appear as noise. On the opposite, all the SSPCA maps segment regions, but their spatial regularity may be considered as forced. MSDL maps are more appealing: less noisy and less geometric.

Fig. 2 .

 2 Fig. 2. (a) Functional connectivity learned with population-sparse inverse covariance; (b) corresponding graph communities; (c) inter-communities interactions [16].


Vincent Michel

Alexandre Gramfort

Evelyn Eger

Gaël Varoquaux

Bertrand Thirion

A comparative study of algorithms for intra-and inter-subjects fMRI decoding

provides an unique opportunity to study brain functional architecture, while being minimally invasive. Reverse inference, a.k.a. decoding, is a recent statistical analysis approach that has been used with success for deciphering activity patterns that are thought to fit the neuroscientific concept of population coding. Decoding relies on the selection of brain regions in which the observed activity is predictive of certain cognitive tasks. The accuracy of such a procedure is quantified by the prediction of the behavioral variable of interest -the target. In this paper, we discuss the optimality of decoding methods in two different settings, namely intraand inter-subject kind of decoding. While inter-subject prediction aims at finding predictive regions that are stable across subjects, it is plagued by the additional inter-subject variability (lack of voxel-to-voxel correspondence), so that the best suited prediction algorithms used in reverse inference may not be the same in both cases. We benchmark different prediction algorithms in both intra-and inter-subjects analysis, and we show that using spatial regularization improves reverse inference in the challenging context of inter-subject prediction. Moreover, we also study the different maps of weights, and show that methods with similar accuracy may yield maps with very different spatial layout of the predictive regions.

Introduction

Reverse inference [START_REF] Dehaene | Inferring behavior from functional brain images[END_REF][START_REF] Cox | Functional magnetic resonance imaging (fMRI) "brain reading": detecting and classifying distributed patterns of fMRI activity in human visual cortex[END_REF], a.k.a. decoding, is an approach for mining fMRI data that uses pattern analysis in order to reveal the information produced by brain activations. The core of this approach is to consider fMRI data analysis as a pattern recognition problem, i.e. using a pattern of voxels to predict a behavioral, perceptual or cognitive variable. In such studies, the accuracy of the prediction can be used to assess whether the pattern of voxels used in the predictive model actually encodes the information about the variable of interest. This approach has been used more frequently in intra-subject settings than in inter-subject analysis. The main interest of inter-subject prediction is to find predictive regions that are stable across subjects, and thus obtain a population-level validation of cognitive hypothesis. The major bottleneck in inter-subject predictions is that such studies are plagued by the inter-subject variability (lack of voxel-to-voxel correspondence) [START_REF] Tucholka | Prise en compte de l'anatomie cérébrale individuelle dans les études d'IRM fonctionnelle[END_REF][START_REF] Tahmasebi | Quantification of Inter-subject Variability in Human Brain and Its Impact on Analysis of fMRI Data[END_REF]. Functional activity localization can vary across subjects due to differences in anatomical structure and in functional organization. As a result, it is challenging to find a common spatial layout of the cognitive substrate across different subjects. In this paper we compare different prediction algorithms in both intra-and inter-subjects settings, in order to investigate the properties required for good inter-subject prediction. We show that using spatial regularization improves the performances in the case of inter-subject studies, by gaining robustness against the spatial variability of the fMRI signal. We also compare the maps obtained by the different methods, and show variability in the spatial support of the predictive regions.

Methods

We briefly introduce the following predictive linear model, in regression settings y = X w + b, where y ∈ R n represents the behavioral variable and (w, b) are the parameters to be estimated on a training set. A vector w ∈ R p can be seen as an image; p is the number of features (or voxels) and b ∈ R is called the intercept. The matrix X ∈ R n×p is the design matrix. Each row is a p-dimensional sample, i.e., an activation map related to the observation. The model performance is evaluated using ζ, the ratio of explained variance (or R 2 coefficient), where

ζ(y t , ŷt ) = var(y t )-var(y t -ŷ t ) var(y t )
. We now detail the different reference methods that will be used in this study.

Non spatially-regularized methods

All the methods are used after an Anova-based feature selection, as this increases their performance. This selection is performed on the training set of each fold in an internal cross-validation loop, and the optimal number of voxels is selected within the range {50, 100, 250, 500}. This feature selection is performed for each method and for each training set separately. It yields different sets of features as the selection is done jointly with the regression within the cross-validation loop. Indeed, some methods such as Elastic Net can perform their own multivariate feature selection, but this step of univariate feature selection allows to reduce the number of features for the regression methods, thus decreasing the computational time.

In this paper, the implementation of Elastic net is based on coordinate descent [START_REF] Friedman | Regularization paths for generalized linear models via coordinate descent[END_REF], while SVR is based on LibSVM [START_REF] Chang | LIBSVM: a library for support vector machines[END_REF]. Methods are used from Python via the Scikit-learn open source package [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF].

Support Vector Regression -SVR -The first prediction function used in reverse inference [START_REF] Cox | Functional magnetic resonance imaging (fMRI) "brain reading": detecting and classifying distributed patterns of fMRI activity in human visual cortex[END_REF] has been Support Vector Machine (SVM) [START_REF] Cortes | Support vector networks[END_REF]. This approach is widely used and has become the reference approach for fMRI reverse inference. Its success comes from its wide availability and good performance on high-dimensional data. In this paper, we use SVR with a linear kernel. The C parameter is optimized by internal cross-validation in the range 10 -3 to 10 1 in multiplicative steps of 10.

Elastic net regularization -Other approaches include built-in feature selection: a generic formulation is given by Elastic net [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF], which uses a combined ℓ 1 (Lasso, parametrized by λ 1 ) and ℓ 2 (Ridge Regression, parametrized by λ 2 ) penalization. While setting many weights to zero, Elastic net, unlike Lasso, can extract more features than samples and correlated features. Elastic net is therefore an attractive approach for reverse inference, as we expect to extract some groups of correlated features, while seeking for an interpretable model (i.e. few selected groups). We use a cross-validation procedure within the training set to optimize the parameters, with λ 1 ∈ {0.2 λ, 0.1 λ, 0.05 λ, 0.01 λ} ( λ = X T y ∞ ), and λ 2 ∈ {0.1, 0.5, 1., 10., 100.}.

Bayesian Ridge Regression -BRR -Bayesian Ridge regression is based on the following Gaussian assumptions p(y|X, w, α)

= n i=1 N (y i |X i w, α -1 ), p(ǫ|α) = N (0, α -1 I n ) (Gaussian noise)
, and p(w|λ) = N (w|0, λ -1 I p ). This Bayesian framework includes the explicit estimation of the α and λ parameters, through the direct maximization of the marginal likelihood. The price to pay is the non-convexity of the whole estimation procedure.

Automatic Relevance Determination -ARD -One can use a more stringent prior on w, known as ARD [START_REF] Mackay | Bayesian interpolation[END_REF][START_REF] Neal | Bayesian Learning for Neural Networks[END_REF], where we assume that the weights w i are drawn from independent Gaussian distributions centered about zero and with a precision

λ i (λ i = λ j if i = j): p(w|λ) = N (0, Λ -1 ) with Λ = diag {λ 1 , ..., λ p } .
This choice of hyper-parameters yields very sparse models.

Multi-Class Sparse Bayesian Regression (MCBR) -MCBR -We also use the MCBR approach [START_REF] Michel | Multiclass Sparse Bayesian Regression for fMRI-Based Prediction[END_REF], which is an intermediate between BRR and ARD. MCBR consists in grouping the features into Q different classes, and regularizing these classes differently; as a consequence, it requires the estimation of fewer parameters than ARD and is far more adaptive than BRR. In this paper, we set K = 9, with weakly informative priors λ 1,k = 10 k-4 , k ∈ [1, .., K] and λ 2,k = 10 -2 , k ∈ [1, .., K], and α 1 = α 2 = 1, following the work in [START_REF] Michel | Multiclass Sparse Bayesian Regression for fMRI-Based Prediction[END_REF]. MCBR can be estimated using Gibbs sampling (Gibbs-MCBR) or Variational Bayes (VB-MCBR).

Spatially-regularized methods

Searchlight The searchlight [START_REF] Kriegeskorte | Information-based functional brain mapping[END_REF] is a brain mapping approach, that yields the amount predictive information conveyed by the voxels in any sub-region about the target variable. This approach is used here for comparing the weights maps obtained by the predictive approaches. In this paper, we use spherical regions with a radius of two voxels, combined with an SVR function (C = 1).

Supervised clustering -SC The supervised clustering algorithm [START_REF] Michel | A supervised clustering approach for fMRI-based inference of brain states[END_REF] is a procedure that creates parcels (i.e. spatially structured group of voxels), while considering the target to be predicted as early as in the clustering procedure. It yields an adaptive segmentation into both large regions and fine-grained information, and can thus be considered as multi-scale. We used the SC with BRR, and we set ∆ = 75 (depth of exploration), following the work in [START_REF] Michel | A supervised clustering approach for fMRI-based inference of brain states[END_REF].

Total Variation regularization -TV TV is defined as the ℓ 1 norm of the image gradient, and has primarily been used for image denoising [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] as it preserves edges. The motivation for using TV for brain imaging [START_REF] Michel | Total variation regularization for fMRI-based prediction of behaviour[END_REF] is that it promotes estimates of the weights with a block structure, therefore outlining the brain regions involved in the target behavioral variable. A particularly important property of this approach is its ability to create spatially coherent regions with similar weights, yielding simplified and informative sets of features. We use TV with a regularization parameter λ = 0.05, following the work in [START_REF] Michel | Total variation regularization for fMRI-based prediction of behaviour[END_REF].

Experiments

We apply the different methods on a real fMRI dataset related to an experiment studying the representation of objects, as detailed in [START_REF] Eger | Graded size sensitivity of object exemplar evoked activity patterns in human loc subregions[END_REF]. During this experiment, ten healthy volunteers viewed objects from one of two categories (each one of the two categories used in equal halves of subjects) with 4 different exemplars each shown in 3 different sizes (yielding 12 different experimental conditions), with 4 repetitions of each stimulus in each of the 6 sessions. We averaged data from the 4 repetitions, resulting in a total of n = 72 images by subject (one image of each stimulus by session). Functional images 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 2 s (echo time, 30 ms; flip angle, 70 • ; 2×2×2-mm voxels; 0.5-mm gap). Realignment, normalization to MNI space, and General Linear Model (GLM) fit were performed with the SPM5 software http://www.fil.ion.ucl.ac.uk/spm/software/spm5. In the GLM, the effect of each of the 12 stimuli convolved with a standard hemodynamic response function was modeled separately, while accounting for serial autocorrelation with an AR(1) model and removing low-frequency drift terms by a high-pass filter with a cutoff of 128 s. In the present work we used the resulting session-wise parameter estimate images. All the analysis are performed on the whole brain volume.

Intra-subject experiment

The four different shapes of objects (from either category) were pulled for each of the three sizes, and we are interested in finding discriminative information between sizes. Each subject is evaluated independently, in a leave-one-condition-out cross-validation (i.e., leave-6-images-out). The parameters of the different methods are optimized with a nested leave-onecondition-out cross-validation within the training set.

Inter-subject experiment The inter-subject analysis relies on subject-specific fixed-effects activations, i.e. for each condition, the 6 activation maps corresponding to the 6 sessions are averaged together. This yields a total of 12 images per subject, one for each experimental condition. We evaluate the performance of the method by cross-validation (leave-one-subject-out). The parameters of the different methods are optimized with a nested leave-one-subject-out crossvalidation within the training set. Spatial correspondence of images within between subjects was assumed after realignment and normalization to the MNI space had been carried out, based on the available anatomical image of each subject.

Results

The results obtained in both intra-subject and inter-subject experiments are given Table . Intra-subject analysis The two sparsity-adaptive approaches VB-MCBR and Gibbs-MCBR outperform the alternative methods, yielding an average explained variance of 0.94 and 0.93 across the subjects. Moreover, their results are more stable across subjects. The SC algorithm and SVR perform poorly.

Inter-subject analysis In this study, TV regression outperforms the non spatially regularized methods, yielding an average explained variance of 0.84, and also more stable predictions. SC also performs well, with an average explained variance of 0.82. ARD, which yields the sparsest model, performs poorly.

Inter-subject analysis -Interpretability of the resulting maps In the case of a linear prediction function, it is easy to look at the voxels weights used in the model. Indeed, these maps can be used to study some large-scale characteristics of the encoding of the cognitive information in brain regions.

In general, one might expect the spatial layout of neural activity to be sparse and spatially structured in the sense that non-zero weights are grouped into connected clusters. Weighted maps showing such characteristics will be called interpretable, as they reflect our hypothesis on the spatial layout of the neural activity. We give Fig. 1 the images obtained with the different methods. From a neuroscientific point of view, the regions are concentrated in the early visual cortex. Indeed, the processings of visual information about sizes are performed in early occipital cortex, with some extent in more parietal regions [START_REF] Eger | Graded size sensitivity of object exemplar evoked activity patterns in human loc subregions[END_REF]. We can see that the SC method creates identifiable clusters, yielding a map similar to the Searchlight procedure, but that it also retrieves additional clusters. TV regression also yields weight maps very similar to the maps obtained by a classical brain mapping approach (such as Anova or Searchlight) but more sparse. We can see that SC and TV benefits from the power of a predictive framework similarly to SVR and Elastic net, while providing brain maps similar to classical SPMs. This confirm the results obtained in terms of prediction accuracy, that spatial regularization is a good way to tackle the spatial variability problem in inter-subjects studies [START_REF] Tahmasebi | Quantification of Inter-subject Variability in Human Brain and Its Impact on Analysis of fMRI Data[END_REF][START_REF] Tucholka | Prise en compte de l'anatomie cérébrale individuelle dans les études d'IRM fonctionnelle[END_REF]. On the contrary, voxel-based methods suffer from the inter-subject spatial variability, and do not yield interpretable maps, even when they achieve high prediction accuracy (e.g. Gibbs-MCBR). We can notice that Gibbs-MCBR, VB-MCBR and SVR yield similar maps, that retrieve a part of the spatial structure obtained in brain mapping approaches. However, Elastic net, while achieving high prediction accuracy, yields very sparse map that is difficult to interpret to retrieve the spatial support of the neural activity.

Discussion

Bayesian versus classical discriminative approaches The methods presented can be roughly classified in two groups: Bayesian approaches (e.g. BRR, ARD, Gibbs-MCBR and VB-MCBR) and classical approaches (e.g. Elastic net, SVC or Total Variation framework). In term of prediction accuracy, the two types of approaches performed similarly, with a slight advantage of the Bayesian methods in the intra-subject analysis. An explanation is that such approaches can more finely estimate the regularization parameter which is not restricted to a predefined grid. A finer grid would be possible but would require more computation. In the inter-subject analysis, classical approaches perform slightly better, because the parameter tuning by internal cross-validation makes them less prone to overfit a particular training subset of subjects.

Sparse versus non-sparse discriminative approaches It is interesting to compare the performance of BRR, ARD and MCBR as they perform a Bayesian regularization with different degree of sparsity. In intra-subject, MCBR performs better than BRR and ARD, as the number of classes may be used to adapt the sparsity between the two extremal cases of BRR (no sparsity), and ARD (high sparsity). In inter-subject, MCBR still performs better than BRR and ARD, but there is an increase in the difference of accuracy between the different methods. Thus, it seems promising to develop methods that are able to adapt their sparsity to the specificity of the dataset, yielding a high regularization (intra-subject), and a less drastic regularization (inter-subject).

Impact of spatial regularization

We can see a clear dissociation between intrasubject and inter-subject analyzes. This can be explained by the different intrinsic resolution of spatial information present in intra and inter-subject settings. Indeed, prediction can rely on relatively sparse and fine-grained patterns at the single subject level. On the opposite, in inter-subject settings, it must be robust to misalignments. Such robustness is obtained through spatial regularization, as in supervised clustering and total variation penalization.

Interpretability of the resulting maps

We have seen that spatially regularized methods yield more interpretable maps that other voxel-based methods. Moreover, compared to a state of the art approach for fine-grained decoding, namely the searchlight, SC and TV yield similar maps, but additionally, take into account non-local information and also have the advantages of a predictive framework (e.g. a prediction score corresponding to whole brain). A joint comparison between the prediction accuracies and the resulting maps also shown that it is difficult to choose a method close to a potential ground truth. A good prediction accuracy and interpretable map do not always come together.

Conclusion

In this paper, we compare different prediction algorithms in both intra-and inter-subject analysis. We show that using spatial information within voxel-based analysis with Total Variation regularization, or by creating intermediate structure as parcels, makes it possible to deal with spatial variability, and yields accurate and interpretable results for reverse inference. We also find that Bayesian approaches, by tuning more precisely the level of sparsity, work well for intra-subject analysis. They might however be trapped into local minima more easily, and did not perform very well in inter-subject experiments.

Fig. 1 .

 1 Fig. 1. Maps obtained for different methods studied, in the inter-subject analysis. First row: the two brain mapping approaches Anova (left) and Searchlight (right) found similar regions. Second row: SVR (left) and Elastic net (right) retrieve part of the spatial structure of the predictive regions, similarly to VB-MCBR (third row, left) and Gibbs-MCBR (third row, right). Fourth row: SC and TV create identifiable clusters.

Table 1

 1 1 (p-values are computed using a paired t-test).

				Intra-subject		Inter-subject
		Method	mean ζ std ζ p-val (VB) mean ζ std ζ p-val (TV)
		SVR	0.82	0.07 < 10 -3 **	0.77	0.11 3.10 -2 *
		BRR	0.92	0.02	10 -3 ***	0.72	0.10	10 -3 **
		Enet	0.90	0.02 < 10 -3 **	0.78	0.10 4.10 -2 *
	Sparse	ARD G-MCBR	0.89 0.93	0.03 < 10 -3 ** 0.01 10 -2 **	0.52 0.79	0.33 0.10 3.10 -2 * 10 -2 **
		VB-MCBR 0.94	0.01	-	0.78	0.10	10 -2 *
	Spatial	SC	0.82	0.07 < 10 -3 **	0.82	0.08	0.6
	regularization TV	0.92	0.02 < 10 -3 **	0.84	0.07	-

. Prediction performance: explained variance ζ for different methods. Sparse methods are in light gray, and spatially-regularized methods in dark gray.


Alexis Lechervy

Philippe-Henri Gosselin

Frédéric Precioso

Linear kernel combination using boosting

HAL is

Linear kernel combination using boosting

Alexis Lechervy 1 * , Philippe-Henri Gosselin 1 and Frédéric Precioso 2 1-ETIS -CNRS/ENSEA/Université de Cergy-Pontoise, 6 avenue du Ponceau, F-95000 Cergy-Pontoise -France 2-I3S -UMR7271 -UNS CNRS 2000, route des Lucioles, 06903 Sophia Antipolis -France Abstract. In this paper, we propose a novel algorithm to design multiclass kernels based on an iterative combination of weak kernels in a schema inspired from the boosting framework. Our solution has a complexity linear with the training set size. We evaluate our method for classification on a toy example by integrating our multi-class kernel into a kNN classifier and comparing our results with a reference iterative kernel design method. We also evaluate our method for image categorization by considering a classic image database and comparing our boosted linear kernel combination with the direct linear combination of all features in a linear SVM.

Context

Recent machine learning techniques have demonstrated their power for classifying data in challenging contexts (large databases, very small training sets, huge training sets, dealing with user interactions...). However, emerging problems are pushing these methods to their limits with sereval hundred of image categories to be classified, with millions of images both in training and testing datasets.

A key component in a kernel machine is a kernel operator which computes for any pair of instance their inner-product in some induced vector space. A typical approach when using kernels is to choose a kernel before the training starts. In the last decade, many researches have been focused on learning the kernel to optimally adapt it to the context and propose a computational alternative to predefined base kernels. Such approaches are particularly interesting in the context aforementioned of huge image databases with hundreds of object categories. In fact, it is pretty much unlikely that a unique base kernel would be adequate to separate all categories while grouping all data from the same category. For instance, Gehler et al. [START_REF] Gehler | On feature combination for multiclass object classification[END_REF] propose to linearly combine several base kernels in order to improve the performance of the major kernel function hence designed for multi-class supervised classification.

In this paper, we propose to design a linear combination of weak base kernels using the boosting paradigm, similarly to [START_REF] Crammer | Kernel design using boosting[END_REF]. However, we focus on a strategy for multi-class learning using many different features.

Before designing a method for the combination of base kernels, it is necessary to define a target kernel K ⋆ that reflects the ideal case on a given training set. In the case of two-class context, this kernel can be defined by K ⋆ (x i , x j ) = 1 if the training samples x i , x j are in the same class, -1 otherwise. This can be expressed on a training set as the Gram matrix K ⋆ = LL ⊤ , where L i is the class label of the i th training sample. Then, the design of the kernel combination K(., .) is driven by the optimization of a criterion between the Gram matrix of kernel combination K and target kernel K ⋆ . Several criteria have been proposed among which class separability [START_REF] Wang | Learning kernel parameters by using class separability measure[END_REF] and data centering [START_REF] Meila | Data centering in feature space[END_REF]. The most classic one is probably the kernel alignment proposed by Cristianini et al. [START_REF] Cristianini | On kernel target alignement[END_REF] which is defined as the cosine of the angle between the two Gram matrices K 1 and K 2 of two kernels k 1 and k 2 .

In the context of multi-class classification, the definition of target kernel is not straightforward. Let L be the n X × n C matrix of annotations that is for n X training samples and n C classes, with L ic = 1 if the i th training sample belongs to c class and -1 otherwise. Then a naive target kernel can be defined by the Gram matrix K L = LL ⊤ . In the following we denote the outer-product XX ⊤ by K X . A first improvement of this target kernel is the introduction of centering, which accounts the high unbalance of multi-class context. Thus, it is recommended to center all Gram matrix (target and combination) using the centering matrix H, and centered Kernel-Target Alignment A H as in [START_REF] Kawanabe | A procedure of adaptive kernel combination with kernel-target alignment for object classification[END_REF]. Furthermore, Vert [START_REF] Vert | Designing a m-svm kernel for protein secondary structure prediction[END_REF] proposed a solution that handles the case of the classes with correlations.

Our learning method based on boosting is presented in section 2. In the following section 3, we present a target for weak learner. The section 4 presents some experiments of classifying both toy data and real data from a standard image database. We then conclude and present the perspectives of this work.

Linear kernel combination using boosting

To overcome the inter-class dependency, we propose to consider the matrix Q of the QR decomposition of HL. We only select the columns where the diagonal element of R is not zero. Thus Q is a n X × n C full rank matrix, assuming that classes are independent. Our target Gram matrix is then defined as K Q = QQ ⊤ . The specific form of this target matrix is further exploited to find the optimal boosting increment (i.e. the kernel evolution direction towards the next best major kernel alignment). Furthermore, as we will see in the next section, the properties of the QR decomposition ensure the convergence of our strategy.

The second contribution is a new boosting method for kernel design. We design a kernel function K(., .) as a linear combination of base kernel functions k t (., .):

K T (x i , x j ) = T t=1 β t k t (x i , x j )
where x i is the feature vector for training sample i, for instance colors and textures in the case of images. We consider base kernel functions defined by k t (x i , x j ) = f t (x i )f t (x j ), where f t (.) is a function built by a weak learner.

In order to build the combination, we work on finite matrix on a given training set X, which leads to the following expression, with f t = f t (X) and

F t = (f 1 f 2 . . . f t )β 1 2 : K t = t s=1 β s f s f ⊤ s = F t F ⊤ t
In the following we mix functions (f , K,...) and their values on the training set (written in bold f , K...).

We select base kernels iteratively in a boosting scheme:

K t = K t-1 + β t f t f ⊤ t ⇔ F t = (F t-1 β 1 2 t f t ) where β t , f t = ζ(F t ) is the result of the problem solved by ζ : ζ(F) = argmax β>0,f A H (FF ⊤ + βff ⊤ , QQ ⊤ )
For a f given, the optimal β can be solved analytically by methods of linear algebra. We build f using least mean squares (LMS) and a target function presented in the following section.

Weak learners optimal target

In order to train weak learners, we need to choose a target function f ⋆ , a function that leads to the best alignment. In the case of two-class context, it can be defined by f ⋆ (x i ) = 1 if training sample i is in the first class, -1 otherwise. However, in the case of multi-class context, this not obvious, since we need to spread each class data around equidistant centers [START_REF] Vapnik | Estimation of Dependences Based on Empirical Data[END_REF][START_REF] Vert | Designing a m-svm kernel for protein secondary structure prediction[END_REF][START_REF] Gosselin | Feature based approach to semi-supervised similarity learning[END_REF].

We propose to consider the centers of (orthonormalized) classes in the space induced by the current combination kernel

K t = F t F ⊤ t : G t = Q ⊤ F t
The rows of G t are coordinates of class centers.

The idea of our method is to move each center to make it equidistant from others. In [START_REF] Vapnik | Estimation of Dependences Based on Empirical Data[END_REF], Vapnik states that the largest possible margin is achieved when the c vertices of (c -1)-dimensional unitary simplex are centered onto the origin. A sufficient means to achieve this properties is to build c orthonormal vertices, whose projection on a (c-1) dimension space is the unitary simplex. In our case, that means that an ideal target set of class centers

G ⋆ t is such that G ⋆ t (G ⋆ t ) ⊤ is proportional to the identity matrix Id c,c .
If we apply the Cauchy-Schwarz inequality to the alignment, we can show a similar observation:

A H (K t , QQ ⊤ ) = 1 ⇐⇒ ξ(G t )ξ(G t ) ⊤ = Id c,c with ξ(G) = QQ ⊤ F HF(HF) ⊤ F G
The aim is to find weak learners that lead to this identity. In other words we are looking for a function f ⋆ such as: As we proceed iteratively, we can only focus on the new column g ⋆ of G ⋆ t = (G t g ⋆ ). A good candidate is:

Id c,c -ξ(G t )ξ(G t ) ⊤ > Id c,c -ξ(G ⋆ t )ξ(G ⋆ t ) ⊤
g ⋆ = 1 -λ QQ ⊤ F HF T (HF T ) ⊤ F v
where λ is the smaller eigen value of G t G ⊤ t and v the eigenvector. Thanks to this result, we can select the target function f ⋆ for weak learners. It can be shown that f ⋆ = Qg ⋆ always leads to the convergence of the whole algorithm.

Experiments and results

In a first experiment (Figure 1) we illustrate the convergence of the method to a regular (#class -1)-simplex.

We consider a toy dataset with 2 classes and 200 examples per class (100 for training and 100 for testing). We use a pool of 10 features of 2 dimensions. For each class c and each feature f , a center C c,f is picked up at uniformly random in [0, 1] 2 . Each example is described by the 10 features with a gaussian with 0.5 standard deviation centered on C c,f . For visualisation we use PCA on F t and select the two first dimensions. After 2 iterations Figure 1 (a), our algorithm separates each class but the distances between classes are uneven. After 45 iterations Figure 1 (b), the classes barycenters converge to the vertex of a regular 2-simplex.

In a second experiment, we compare experimentally our method with a recent Multiple Kernel Learning algorithm proposed by Kawanabe et al. [START_REF] Kawanabe | A procedure of adaptive kernel combination with kernel-target alignment for object classification[END_REF] which computes the optimal linear combination of features. We consider a toy dataset with 10 classes and 60 examples per class (30 for training and 30 for testing). We use a pool of 25 features of 16 dimensions. For each class c and each feature f , a center C c,f is picked up at uniformly random in [0, 1] 16 . Each example is described by the 25 features with a gaussian with 0.5 standard deviation centered on C c,f . Fig. 3: Average precision in % (VOC2006) for linear SVM Figure 1(c) shows that the alignment of our method increases at each iteration on both training and testing data. When comparing the two methods with respect to their alignment results for the same dataset and the same features, their alignment is 0.772469 on train and 0.773417 on test while our method, as seen on Figure 1(c), reaches after 180 iterations an alignment of 0.8507977 on train and of 0.8034482 on test.

Both methods have linear complexity in the number of training samples but Kawanabe et al. approach [START_REF] Kawanabe | A procedure of adaptive kernel combination with kernel-target alignment for object classification[END_REF] has a quadratic complexity in the number of features while our method is linear.

We have also compared our features and Kawanabe features in a multi-class classification context. On the same second synthetic dataset, we classified test data with k-nearest neighbor classifier (kNN) for different k (Figure 2). Our method outperforms Kawanabe et al. when considering fewer neighbors. It aggregates more examples per class.

In a third experiment, we have also compared our method on real data. We evaluate the performance of our algorithm on Visual Object Category (Voc)2006 dataset. This database contains 5304 images provided by Microsoft Research Cambridge and Flickr. Voc2006 database contains 10 categories (cat, car, motorbike, sheep ...). All images can belong to several categories. There are two distinct sets, one for training and one for testing with 9507 annotations. We create our weak kernels from 8 initial features: normalized (L2) histograms of 16, 32, 64, 128-bins for both color CIE L*a*b and quaternion wavelets.

Then we use linear SVM (normalized with L2) to compare the features extracted from the final F matrix with the initial features. We have also evaluated the performance of each extracted feature form F again a feature concatenating all 8 initial features (Figure 3). For all classes our methods reaches higher average precision.

We numerically assess the performance of our method on Oxford Flowers 102 [START_REF] Nilsback | Automated flower classification over a large number of classes[END_REF]. As the authors of this base [START_REF] Nilsback | Automated flower classification over a large number of classes[END_REF], we use four different χ 2 distance matrices to describe different properties of the flowers. Results show that our method improves the performance from 72.8% [START_REF] Nilsback | Automated flower classification over a large number of classes[END_REF] to 77.8%.

Conclusion

In this paper, we propose a new algorithm to create a linear combination of kernels for multi-class classification context. This algorithm is based on an iterative method inspired from boosting framework. We thus reduce both the computation time of final kernel design and the number of weak kernels used.

Considering the QR decomposition leads to a new solution to address the problem of inter-class dependency and provides quite interesting properties to develop an interactive method.

The proposed solution is linear in the number of training samples.

Our method shows good results both on a toy dataset when compared to a reference kernel design method and on a real dataset in image classification context. We are currently working on a generalization of our method to collaborative learning context. Indeed, the same algorithm can target a kernel matrix for collaborative learning by considering that initial annotation matrix stores all previous retrieval runs.

Fig. 1 :

 1 Fig. 1: Converge of the method to regular 2-simplex in 3 classes case after 2 steps(a) and 45 steps(b). Evolution of alignment for 10 classes(c)

  k

* Thanks to DGA agency for funding.


Alexis Lechervy
email: alexis.lechervy@ensea.fr

Philippe-Henri Gosselin
email: gosselin@ensea.fr

Frédéric Precioso
email: precioso@polytech.unice.fr

BOOSTING KERNEL COMBINATION FOR MULTI-CLASS IMAGE CATEGORIZATION

Keywords: Image databases, Machine learning algorithms, Boosting

In this paper, we propose a novel algorithm to design multiclass kernel functions based on an iterative combination of weak kernels in a scheme inspired from boosting framework. The method proposed in this article aims at building a new feature where the centroid for each class are optimally located. We evaluate our method for image categorization by considering a state-of-the-art image database and by comparing our results with reference methods. We show that on the Oxford Flower databases our approach achieves better results than previous state-of-the-art methods.

INTRODUCTION

Recent machine learning techniques have demonstrated their power for classifying data in challenging contexts (large databases, very small training sets or on the opposite, huge training sets, dealing with user interactions...). However, emerging problems are pushing these methods away their limits with several hundred of image categories to be classified, with millions of images both in training and testing datasets.

Different frameworks have been proposed to address image categorization task. One can consider kernel frameworks based on SVM [START_REF] Chang | Statistical learning for effective visual information retrieval[END_REF] or based on boosting [START_REF] Freund | An efficient boosting algorithm for combining preferences[END_REF]. In this article, we propose a new method combining these two frameworks: we design a scheme inspired from boosting framework to build kernel functions which are then integrated in a SVM.

Different techniques have been proposed, in the literature, to linearly combine several (minor) kernels leading to improved performances of the hence designed major kernel. Among these methods, the last decade has seen an important increase of interest of methods which try to learn the combination of kernels in a so-called Multiple Kernel Learning framework (MKL) [START_REF] Jordan | Multiple kernel learning, conic duality, and the smo algorithm[END_REF]. In the the classic MKL framework, methods aim at jointly optimizing the weight of each kernel Thanks to DGA agency for funding. in the combination. However, the two currently most powerful MKL methods are not following this exact scheme: Gehler and Nowozin [START_REF] Gehler | On feature combination for multiclass object classification[END_REF] proposed a boosting-oriented scheme optimizing alternately the combination weights and the combined kernels which reached the best performance on classic MKL evaluation databases. But more recently, Orabona et al. [START_REF] Orabona | Ultra-fast optimization algorithm for sparse multi kernel learning[END_REF] achieve the best current performances by proposing to change the MKL problem formulation.

In the proposed method, we are inline with the works done by Crammer et al. [START_REF] Crammer | Kernel design using boosting[END_REF], Kawanabe et al. [START_REF] Kawanabe | A procedure of adaptive kernel combination with kernel-target alignment for object classification[END_REF] and Gehler et al. [START_REF] Gehler | On feature combination for multiclass object classification[END_REF]. We propose to design a linear combination of base kernels using the boosting paradigm, similarly to [START_REF] Crammer | Kernel design using boosting[END_REF]. However, we focus on a strategy for multi-class learning using many different features. We also take inspiration from [START_REF] Kawanabe | A procedure of adaptive kernel combination with kernel-target alignment for object classification[END_REF] taking advantage of kernel centered alignment to compute the combination weights but they compute all the combination weights at once while we do it iteratively. We consider at each iteration their combination weight computation approach in order to determine the weight in the combination of the major kernel at previous iteration and the weak kernel to be added.

In section 2, we introduce the context of linearly combine kernels inside a boosting framework. In section 3 we develop in details our method. Section 4 presents experiments on image categorization in referenced database for the MKL context and show that our method outperforms the current stateof-the-art approaches. Then we conclude and present of some perspectives of this work.

LINEAR KERNEL COMBINATION USING BOOSTING

Linear combination

The aim of this paper is to design a kernel function K(., .) as a linear combination of base kernel functions k t (., .):

K T (x i , x j ) = T t=1 β t k t (x i , x j ),
where x i is the features for image i.

This combination is performed to match a target kernel function K ⋆ , whose values K ⋆ are sampled on a given training set. In the following, we mix functions (f , K,...) and their values on the training set (written in bold f , K...). The combination is evaluated thanks to Kernel Target Alignment, a popular criterion for kernel optimization. More specifically, we use the centered alignment which allows kernel matrices with unbalanced classes [START_REF] Kawanabe | A procedure of adaptive kernel combination with kernel-target alignment for object classification[END_REF]:

A H (K T , K ⋆ ) = HK T H, HK ⋆ H HK T H HK ⋆ H ,
with H = I -1 n 11 ⊤ and 1 the vector full of ones.

Target kernel

The target kernel matrix K ⋆ is defined from the n×n c ′ matrix L of labels on the training set (1 for positive image and -1 otherwise):

L i,c = label of image i for category c.
To overcome the dependence between classes, we propose to consider the matrix Q of the QR decomposition of HL. We select only the columns where the diagonal element of R is not null. Thus Q is a n × n c full rank matrix, assuming that all classes are independent. Our target Gram matrix is then defined as K ⋆ = QQ ⊤ . The specific form of this target matrix is further exploited to find the optimal boosting increment (i.e. the kernel evolution direction towards the next best major kernel alignment). Furthermore, it can be shown that the orthonormality of Q ensures the convergence of the learning algorithm.

One rank base kernels and weak functions

We consider base kernel functions defined by k t (x i , x j ) = f t (x i )f t (x j ), with f t (.) a function built by a weak learner. In order to build the combination, we work on a finite matrix on a given training set X, which leads to the following expression::

K t = t s=1 β s f s f ⊤ s = F t F ⊤ t , with f t = f t (X) and F t = β 1/2 ⊙ (f 1 f 2 . . . f t )
where ⊙ is the Hadamard product. The values of weak functions f t (.) will be referred as supervised or semantic features built by the method.

We select base kernels iteratively in a boosting scheme:

K t = K t-1 + β t f t f ⊤ t ⇔ F t = (F t-1 β 1 2 t f t ),
where (β t , f t ) is the solution of the following problem:

(β t , f t ) = arg max β>0,f A H (FF ⊤ + βff ⊤ , QQ ⊤ ). (1)

Weak targets

In order to train weak learners, we need to choose a target function f ⋆ , which leads to the best alignment A H (FF ⊤ , QQ ⊤ ).

In a two-classes context, it can be defined by f ⋆ (x i ) = 1 if sample i is in the first class, -1 otherwise. However, in the case of multi-class context, this is not obvious, since we need to spread each class data around equidistant centers [START_REF] Vapnik | Estimation of Dependences Based on Empirical Data[END_REF].

We propose to consider the centers of (orthonormalized) classes in the space induced by the current combination kernel

K t = F t F ⊤ t : G t = Q ⊤ F t .
The rows of G t are the coordinates of each class center. The idea of our method is to move each center to make it equidistant from others. In [START_REF] Vapnik | Estimation of Dependences Based on Empirical Data[END_REF] section 10.5, Vapnik states that the largest possible margin is achieved when the n c vertices of (n c -1)-dimensional unitary simplex are centered onto the origin. A sufficient mean to achieve this properties is to build c orthonormal vertices, whose projection on a (n c -1) dimension space is the unitary simplex. In our case, that means that an ideal target set of class centers

G ⋆ t is such that G ⋆ t (G ⋆ t ) ⊤ is proportional to the identity matrix Id c,c .
In next section we present an algorithm to compute these target functions.

ALGORITHM

Initialization

The first weak function target f ⋆ 0 is defined as the first column of Q. It is equivalent to set the first target centroids g ⋆ as one of the column of identity matrix. Then, we build a set of weak function candidates using least mean squares (LMS). The first weak function f 0 is the candidate that maximizes the centered alignment A H (ff ⊤ , QQ ⊤ ).

Iterations

For each iteration t, we need to define a weak function target f ⋆ t . The first t < n c target function are defined as columns of matrix Q, f ⋆ t = q t . Then target functions for t ≥ n c is defined as f t ⋆ = Qg ⋆ t , with g ⋆ a vector based on the smallest eigenvalue λ and corresponding eigenvector v of G t G ⊤ t :

g ⋆ = 1 -λ n c HF t F ⊤ t H v.
Then, we build a set of weak function candidates using least mean squares (LMS), and the weak function f t is the one that maximize Eq. (1).

For a better generalization, we first evaluate the weak candidate function built using low dimension features. In where no one of them increases enough the centered alignment, we evaluate weak function candidates built using features with a larger dimension. This process is repeated until a relevant one is found.

Final supervised features

Final supervised features F T on the training set is then the concatenation of weighted weak functions values, F T = (β

1 2 t f t ) t .
Let us notice that we can compute these supervised feature for any image, as we learned functions f t (.) that can be evaluated outside the training set.

EXPERIMENTS

In this paper, we test our method on Oxford Flower17 and Flower102 databases [START_REF] Nilsback | Automated flower classification over a large number of classes[END_REF]. In the following experiments, we use the same χ2 distance matrices D provided by the authors [START_REF] Nilsback | Automated flower classification over a large number of classes[END_REF]. For each matrix, we compute a new matrix X = e -µ f D , where µ f is set to be 1 over the mean value of the χ 2 distances over all the training images.

Then we reduce the dimension of X by applying a PCA. Our method then produce supervised features thanks to the training set and these supervised features are further used in a Gaussian Kernel SVM for each category.

Performances are measured with the same protocol as in [START_REF] Nilsback | Automated flower classification over a large number of classes[END_REF].

Oxford Flower 17

Oxford Flower 17 is a database of 17 categories with 80 flower images per category. The dataset is made of 3 predefined random splits and each split consists in 3 sets: training (40 images per class), validation (20 images per class) and test (20 images per class). There are both large intra-category variance owing to appearance variations and small intercategory variance owing to partial appearance similarities between categories.

In our experiments, we use 7 χ 2 -distance matrix provided by the authors of this base1 to describe different properties of the flowers: colour, shape and texture histogram, HSV histogram, SIFT inside the foreground region, SIFT on the foreground region boundary and Histogram of Oriented Gradients.

We run our method with different settings. We build our classifier with features of different granularity. Initially, the coarsest features (i.e in our case the feature with low dimensions) are tried first. If the boosting algorithm fails to find a good kernel from these, then we consider more refined features with more information and detail. In practice, after PCA we keep features of 16, 32, 64, 128, 256, 512 and 1024 dimensions and we start our selection with low dimensions features. The results of this strategy are given in Table 2. This table gives the classification rate for several T and several max dimensions of feature. For example in case T=300 and dim feature = 128 , we have selected up to 300 weak functions and the features have 16, 32, 64 or 128 dimensions. The method may stop by itself before reaching the desired number of iterations because no available feature increases enough the alignment. We can see that the average precision gets higher as the dimension of initial features and the number of iterations both increase.

We reported results from the literature [START_REF] Mikolajczyk | Augmented kernel matrix vs classifier fusion for object recognition[END_REF] in Table 3. Our method outperforms these reference results, and gets lower standard deviation. This standard deviation is computed thanks to the 3 data split predefined by the authors of the dataset.

Oxford Flower 102

Oxford Flower 102 database contains 8189 images divided into 102 flower classes (cf. Fig. 1). Each class consists in about 40 up to 250 images and the distribution of the images between classes is none uniform. The dataset is divided into three sets: training, validation and test. The training and validation sets consist in 10 images per class (1030 images) while the test set has 6129 images with at least 20 images per class.

Like [START_REF] Nilsback | Automated flower classification over a large number of classes[END_REF] we use four different χ 2 distance matrices provided 2 to describe different properties of the flowers: HSV histogram, SIFT inside the foreground region, SIFT on the foreground region boundary and Histogram of Oriented Gradients.

In Fig. 4, we compare our method with the methods from the literature and mainly to the reference method of the stateof-the-art [START_REF] Mikolajczyk | Augmented kernel matrix vs classifier fusion for object recognition[END_REF]. As one can see, our method provides a significant gain.

CONCLUSION

In this paper, we propose a novel algorithm to design multiclass kernel functions based on an iterative combination of weak kernels in a scheme inspired from boosting framework. The proposed method aims at building supervised features where the centroid for each class are optimally located. Experiments have been carried on state-of-the-art image databases Flower17 and Flower102. We compare our results with the reference methods, and showed that better results can be achieved using our method. We are currently working on a generalization of our method to collaborative learning context. Indeed, the same algorithm can target a kernel matrix for collaborative learning by considering that initial annotation matrix stores all previous retrieval runs.

  ± 2.70 78.24 ± 2.70 78.24 ± 2.70 78.24 ± 2.70 78.24 ± 2.70 78.24 ± 2.70 78.24 ± 2.70 100 83.73 ± 2.45 83.73 ± 2.45 83.73 ± 2.45 83.73 ± 2.45 83.73 ± 2.45 83.73 ± 2.45 83.73 ± 2.45 200 84.12 ± 1.35 84.8 ± 1.80 85.78 ± 2.67 85.88 ± 2.56 85.88 ± 2.56 85.88 ± 2.56 85.88 ± 2.56 300 -84.9 ± 1.67 86.27 ± 2.25 85.88 ± 2.33 86.47 ± 2.22 86.47 ± 2.22 86.47 ± 2.22 400 --86.67 ± 2.47 86.27 ± 2.64 87.55 ± 1.80 87.55 ± 1.80 87.55 ± 1.80 500 ---86.47 ± 2.84 87.94 ± 1.47 88.33 ± 1.11 88.33 ± 1.11 Fig. 2: Classification rate on Flower 17 for several T and several dimensions of feature

	Blackberry lily Bolero deep blue Bougainvillea	Bromelia	Camellia	Columbine	Orange dahlia	Gazania
				Fig. 1: Examples of categories for flower 102
		16	32	64	128	256	512	1024
	50	78.24 MKL [9]	72.8
		MKL [10]	73.4
		NLP-β	75.7
		NLP-νMC 73.4
		NLP-B	73.60
		MKL-prod 73.8
		MKL-avg	73.4
		Proposed	77.8
	Fig. 4: Classification rate on Oxford Flower 102

http://www.robots.ox.ac.uk/~vgg/data/flowers/17/index.html

http://www.robots.ox.ac.uk/~vgg/data/flowers/102/index.html


V F Mota

E A Perez
email: eder.perez@ice.ufjf.br

M B Vieira

L M Maciel

F Precioso
email: precioso@polytech.unice.fr

P H Gosselin
email: gosselin@ensea.fr

A tensor based on optical flow for global description of motion in videos

Keywords: Global descriptor, Orientation tensor, SVM, Optical Flow, Motion modeling

Motion is one of the main characteristics that describe the semantic information of videos. In this work, a global video descriptor based on orientation tensors is proposed. This descriptor is obtained by combining polynomial coefficients calculated for each image in a video. The coefficients are found through the projection of the optical flow on Legendre polynomials, reducing the dimension of per frame motion estimations. The sequence of coefficients are then combined using orientation tensors. The global tensor descriptor created is evaluated by a classification of the KTH video database with a SVM classifier.

I. INTRODUCTION

Human action recognition in videos has several applications. Most works approach this problem by a motion analysis and a representation step. Starting from the interesting work of modeling the apparent motion by a projection onto a polynomial basis [START_REF] Druon | Modélisation du mouvement par polynômes orthogonaux : application à l'étude d'écoulements fluides[END_REF], we propose a new global motion descriptor, based on orientation tensors, for the problem of action recognition in videos. Tensors have been largely used to compute optical flow, but its use is still growing on motion analysis and there are few works that use tensor as a descriptor for human action recognition [START_REF] Kihl | Human activities discrimination with motion approximation in polynomial bases[END_REF], [START_REF] Jia | Tensor analysis and multi-scale features based multi-view human action recognition[END_REF], [START_REF] Khadem | Appearance-based action recognition in the tensor framework[END_REF].

The main contribution of this work is the introduction of a video motion indexing scheme based on the modeling of optic flow in video stream. Thus, using statistics on the model hyperparameters, we are able to build a signature (or feature vector) whatever is the width of frames duration of the video. More specifically, we use Legendre polynomial coefficients (which only depends on the chosen degree) to code a per frame optical flow in a tensor that is then accumulated in time defining a global descriptor for video motion.

A. Related work

Through the use of spectral features of a tensor, Jia et al [START_REF] Jia | Tensor analysis and multi-scale features based multi-view human action recognition[END_REF] proposes a method for action recognition based on multiresolution features. A series of silhouettes are transformed into an image called Serial-Frame from which are extracted features for the construction of an eigenvalues and eigenvector space called Serials-Frame Tensor. Analyzing this space, they can separate useful information to recognize different types of action. The video dataset created by authors has 4 different types of actions: running with a tool lifted, walking, running while hitting with bare hand and walking while waving hand.

Also using the idea of silhouettes, Khadem and Rajam [START_REF] Khadem | Appearance-based action recognition in the tensor framework[END_REF] create a set of silhouettes to form a third order tensor comprising three modes: pixels, actions and people. Then they apply a single value decomposition which allows the descriptor to be used to compare the videos of Weizmann database. They argue that this approach is better than the classic principal component analysis.

We can see that these tensor descriptors are all local ones. In general, local descriptors for human action recognition are more explored and achieve greater recognition rates. Hence, there are few references about global descriptors for human action recognition.

Zelnik et al presents a global descriptor based on histogram of gradients [START_REF] Zelnik-Manor | Event-based analysis of video[END_REF]. This descriptor is applied on the Weizmann video database and is obtained extracting multiple temporal scales through the construction of a temporal pyramid. To calculate this pyramid, they apply a lowpass filter on the video and sample it. For each scale, the intensity of each pixel gradient is calculated. Then, a histogram of gradients is created for each video and compared with others histograms to classify the database.

In order to apply a global descriptor on the KTH database, Laptev et al [START_REF] Laptev | Local velocityadapted motion events for spatio-temporal recognition[END_REF] apply the Zelnik descriptor [START_REF] Zelnik-Manor | Event-based analysis of video[END_REF] in two different ways: using multiple temporal scales as the original and using multiple temporal and spatial scales.

A more recent approach is presented in Kihl et al [START_REF] Kihl | Human activities discrimination with motion approximation in polynomial bases[END_REF] and is based on vector fields comparison. The vector fields are obtained through the projection of optical flow on an orthogonal polynomial basis of a given degree which gives a polynomial approximation of it [START_REF] Druon | Modélisation du mouvement par polynômes orthogonaux : application à l'étude d'écoulements fluides[END_REF]. From the similarity measure, they can retrieve a vector field within video sequences. This similarity measure is based on the covariance of the highest eigenvalues of a tensor created from the vector field. This descriptor is tested on videos of motions and can be used for interactive user interfaces. Our work is based on the same idea of polynomial approximation, however our work is not related to the descriptor created by them. They propose a per frame motion descriptor. We propose a tensor descriptor for the whole video or any sequence of frames. Another difference is that they propose a similarity that is based in a tensor decomposition. In our approach, the tensor itself is the descriptor and they are compared by a simple L2 norm.

II. PROPOSED METHOD A. Optical flow modeling using polynomials

The basic idea of a polynomial based model is to approximate a vector field with a linear combination of orthogonal polynomials [START_REF] Druon | Modélisation du mouvement par polynômes orthogonaux : application à l'étude d'écoulements fluides[END_REF], [START_REF] Kihl | Human activities discrimination with motion approximation in polynomial bases[END_REF]. Let us define F an optical flow:

F : Ω ⊂ R 2 → R 2 (x 1 , x 2 ) → (V 1 (x 1 , x 2 ), V 2 (x 1 , x 2 )) (1)
where the functions V 1 (x 1 , x 2 ) and V 2 (x 1 , x 2 ) corresponds to the horizontal and vertical displacement of the point

(x 1 , x 2 ) ∈ Ω.
This optical flow is then approximated by projecting the displacement functions onto each polynomial P i,j , which belong to an orthogonal basis, as such Legendre basis.

In that way, it reduces the dimension of the optical flow field. Thus, we can express F = ( Ṽ 1 (x 1 , x 2 ), Ṽ 2 (x 1 , x 2 )), using a basis of degree g, as:

Ṽ 1 (x 1 , x 2 ) = g i=0 g-1 j=0 ṽ1 i,j P i,j Ṽ 2 (x 1 , x 2 ) = g i=0 g-1 j=0 ṽ2 i,j P i,j (2)
where

ṽ1 i,j = Ω V 1 (x 1 , x 2 )P i,j ω(x 1 , x 2 )dx 1 dx 2 ṽ2 i,j = Ω V 2 (x 1 , x 2 )P i,j ω(x 1 , x 2 )dx 1 dx 2 (3)
It is important to note that the number of polynomials which composes a basis of degree g is:

n g = (g + 1)(g + 2) 2 (4)

B. Orientation tensor: coding frame coefficients

An orientation tensor is a representation of local orientation which takes the form of an N × N real symmetric matrix for N -dimensional signals [START_REF] Westin | A tensor framework for multidimensional signal processing[END_REF].

Given the vector v with N elements, it can be represented by the tensor T = v v T . It is desired that the eigenvector with the largest eigenvalue of the tensor points out the dominant direction of the signal. A signal with no dominant direction is represented by an isotropic tensor, which is a sphere on three dimensions. It is important to note that the well known structure tensor is a specific case of orientation tensor [START_REF] Johansson | A theoretical comparison of different orientation tensors[END_REF].

An individual tensor for each frame is created using the coefficients ṽ1 i,j and ṽ2 i,j (Eq. 3) calculated through the polynomial approximation of optical flow.

From the optical flow approximation, we create a vector ṽf for each frame f of the video:

ṽf = [ṽ 1 0,0 , ..., ṽ1 g,0 , ṽ2 0,0 , ..., ṽ2 g,0 ] ∈ R m (5)
Using the vector ṽf , we create an orientation tensor T f = ṽf ṽf T for each frame f of the video, which is a matrix 2n g × 2n g . This orientation tensor captures the covariance information between ṽ1 i,j and ṽ2 i,j . It carries only the information of the polynomial of frame f .

C. Global tensor descriptor: series of frame tensors

The main idea of our work is to present a new global motion descriptor. We have to express the motion average of consecutive frames using a series of tensors. This can be achieved by b a T f using all video frames or an interval of interest. By normalizing T f with a L2 norm, we are able to compare different video clips or snapshots regardless their length or image resolution.

If the accumulation series diverges, we obtain an isotropic tensor which does not hold useful motion information. But, if the series converge as an anisotropic tensor, it carries meaningful average motion information of the frame sequence. The conditions of divergence and convergence need further studies. In our experiments (Sec. III), this basic global descriptor is called ND.

D. Global tensor descriptor: adding the variation of the polynomial coefficients

It is also possible to add more information to the global descriptor. In order to capture the motion variation in time, we can use both the polynomial coefficients (Eq. 3) and an approximation of their first temporal derivative:

∂t ṽf = ṽ1 i,j (f ) -ṽ1 i,j (f -1) ∆t , ṽ2 i,j (f ) -ṽ2 i,j (f -1) ∆t i+j<g (6)
Hence, the new coefficient vector is

ṽf new = [ ṽf , ∂ t ṽf ].
The orientation tensor T f for each frame is computed as above, resulting in a matrix of 4n g × 4n g . The global motion descriptor obtained by accumulation and normalization now captures the rate of motion changes in time.

It is important to note that the accumulated tensor created is symmetric, so we can use only a triangular superior (or inferior) matrix to represent the video, which reduces the number of coefficients of the final tensor descriptor.

In our experiments (Sec. III), this descriptor enhanced with temporal derivative of coefficients is called WD.

III. EXPERIMENTAL RESULTS

Validation set.: To validate our tensor descriptor, we use the KTH video dataset [START_REF] Schuldt | Recognizing human actions: A local svm approach[END_REF]. Although it is a base with relatively simple actions, there are so few studies suggesting global descriptors that it is still a good reference for comparison.

Experimental protocol: . The optical flow is computed by a method described in [START_REF] Augereau | Vectorial computation of the optical flow in color image sequences[END_REF]. This method was chosen because we found experimentally that it computes a more regular optical flow than the one computed by the standard Lucas-Kanade [START_REF] Augereau | Vectorial computation of the optical flow in color image sequences[END_REF]. We run a multiclass classifier using a oneagainst-all strategy and a Bayes criterion for model selection. Each class is modeled using a SVM classifier with a triangular kernel function with Euclidian distance.

Results.: In order to study the performance of our descriptor, we evaluate the two types of descriptors as described in the previous section: basic tensor (ND) (Sec. II-C) and with derivative information (WD) (Sec. II-D). Instead of using the entire optical flow of the video frames, it is also possible to use only the optical flow from a region with most representative motion. Then, we tested a sliding window with fixed dimensions put around the subject who is doing the action. The center of mass of global optical flow gives the center of the window. It works for KTH scenes because they have only one person acting and a nearly static background. Both descriptors can be computed inside this window, thus we have two more variations: basic tensor using window (ND-W), and tensor with derivative using window (WD-W). Table I shows the recognition rates of the basic descriptor using several polynomial degrees (Sec. II-C). Note that higher degree is, better is the recognition rate. The best recognition rate was 81.13% with polynomial degree of 22. The resulting confusion matrix is shown on Table II. As expected, the worst recognition rate is found for jogging. This can be by high similarities between the optical flow of this class and the ones of walking and running. Table III shows the recognition rates for the descriptor coding the derivative of polynomial coefficients. The best recognition was 84.83% with a basis of degree 4. The confusion matrix is shown on Table IV. We can see that frame coherence reduces significantly the mislabelling of jogging, walking and running actions. The speed of the motion (captured by the added derivative) is the main difference between them. The high recognition rate with low polynomial degree is remarkable. With degree 4, the final descriptor for (WD) has only 1830 elements. This is much better compared to the best result for (ND) with 609960 elements for degree 22 (Table II). Table V shows the results for both descriptors using the sliding window. We have tested several sizes of window and the best for KTH database is 60x100 pixels. This size can capture the most representative motions of the scenes.

We can see that the best result was found for the descriptor with derivative and sliding window, which has a recognition rate 86.44% with a basis degree of having 234270 elements. Its confusion matrix is shown in Table VII. A similar performance of 84.83% is obtained by the descriptor (WD) (Sec. II-D) without using the sliding window on the degree 4 with 1830 coefficients. The large amount of coefficients is not advantageous in terms of time processing. The first result took about 10 minutes and the former 89 minutes in a platform based on Intel Core i7-870 2930Mhz with 8Gb of 1333Mhz/DDR3. An overview of recognition rates for all descriptor variants (ND, WD, ND-W, WD-W) in function of the polynomial degree is shown in Figure 1. Table VI shows the confusion matrix for the descriptor without derivatives and with sliding window, which has a recognition rate of 84.71% with a basis of degree 23. Finally, the descriptor proposed does not outperform the best recognition rates with local information found for KTH database as shown in Table VIII ( [START_REF] Kovashka | Learning a hierarchy of discriminative space-time neighborhood features for human action recognition[END_REF], [START_REF] Kim | Tensor canonical correlation analysis for action classification[END_REF]). However, its performance is close to the state-of-the-art methods, without tracking salient objects. Note that it outperforms the best global descriptor we have found in the literature [START_REF] Laptev | Local velocityadapted motion events for spatio-temporal recognition[END_REF]. [

IV. CONCLUSION

In this work, we proposed a global motion tensor descriptor using a polynomial representation of the optical flow. We use Legendre polynomial coefficients to code a per frame optical flow in a tensor that is then accumulated in time. Although our descriptor does not reach the recognition rate found by local descriptors on KTH, we argue that global descriptors can achieve a good balance between descriptor size, recognition rate and time complexity, and should be more investigated.

Our method beats with a 86.83% recognition rate the global descriptor found in literature [START_REF] Laptev | Local velocityadapted motion events for spatio-temporal recognition[END_REF], that has almost 72% of rate based on histogram of gradients, and is a promising motion representation that can be further improved.

The drawback of our method is that larger and complex video datasets require higher degree polynomials to give good classification results. As a consequence, the number of coefficients increases exponentially leading to high time complexity.

In order to improve the recognition rate of our descriptors, we intend to further analyze the spectral characteristics of the proposed orientation tensor. Moreover, we need to study the conditions of divergence and convergence of the tensor accumulation.

Fig. 1 .

 1 Fig. 1. Overview of recognition rates for all four types of descriptors.

TABLE II CONFUSION

 II MATRIX FOR THE ND DESCRIPTOR AND BASIS DEGREE OF 22: FINAL TENSOR DESCRIPTOR WITH 609960 ELEMENTS.

TABLE IV CONFUSION

 IV MATRIX FOR THE WD DESCRIPTOR AND BASIS DEGREE OF 4: FINAL TENSOR DESCRIPTOR WITH 1830 ELEMENTS.

TABLE VIII COMPARISON

 VIII BETWEEN OUR GLOBAL DESCRIPTORS AND OTHER STATE-OF-THE-ART APPROACHES.

ACKNOWLEDGMENT

The authors would like to thank to CAPES and FAPEMIG for funding.