extfieldsofstudy:
Computer Science
Mathematics

text:
Ensemble learning with 3D convolutional neural networks for connectome-based prediction

The specificty and sensitivity of resting state functional MRI (rs-fMRI) measurements depend on pre-processing choices, such as the parcellation scheme used to define regions of interest (ROIs). In this study, we critically evaluate the effect of brain parcellations on machine learning models applied to rs-fMRI data. Our experiments reveal a remarkable trend: On average, models with stochastic parcellations consistently perform as well as models with widely used atlases at the same spatial scale. We thus propose an ensemble learning strategy to combine the predictions from models trained on connectivity data extracted using different (e.g., stochastic) parcellations. We further present an implementation of our ensemble learning strategy with a novel 3D Convolutional Neural Network (CNN) approach. The proposed CNN approach takes advantage of the full-resolution 3D spatial structure of rs-fMRI data and fits non-linear predictive models. Our ensemble CNN framework overcomes the limitations of traditional machine learning models for connectomes that often rely on region-based summary statistics and/or linear models. We showcase our approach on a classification (autism patients versus healthy controls) and a regression problem (prediction of subject's age), and report promising results.

Introduction
Functional connectivity, as often captured by correlations in resting state functional MRI (rs-fMRI) data, has produced novel insights linking differences in brain organization to individual or group-level characteristics. Recently, machine learning models are being increasingly applied to study and exploit individual variation in functional connectivity data [1,2,3]. These models often employ hand-engineered features, such as pairwise correlations between regions of interest (ROIs) and network topological measures of clustering, modularity, small-worldness, integration, or segregation [4,5,6]. The ROIs are usually computed based on a pre-defined atlas or a parcellation scheme. The choice of the ROIs can have a significant impact on downstream analyses [7,8,9].
Brain ROIs can be defined based on macro-anatomical features, cytoarchitecture, functional activations, and/or connectivity patterns [10,11,12,13]. A common approach is to derive the ROIs either based on input from experts and/or using a data-driven strategy on a small number of subjects. Expert-defined ROIs are challenging to standardize across studies [14] and often rely on arbitrary decisions. Data-driven ROIs, on the other hand, can be biased by the selection of the subjects, especially for regions that exhibit large variability across the population. Popular data-driven techniques include clustering, dictionary learning and Independent Component Analysis (ICA) [15,16,9]. Such methods can be sensitive to confounds such as motion, while initialization, optimization, and other algorithmic choices can also significantly influence the results [17]. A parcellation scheme not only defines the boundaries of ROIs, but also restricts the analysis to a certain spatial scale.
Given the arbitrary nature of a chosen parcellation scheme and its impact on predictive models, we hypothesized that machine learning models can benefit markedly from an ensemble strategy that integrates across different scales and ROI definitions. Figure 1 shows a general schematic of our proposed framework. In this work, we conducted a thorough empirical evaluation of different choices for brain parcellations.
Another important factor in connectome-based machine learning pertains to the choice of the classification algorithm. A large body of related work in the literature has focused on simple linear predictive models using vectorized connectivity data. A relatively recent trend is to represent connectome data as a graph (with subjects as nodes in the graph) and harness graph kernels or specialized neural network architectures to build predictive models. Ktena et al. [18] applied spectral graph convolutions in a distance-metric learning framework to train a k-nearest neighbor classifier on connectivity data. In a similar vein, Kawahara et al. [19] proposed the BrainNetCNN architecture that extends convolutional neural networks (CNNs) to handle graph-structured data. While CNNs are motivated via the translation-invariance property of image-based classification problems and thus have achieved tremendous success, the neuroscientific basis of the invariance property exploited by BrainNetCNN remains elusive. Furthermore, this approach works directly with an adjacency matrix derived from the connectome data, while disregarding spatial information. As we discuss below, we propose an alternative representation of connectivity data, which allows us to leverage modern deep learning architectures, like CNNs, to build a prediction model that exploits spatial information.
In this work, we consider two applications: discrimination of autism patients and healthy controls; and regression of age. The first problem is a particularly challenging one. Several previous studies have reported altered functional connectivity patterns in Autism Spectrum Disorder (ASD) pa-tients [20,21,22,23]. While studies using small samples have reported classification accuracies over 75% [24], application of similar models on large heterogeneous datasets, such as ABIDE [25], have shown more modest performance levels over a wide range of connectome preprocessing schemes (accuracies that range 60-67%) [26].
Our main contributions in this paper are: • An extensive evaluation of the influence of brain parcellations on functional connectome-based machine learning models • An ensemble learning strategy for combining predictions from multiple classifiers corresponding to different brain parcellations • An easy-to-implement 3D CNN framework for connectome-based classification • A technique to visualize and compare trained CNN models using saliency maps that reflect the importance of individual voxels for prediction

Dataset
Autism Brain Imaging Data Exchange (ABIDE) is a multi-site consortium aggregating and openly sharing anatomical, functional MRI and phenotypic datasets of individuals diagnosed with ASD, as well as healthy controls (HC) [25]. The first phase of ABIDE (ABIDE-I) collected data from 1,112 individuals, comprising 539 individuals diagnosed with ASD and 573 typical controls across 17 sites. The second phase (ABIDE-II) aggregated 1,114 additional datasets, comprising 521 individuals with ASD and 593 healthy controls across 19 sites.

Preprocessing of fMRI Data
The Preprocessed Connectomes Project (PCP) released preprocessed versions of ABIDE-I using several pipelines [27]. We used the data processed through the Configurable Pipeline for the Analysis of Connectomes (CPAC). This pipeline performs motion correction, global mean intensity normalization and standardization of functional data to MNI space (3x3x3 mm resolution) before the extraction of ROI time series. Among the different strategies in the release, our analysis used data de-noised by regression of nuisance signals including motion parameters, CompCor WM+CSF components, and global signal, followed by band-pass filtering (0.01-0.1Hz).
We preprocessed the ABIDE-II dataset following the same sequence of steps listed for ABIDE-I in CPAC (using the version v1.0.2a). Since manual quality control (QC) was not yet available for ABIDE-II, we performed an automatic QC by selecting those subjects that retained at least 100 frames or 4 minutes of fMRI scans after motion scrubbing [28]. Motion scrubbing was performed based on Framewise Displacement (FD), discarding one volume before and two volumes after the frame with FD exceeding 0.5mm [29].

Cohort selection
In our experiments, we used ABIDE-I subject data that passed manual QC by all the functional raters. This yielded a final sample size of 774 ABIDE-I subjects, comprising 379 subjects with ASD and 395 typical controls. As an independent test dataset, we employed ABIDE-II subjects from sites that participated in ABIDE-I and used the same MRI sequence parameters for data collection. After automatic QC, we ended up with a final ABIDE-II sample size of 163 individuals with ASD and 230 healthy controls. For age prediction, we only considered healthy controls. Furthermore, subjects whose age were more than 3.5 standard deviations away from the median were excluded from the task of age prediction. Table 1 summarizes the dataset characteristics for the two prediction tasks considered in this study.
For our 3D CNN model, described below, the parcellated regions were used as target ROIs to derive the input connectivity features at the voxel level. For the non-CNN benchmark models, also described below, each atlas was used to define a corresponding connectivity matrix which was fed as input to each model after collapsing into a vector. We report results for ensemble learning strategies as well, where we combined the predictions of models corresponding to individual atlases.

Creating stochastic parcellations
Stochastic parcellations were created by Poisson Disk Sampling using the method described in [39]. Given a number of ROIs, this approach divides the gray matter voxels (as defined by a given mask) into roughly equal-sized parcels while ensuring that the parcels do not cross hemisphere boundaries. Stochasticity is introduced in the ROI center locations, and all the remaining voxels are assigned to the closest region center. These centers are kept a minimum distance apart based on the desired number of regions in the parcellation. All parcellations were created in the MNI152 template at a 3mm resolution, same as the resolution of the preprocessed functional data. For creating these parcellations, we relied on a whole brain gray matter mask including sub-cortical structures. To create the mask, we took the union of the gray matter tissue prior provided in the standard MNI152 template and the cortical mantle mask used in [16].

3D Convolutional Neural Network Approach
Here, we present our novel strategy to adopt a 3D CNN architecture for use with connectomic data. The input to the CNN is formed by concatenating voxel-level maps of "connectivity fingerprints", which are represented as a multi-channel 3D volume. Each channel is a connectivity feature, such as the Pearson correlation between each voxel's time series and the average signal within a target ROI. In our implementation, we use both atlas-based and stochastic brain parcellation schemes to define target ROIs. The total number of input channels thus represents the number of ROIs used for creating voxel-level fingerprints. For each parcellation scheme (atlas-based or stochastic), we trained a separate model.
In our experiments, we employed a simple CNN architecture, illustrated in Fig. 2. Our architecture has several convolutional layers, interspersed with max-pooling based down-sampling layers, followed by a couple of densely connected layers. The models were trained with a mini-batch size of 64, until convergence of validation loss. For classification, we used binary crossentropy, whereas for regression we adopted absolute difference as the loss function.
The neural network weights were then optimized via stochastic gradient descent (SGD). The learning rate and momentum for SGD were set to 0.001 and 0.9 respectively. The same architecture and settings were used for all atlases and stochastic parcellations. We note that each atlas is defined on a unique gray matter mask. To ensure that all prediction models (benchmark and proposed) relied on information from the same voxels, the atlas-specific gray matter mask was applied to the voxel-level connectivity fingerprint data before feeding into the proposed convolutional architecture. For stochastic parcellations, the custom gray matter mask as described above was used for masking the fingerprints.

Benchmark Methods
In our experiments, we implemented following benchmark methods.

Ridge Regression
A linear regression model was trained with squared loss and α times the squared norm of the weight vector (See Appendix). For classification, the ground truth labels were encoded as ± 1 for the two output categories. We tested 10 linearly spaced values for the hyper-parameter α in the range [0.1,10] and report the highest cross-validation accuracy. Thus this baseline result reflects an optimistic estimate of performance.

Support Vector Machine
We implemented a standard SVM as a benchmark (See Appendix). We found that a radial basis function (RBF) kernel performed better than a linear model. Thus we report results for the RBF-kernel SVM. The two hyper-parameters (RBF kernel width γ and and misclassification cost weight C) were fine-tuned by maximizing cross-validation accuracy via a grid search. Therefore, as with the ridge classifier, this should be considered as an upper bound on generalization performance. For regression, we implemented the standard SVR scheme with an -insensitive loss function, optimizing for the -tube and penalty parameter of the error term via grid search.

Fully Connected Architecture
The fully-connected neural network (FCN) architecture takes as input functional connectivity estimates between pairs of ROIs, which is vectorized and processed by a feed-forward network. We implemented following architecture, which performed best on ABIDE-I cross-validation: 4 fully connected hidden layers, with 800, 500, 100 and 20 numbers of features and each linear layer followed by an elementwise Exponential Linear Unit (ELU) activation. Dropout regularization parameter was set to 0.2 and applied to each layer during training. For classification, the output node was a sigmoid, and crossentropy loss was used. For age prediction, the sigmoidal output was replaced with a linear activation and absolute difference was used as the loss function. The models were trained with a mini-batch size of 64, until convergence of validation loss.
SGD was used as the optimizer with learning rate and momentum set to 0.01 and 0.9 respectively for classification. For age prediction, a smaller learning rate of 0.002 was used.

BrainNet Convolutional Neural Networks
BrainNet CNN, originally proposed in [19], utilizes specialized kernels to handle connectomic data. Their work described novel edge-to-edge, edgeto-node and node-to-graph convolutional layers that can potentially capture topological relationships between network edges. For BrainNet CNN, we implemented the following architecture that worked best on ABIDE-I crossvalidation: 1 edge-to-node layer with 256 filters, followed by a node-to-graph layer with 128 output nodes and finally a dense layer with single output. A leaky ReLU non-linearity with alpha equal to 0.33 was applied to the output of each layer except the last layer. The activation of the last layer was set to linear and sigmoid for the regression and classification tasks, respectively. Dropout regularization with rate 0.2 was used for the edge-to-node layer. Similar to [19], Euclidean loss was minimized for age regression, whereas cross-entropy loss was used to optimize the classification models. The models were trained for 1000 iterations using SGD with momentum equal to 0.9. The learning rate was set to 0.0005 for age prediction and 0.008 for ASD/Healthy classification. The training curves were monitored for atlases to ensure convergence.

Ensemble Learning
In our experiments, we explored two ensemble learning strategies. The first one is what we call multi-atlas ensemble (or MA-Ensemble). MA-Ensemble averages the predictions of the models of a specific method (e.g., BrainNet CNN) computed using each one of the seven atlases. For classification, the final prediction is computed as the majority vote of the individual binary class predictions. For regression, the ensemble prediction is simply the mean. The second ensemble strategy (SP-Ensemble) averages across the models of a specific method computed using stochastic parcellations. In our experiments, unless stated otherwise, we used 30 stochastic parcellations at each of the following four spatial scales: 110, 160, 200 and 400 ROIs. These scales were chosen in accordance with existing atlases. Thus the SP-Ensemble's prediction was computed based on fusing 120 (30 × 4 scales) models. We also implemented single-scale SP-Ensemble models, which averaged over the 30 parcellations at the same spatial scale.

Visualizing the CNN model
In order to understand the connectivity features captured by the CNN model, we employed the saliency map approach of [40]. This visualization  technique computes the gradient of the output prediction with respect to the input image voxel values, i.e., the 3D volume, using a single backward pass through the trained neural network. We then computed voxel-level saliency as the maximum absolute gradient value across all input channels corresponding to different target ROIs. More formally, consider an input image I, representing the connectivity fingerprints of v voxels with R ROI signals. The saliency weights w R v×R are computed by taking the absolute value of the gradient of neural network output O with respect to the input image, i.e., w = | ∂O ∂I |. In order to obtain the saliency at the voxel level S R v , we take the maximum across all the ROIs, i.e., S i = max 1≤j≤R w ij . Finally, to visualize an ensemble model, we averaged the individual saliency maps that made up the ensemble.

Experiments
In our experiments, we considered two tasks: i) binary classification of autism vs healthy, and ii) age prediction. For each task, we implemented two evaluation schemes. First, we conducted 10-fold cross-validation on the ABIDE-I dataset, so that we could present results that were consistent with previously reported classification results such as [1,26]. Second, we trained each model on the entire ABIDE-I dataset and computed test performance on the independent ABIDE-II set. We report classification accuracy and  the receiver operating curves (ROC), along with corresponding area under the curves (AUC) for each of these scenarios under various combinations of parcellation schemes and prediction algorithms. For age prediction, we report the Mean Absolute Error (MAE). Table 2 shows the cross-validation and independent test performance for different models on the classification problem. Clearly, the proposed 3D CNN approach shows superior prediction performance compared to the benchmark methods, including fully-connected deep neural network (FCN) and BrainNet CNN. In particular, the 3D CNN approach yields the best ABIDE-I crossvalidation accuracy for all parcellation schemes, including the ensembles. Similarly, the SP-Ensemble achieves the best ABIDE-I cross-validation for all algorithms, including the 3D CNN. The ABIDE-II results are in general compatible with the cross-validation results, where the 3D CNN and SP-Ensemble techniques mostly outperform the competition. Figure 3 shows the Receiver Operating Characteristic (ROC) curves for SP-Ensemble models for the different algorithms on the independent ABIDE-II test dataset. We observe that the 3D CNN achieves the largest AUC of 76.86%. Table 3 lists cross-validation and independent test results for the age prediction task. The 3D CNN approach consistently shows superior performance, yielding the best results for all but two parcellation schemes. Sim- ilar to the classification scenario, SP-Ensemble also yields the best crossvalidation and independent test performance values for the majority of the algorithms, including 3D CNN. Overall, the best accuracy is achieved by SP-Ensemble 3D CNN, which yields a mean absolute error of 2.41 years on ABIDE-I cross-validation and 1.94 years on the independent ABIDE-II dataset.

Comparison of stochastic parcellations and atlases
Here, our objective is to conduct a detailed investigation of how the choice of ROIs affects prediction performance for different machine learning (ML) algorithms. For each ML algorithm and each parcellation we have a model trained on the ABIDE-I data, which we then used on the independent ABIDE-II data to quantify prediction accuracy. Figure 4 shows the distribution of accuracy values (estimated with a kernel density model) obtained using stochastic parcellations , while also illustrating the results for each of the atlases and the scale-specific SP-ensembles. The scale-specific SP-Ensemble strategy, as the name implies, averaged the models corresponding to the 30  stochastic parcellations in each scale. We observe that the atlas-based models performed no better than typical stochastic parcellation models, independent of scale and algorithm.
This result presents a revelatory perspective: perhaps we do not need anatomically or functionally derived brain parcellations for machine learning. Furthermore, the proposed SP-Ensemble CNN strategy yielded accuracy results that were about as good as the best scale-specific SP-Ensemble model. Finally, the ensemble models were almost always better than the atlas-based models and they compared favorably against the individual stochastic parcellation models. The same observations can be made for ABIDE-I crossvalidation (see Supplementary Figure S.1).
In above analysis, one potential confound was the different gray matter masks of atlases and stochastic parcellations (SPs). In order to account for this confound, we conducted following analysis. For each of the atlases, we generated 100 SPs using the same gray matter mask as the atlas. We excluded DOS160 because it does not rely on a well-defined gray matter mask and places discontiguous 4.5 mm spherical regions over fixed coordinates in the brain (sampling only 5% of brain voxels). We then trained on each of these SPs using the same hyper-parameters that were found to be optimal for the corresponding atlas. Here, we show the results for ridge regression (the model that was fastest to train), but we obtained similar results for all  other algorithms as well. As can be seen from Figure 5, for most atlases and corresponding gray matter masks, the model trained on the atlas ROIs performed no better than an average SP model. Furthermore, and importantly, the SP-Ensemble (computed by averaging across SPs on the atlas-specific mask) yielded better performance than the atlas models for all atlases.

Visualization
An important goal of machine-learning tools in neuroimaging is to generate novel insights linking imaging biomarkers with disease or phenotypic traits. Visualization techniques for CNNs can help reveal important features used by the model for discriminating between output classes. Figure 6 shows the saliency maps computed for the SP-Ensemble CNN ASD classification and age prediction models. As can be seen from these maps, the precuneus, often considered a core node of the default mode network [41], seems to play a significant role for both prediction problems. However, there are also salient regions that are unique to each problem. For example, the anterior cingulate/ventromedial prefrontal cotex, a region that has been linked to autism [42], was distinctly highlighted for the ASD classification problem. The left parietal cortex was also emphasized for ASD prediction, which is consistent with the laterilized activation observed in this region in Autism patients [43]. On the other hand, for age prediction, the left dorsolateral prefrontal cortex (dlPFC) is a uniquely salient region. The dlPFC is associated with executive functions, such as working memory and abstract reasoning. For working memory, dlPFC's function seems to be age-associated and more lateralized in younger adults [44].

Discussion
In this study, we presented a detailed empirical analysis of how the choice of ROIs can impact the performance of machine learning models trained on functional connectomes. We considered several machine learning algorithms, together with a range of spatial scales and parcellation schemes, including the popular atlas-based techniques and a stochastic approach. Our analysis suggests that using a single atlas for summarizing the connectome data is often sub-optimal for training machine learning models, and significantly more accurate predictions can be achieved with an ensemble approach that averages across models trained with different parcellation schemes. Furthermore, we demonstrated that averaging across stochastic parcellations can achieve very high accuracy values, often surpassing atlas-based models.
Another main contribution of this study is a novel approach to employ a 3D CNN architecture on functional connectivity data. Convolutional neural networks achieve state-of-the-art performance on many image-based prediction tasks, as they take advantage of the full spatial resolution of the data and the translation invariance property of the problem. Our proposed approach treats voxel-level connectivity fingerprints as input channels to a conventional 3D CNN framework. This strategy contrasts with widely used machine learning models that work with functional connectivity data, while ignoring the spatial structure. Our results demonstrate that when tailored for connectomes, CNNs offer a promising opportunity to probe brain networks in disease.
Machine learning practitioners have to make a number of preprocessing choices in extracting connectomic features to analyze. While there is no onesize-fits-all solution across different tasks, in the context of machine learning models of functional connectivity, we present some interesting empirical observations below.

Ensemble learning
The motivation behind using multiple stochastic parcellations for prediction is grounded in the concept of ensemble learning. The core idea is to integrate out a latent variable (i.e., parcels or ROI definitions) from the learning problem [45]. This approach also makes the predictions more robust to the precise parcellation scheme. As shown above, the performance of atlas-based models can vary significantly (∼5-10% for parcellations at the same scale). In such a scenario, ensemble learning over multiple stochastic parcellations can be a robust strategy that yields reliable predictions.

Network granularity
We explored the impact of network granularity on prediction performance of machine learning algorithms for connectomes. Our analysis suggests that better prediction performance can be expected with parcellations at higher granularity. However, no significant differences were observed beyond 400 ROIs. Our evaluations contradict with a previously reported result that a coarser network scale (∼ 100-150 ROIs) is more suitable for autism classification [46]. In their paper, these conclusions were drawn by comparing the performances achieved with a few atlases. However, inferring trends from a small number of atlases can be misleading, since factors like the boundary definitions of structures (cortical/subcortical) or the particular gray matter mask used, will effect results. Stochastic parcellations can control for these confounds and depict unbiased trends across network scales.

Number of gray matter voxels
Our empirical study suggests that there is no direct correlation between the number of voxels in the gray matter mask and a model's prediction performance. However, we do observe that the choice of gray matter mask can impact results. For example, the DOS160 atlas with as few as ∼ 3,039 voxels shows performance no worse than other atlases at the same resolution (HO, EZ, TT and AAL) with ∼ 20x more voxels.

Visualization
Saliency maps provide a valuable visualization strategy to probe deep neural network models. We visualized the saliency maps from 3D CNN models trained on ROIs extracted using both atlases and stochastic parcellations. As shown in Figure 6 and Supplementary Figures S.2 and S.3, these maps are remarkably consistent.
These maps reveal that the precuneus, which is a hub of the default mode network and associated with ASD and age, plays an important role for both prediction problems. There were also uniquely highlighted regions, such as the anterior cingulate/ventromedial prefrontal cotex for ASD classification and the left dorsolateral prefrontal cortex (dlPFC) for age prediction.

Influence of motion
Several studies have shown differences in head motion parameters during fMRI between healthy controls and diseased populations, or between subjects from different age groups [47,48]. This, in turn, can manifest as artifacts in the derived resting-state connectivity [49]. Although our independent test data was motion scrubbed, we performed additional analyses to rule out the confounding effect of motion in classifier decisions. We selected a cohort of 151 ASD subjects with motion-matched healthy controls from our independent dataset and analyzed the correlation of 4 motion parameters with classifier predictions. These include the root-mean-square framewise displacement, mean relative displacement, maximum absolute displacement and the number of micro-movements greater than 0.5mm. These summary statistics were chosen in accordance with previous reports of motion artifacts in rs-fMRI [50]. As shown in Figure 7, no significant correlations were observed between motion variables and the predictions of SP-Ensemble (model average over all atlases). In this motion-matched cohort, classification accuracy of 71.8% was obtained using 3D-CNN.
For our regression task, there was no significant correlation between a subject's age and any of these motion parameters in our cohorts.

Limitations and future work
Throughout our analysis, Pearsons correlation was chosen to measure functional connectivity strength between different brain regions. Several other correlation metrics, including tangent-based and partial correlation have been shown to yield superior classification performance in prior studies [9,26]. While we do not expect this to affect the general conclusions and findings of our study, the choice of the correlation metric still remains an arbitrary decision in any machine learning pipeline for connectomes.
Due to the heavy computational burden required for training multiple deep learning models, we only considered one particular scheme for creating stochastic parcellations, i.e., Poisson Disk Sampling. Alternative strategies for creating random parcellations have also been proposed, for instance, through stochastic sub-division of anatomically derived ROIs into smaller parcels [51]. It is also possible to randomize several other more popular schemes for parcellating the brain, such as, using Ward's clustering on functional data from sub-samples of the population [45] or creating Geometric parcellations with different initializations [13].
While the proposed CNN approach achieves promising accuracy on autism detection and age prediction, there is room for further improvement. We have not yet conducted a comprehensive optimization of the convolutional architecture. Furthermore, there are likely more optimal choices than target ROI-based correlations that are used as input to the model. An interesting alternative would be select random gray matter vertices for connectivity profiling, as proposed in [16]. We envision an end-to-end learning strategy that can enable the optimization of these connectomic features.
Saliency maps provide an appealing visualization technique by mapping the neural network activations back to input voxel space. Several modifications to gradient-based back-propagation have been reported in literature that can potentially highlight more informative features learnt by the model [52,53]. Further, the use of saliency maps need not be restricted to depicting group-averaged discriminative features. Unsupervised learning on saliency maps can provide novel insights into clinical subtypes of disease. It is also important to note that machine learning techniques do not unequivocally provide evidence for the salient features being directly associated with the disease or other target variables. However, when combined with detailed future investigations, they can spur clinical discoveries.

Conclusion
The results presented in our paper showcase the utility of ensemble learning for connectomes. Functional network based prediction models are impacted by several a priori choices, the most pivotal of which is the ROI definition. We demonstrate that ensembles of stochastic parcellations yield predictions that are significantly more robust and accurate compared to single atlas-based approaches. Further, our experiments highlight the potential of convolutional neural network models for connectome-based classification.

Linear Classifiers 7.2.1. Ridge Classifier
Given feature vectors x i for n subjects and the corresponding prediction variables denoted by y i , we approximate the fit using a linear regression model. An L 2 regularization for the weights (w) is added to the mean squared error to yield the following loss function of ridge regression: During classification, the output labels y are encoded as ± 1 for the two output categories to minimize the above loss.

Support Vector Machines
(a) Classification Support Vector Machine Classifiers optimize for a hyperplane with maximum margin between the output classes. This results in a decision function of the form, f (x)= sign(w T x + b). The weights {w, b} are obtained by minimizing the following convex loss function consisting of a data loss component (L D ) and a regularization loss for the weights (L W ), L D is modeled using a hinge loss function, n i=1 max(0, 1 − y i (w T x i + b)) over all n training samples {(x 1 ,y 1 ),...,(x n ,y n )}. L W is modeled using a Euclidean 27 norm, i.e., w 2 . Here, C is a tuning parameter that controls the trade-off between regularization and data loss.
The -Support Vector Regression (SVR) scheme optimizes for a decision function of the form, f (x)=w T x + b, that has at most deviation from the true prediction variables y (allowing for errors when the problem is infeasible). The loss function (L SVR ) can be formulated as, L is traditionally referred to as the -insensitive loss function, and is formulated as n i=1 max(0, |w T x i + b − y i | − ) over all n training samples {(x 1 ,y 1 ),...,(x n ,y n )}. The regularization term (L W ) is modeled using a Euclidean norm, i.e., w 2 . The tuning parameter C controls the trade-off between the regularization (i.e., the flatness of the decision function) and the amount up to which deviations beyond are tolerated. Both the classification and regression problems yield weights w that can be represented completely as a linear combination of the training inputs x i . Thus, w is represented as n i=1 α i x i , and the decision function becomes f (x)= n i=1 α i x T i x + b. This makes it easier to extend SVMs for non-linear decision functions using the kernel technique, i.e., by applying transformations φ(x) that map x to a high-dimensional space and replacing the inner product x i , x with the kernel K(x i , x)= φ(x i ), φ(x) . For our experiments, we observed that the radial basis function kernel, K(x i , x) =exp(-x i −x 2 2σ 2 ), yields the best results among linear, sigmoid and polynomial kernels up to degree 4.

ABIDE-I cross-validation results
In order to ensure a fair comparison with other studies that report 10fold cross-validation performance on ABIDE-I, we report the performance obtained using our benchmark and proposed models (along with the ensemble learning strategy) for both stochastic parcellations and atlases in the form of kernel density plots ( Figure S.1). Clearly, the results and conclusions on ABIDE-I remain consistent with ABIDE-II, with the 3D-CNN ensemble strategy outperforming all the baseline methods.

Saliency maps for individual parcellations
Visualizing the saliency maps for models trained on different brain parcellations can reveal interesting differences in the features captured by these models. We visualized the saliency maps of the 3D-CNN model for individual stochastic parcellations at multiple scales for the task of ASD/HC Classification. As shown in Figure S.2, models trained using distinct parcellation schemes are relying on the same basic underlying connectivity patterns for prediction, with small differences in their information content, that can be utilized efficiently by the ensemble learning scheme. Further, the saliency maps of atlas-based (see Figure S.3) and stochastic parcellation-based models are remarkably similar, suggesting that the connectivity patterns of the same set of voxels are guiding the classifier predictions, irrespective of the precise scheme of ROI extraction.


text:
Manifold Adversarial Learning

Recently proposed adversarial training methods show the robustness to both adversarial and original examples and achieve state-of-the-art results in supervised and semi-supervised learning. All the existing adversarial training methods consider only how the worst perturbed examples (i.e., adversarial examples) could affect the model output. Despite their success, we argue that such setting may be in lack of generalization, since the output space (or label space) is apparently less informative.In this paper, we propose a novel method, called Manifold Adversarial Training (MAT). MAT manages to build an adversarial framework based on how the worst perturbation could affect the distributional manifold rather than the output space. Particularly, a latent data space with the Gaussian Mixture Model (GMM) will be first derived.On one hand, MAT tries to perturb the input samples in the way that would rough the distributional manifold the worst. On the other hand, the deep learning model is trained trying to promote in the latent space the manifold smoothness, measured by the variation of Gaussian mixtures (given the local perturbation around the data point). Importantly, since the latent space is more informative than the output space, the proposed MAT can learn better a robust and compact data representation, leading to further performance improvement. The proposed MAT is important in that it can be considered as a superset of one recently-proposed discriminative feature learning approach called center loss. We conducted a series of experiments in both supervised and semi-supervised learning on three benchmark data sets, showing that the proposed MAT can achieve remarkable performance, much better than those of the state-of-the-art adversarial approaches. We also present a series of visualization which could generate further understanding or explanation on adversarial examples.

I. INTRODUCTION
A DVERSARIAL examples refer to augmented data points generated by imperceptible perturbation of input samples. Being difficult to distinguish from real examples, such adversarial examples could however change the prediction of many of the best learning models including the state-ofthe-art deep learning models [1] [2] [3]. To alleviate such problems, researchers have proposed adversarial training, able to certify both the robustness on adversarial examples and the generalization on original examples. Adversarial training could be used in both supervised and semi-supervised training. In supervised adversarial learning, the data labels are needed to derive the worst perturbation against loss function [4] [5]; in semi-supervised adversarial learning, a virtual adversarial training (VAT) is used by smoothing the output distribution with penalizing the KL-divergence between outputs of adversarial and original examples. VAT achieves state-of the-art performance on both image and text classification [6] [7] [8].
Previous adversarial training methods simply consider how to make the results of prediction worse (in the output space) without considering how the data are robustly represented in a latent space. In general, the latent space is much more informative than the output space. It is hence meaningful if we can design the adversarial learning in the latent space rather than the output space. In this work, we develop a novel model called Manifold Adversarial Training (MAT) in the latent space. We engage an information based regularization, i.e., Maximum Mutual Information (MMI) [9] [10] so as to define a distributional manifold in the latent space. We then apply the adversarial training to smooth such manifold by penalizing the KL-divergence between the distributions of latent features of the adversarial and original examples. The novel framework is trained in an adversarial way: the adversarial noise is generated to rough the distributional manifold, while the model is trained to smooth it to make the latent space more representative. It is similar to traditional Laplacian regularization methods with locality-preserving properties [11] [12]. However, our approach is based on information geometry and the information metric KL-divergence.
To our best knowledge, this is one novel work that learns adversarially both a robust and compact representation in the latent space. It also presents a unified framework in that a simplified MAT could derive a famous recently-proposed discriminative feature learning model [13] [14]. Though developed in the framework of supervised learning, it is straightforward and much easier to be extended in semi-supervised learning. We develop a feasible and efficient training algorithm capable of obtaining remarkably better performance than the existing adversarial training approaches. In particular, we implemented our proposed method on benchmark datasets MNIST, CIFAR-10, and SVHN. Our method achieves in supervised and semi-supervise learning the state-of-the-art performance, much better than the best of the existing counterpart methods.

II. RELATED WORK
It has a long history to use the perturbed examples to regularize the output [15]. Bishop et al. proposed a method to add the Gaussian noise to input samples and showed that it is equivalent to adding the penalty term to original objective function [16]. Dropout can also be treated as random perturbation to prevent from over fitting [15] [17]. The Unified Gradient Regularization Family is proposed to find the worst perturbation to increase the objective function [5]. It approximates the non-convex problem with Taylor series and applies the Lagrange multiplier method to evaluate the worst perturbation. Another similar work proposed by Aman et al. is to perturb the underlying data distribution in a Wasserstein ball [4]. Virtual Adversarial training proposed by Takeru [6] [8] is perhaps most related to our work. It developed a method extending the adversarial training to semi-supervised task by promoting the local smoothness of the output distribution.
Another type of semi-supervised learning methods are based on generative models. Ladder network combined the deep network and auto encoder with connections between two networks at each layer and achieves encouraging results [18]. Triple generative Adversarial Network is proposed to combine Generative Adversarial Network (GAN) with classifier. There are three players, generator, discriminator and classifier playing against with each other [19]. Some Bayesian methods employ variational methods with deep learning [20].
This work is also related to some traditional manifold regularization. Laplacian Eigenmaps was developed to consider how to construct a representation for data lying on a low dimensional manifold embedded in a high-dimensional space [11]. A geometric framework was also proposed to exploit the geometry of the marginal distribution using the information of both labeled and unlabeled examples [12]. Both of these two methods tried to smooth the low dimensional manifold with preserving local neighborhood information.

III. MAIN METHOD
In this section, we first introduce two previous adversarial training methods, adversarial training with l p norm constraint and Virtual Adversarial Training (VAT) which are close related to our proposed approach. Then we describe our proposed Manifold Adversarial Training (MAT). Specifically, we first introduce how to regularize the latent space with Maximum Mutual Information (MMI) and represent the latent features with Gaussian Mixtures. We then define the notion of the low dimensional submanifold of latent space and the smoothness of statistic manifold. After that we describe the major framework and present a series of theory to derive the practical optimization algorithm. Finally, we also briefly conduct a computational analysis. We now give a set of notations. The input set is represented as . Let x i ∈ R I and y i ∈ R O where, I is the dimension of input space and O denotes the output dimension. The model distribution is denoted by p m (y|x, θ) with model parameters θ. In this paper, labeled dataset where q(y|x l ) is the true distribution of output label given input and p m (.) denotes the model distribution. ǫ is a perturbation added on the input within a small range. σ is a small value. To implement the adversarial training, the inner maximization problem need to be firstly solved to obtain the worst perturbation. For simplifying optimization, the first Taylor expansion of D(.) is used as the approximation for original objective function. With the l p constraint, the worst perturbation can be approximated as: where p * is the dual of p, i.e, 1 p * + 1 p = 1 and ∇L the first derivative of loss function with respect to input x. When p = ∞, this method can be degraded to Fast Gradient Sign Method (FGSM). And the worst perturbation becomes: For finding the worst perturbation, it is easy to compute ∇L with backpropagation. Moreover it does not need a large amount of computation since it just needs one more backward process. After computing the adversarial perturbation, the model distribution p m (y|x, θ) is trained to approximate the true distribution q(y|x) on both natural and adversarial examples. It has been proved to be able to achieve the better generalization performance than the traditional training methods of deep neural network.

B. Virtual Adversarial Training
Another work similar to our method is Virtual Adversarial Training (VAT). In most cases, the full labeled information is not provided and it is important to make use of the information of unlabeled data. Different from the traditional adversarial training methods, Virtual Adversarial Training method can utilize both labeled and unlabeled data. The optimization problem of VAT can be formulated as: where, x is either labeled data or unlabeled data. Different from traditional adversarial training, Virtual Adversarial Training minimize the divergence between two model distributions instead of the divergence between true distribution and model distribution. The objective of VAT is to smooth the output distribution around x through the method of adversarial training. Specifically, inner optimization problem is to find the worst perturbation ǫ within a small range to make the model distributions of x and x + ǫ the most different, oppositely, the outer one is to try to make two model distributions the same. Therefore, VAT does not require the label information and can be applied on unsupervised and semi-supervised tasks.

C. Manifold Adversarial Training
Both the traditional adversarial training and Virtual Adversarial Training methods attempt to find the worst perturbation to make the outputs of models worst. However, these two methods do not consider how the data is represented in latent space. Therefore, we design a new adversarial method Manifold Adversarial Training (MAT) to learn a good geometric structure of the data distribution for improving the robustness of representation in latent space. In following subsections, we first introduce how to regularize the latent space with Maximum Mutual Information (MMI). We then define the notion of smoothness of distributional manifold. After that we describe the major framework and present a series of theory to derive the practical optimization algorithm. Finally, we also briefly conduct a computational analysis.

D. Modeling the Latent Space
In the probability theory, the mutual information is to measure the mutual dependence between two variables. In other words, it measures how much knowing one of these variables reduces uncertainty about the other. For learning discriminative features with more information of labels, MMI is proposed to increase the mutual dependence between the latent features and the class label which can be seen as an information regularization. We define the label set {c j } M j=1 and the latent features where M and N denote the number of classes and samples respectively. D is the dimensionality of the latent space. The mutual information between the event C = c j and the event Z = z i with respect to θ is given by: where C and Z are two random variables representing class label and latent feature respectively. In this work, we assume the conditional probability P θ (z i |c j ), the multivariate Gaussian distribution, as Eq. (2). We also assume that the latent features are generated by Gaussian mixtures. The number of components of Gaussian mixtures is given as the class number in this paper.
In Eq.
(2), µ j is the class center for class c j and Σ j is corresponding covariance matrix. D represents the dimension of latent space. In this work, we assume the class prior P (c j ) is a constant 1/M (M is the class number). Then, the latent feature can be represented by Gaussian mixtures, Again,the number of the components is set to the class number M . Each component of Gaussian mixtures denotes the probability of the latent feature z i assigned to class c j . We reformulate the mutual information I θ (c j , z i ) as follows: In this case, maximizing the mutual information is equivalent to maximizing the log posterior logP (c j |z i ). The last term of Eq. (7) is the log marginal distribution which can be seen as the normalization term. Therefore, we can just maximize the first term logN (z i |µ j , Σ j ). It is equivalent to making the latent features more compact or discriminate with respect to class centers. Then we define the information regularization in this paper: One of our objectives is to maximize the above information regularization term so as to model a good latent space. In this paper, though we assume GMM over the latent space, it presents a generalized version to a recently proposed famous center loss model [13]. This can be readily obtained in Proposition 0.1.
Proposition III.1. The center loss [13] can be viewed as the special case of our proposed information regularization.
Proof. We reformulate the information regularization as: If we make the covariance matrix as identity, we can have: where C is a constant. It is clear that the above last term changes into the center loss as defined in [13].
Remarks. Center loss [13] is also proposed to penalize latent features to be closer to class centers. However, it implicitly assumes an identity covariance matrix as shown in Proposition 0.1. In comparison, we propose a more generalized regularization. With such regularization, we can represent the latent features with Gaussian mixtures which can be easily employed on information geometry and information metric KL-divergence, which we show in the next subsection. Note that, in this work, we just employ such regularization term in the last hidden layer.

E. Defining Low Dimensional Statistical Manifold
An n-dimensional manifold M is defined as a set of points such that each point has n-dimensional extensions in its neighborhood and such a neighborhood is topologically equivalent to an n-dimensional Euclidean space. Intuitively speaking, it can be viewed as a deformed Euclidean space.
In this paper, the dimension of the entire manifold S 0 of latent space is denoted as D and the coordinate system is defined as ζ = (ζ 1 , ζ 2 , ..., ζ D ). In the previous section, the feature points are modeled by Gaussian mixtures. Hence, the set of all the Gaussian mixtures is an M -dimensional statistical submanifold S 1 , where a point on manifold S 1 denotes a Gaussian mixture function and u = (u 1 , u 2 , ..., u M ) describes a coordinate system (u i is the i th component of Gaussian mixture) as illustrated in Figure 1. The degree of separation of two points on manifold S 1 is measured by the KL-divergence between two Gaussian mixtures. In this paper, our aim is to smooth the manifold S 1 . Fig. 1. The region of circle represents the entire manifold S 0 of latent space defined by the coordinate system ζ = (ζ 1 , ζ 2 , ..., ζ D ) and the grey region in circle denotes a submanifold S 1 defined by the coordinate system u = (u 1 , u 2 , ..., u M ). The degree of separation between two points p 1 and p 2 is described by KL-divergence.

F. Defining Statistical Manifold Smoothing
In this section, we consider how to evaluate the smoothness of statistical manifold. First, we define some notations. The input set is denoted as {x i } N i=1 and the output of model is defined as The training set is given as follows: We then use D to train the model distribution P θ (y|x, θ) with the regularization term I θ . In the previous subsection, the latent features are represented by the Gaussian mixtures u(z i ). We can then easily define the distance between the two points on the statistical manifold with KL-divergence KL[p||q] between mixtures p and q. For learning a good manifold in the latent space, we try to preserve the local information in the embedding as in [11] [12]. Specifically, we add the small adversarial perturbation to input samples and try to make close the latent representation of perturbed samples and original one. Since the latent features are represented by Gaussian mixtures in this work, the information metric KL-divergence can be readily used to measure the similarity. The notion of smoothness of statistical manifold can then be defined as the variation of the latent features on the statistical space caused by the adversarial perturbation ǫ as follows: where f (x, θ) denotes the latent feature of x with the model parameters θ. P GM (f (x)) = 1 M M j=1 N (f (x)|µ j , Σ j ) which is defined as the Gaussian mixture representation for latent features f (x). ǫ represents the adversarial perturbation.
The smaller value of M KL (x, θ) means the statistical manifold is more smooth at x even when the perturbation is imposed. We also define the smoothness in the output or label space as in [6] [8]: In more details, L KL (x, θ) describes the KL-divergence between the output space of adversarial example and the original examples. The overall smoothness for both the output and manifold in latent space are then given as follows:

G. Final Optimization Problem of MAT
In the previous subsections, we have defined the information regularization, manifold smoothness, and label space smoothness. We can now obtain the final optimization problem of our proposed framework as follows: where the hyperparameters λ 1 > 0, λ 2 > 0 and σ > 0.
In the above objective function, the first term represents the data log likelihood, the second term defines the overall data smoothness in both the manifold and output space, and the last term presents the information regularization. On one hand, the proposed MAT framework tries to optimize the model parameter θ so as to find the best latent space (by maximizing the information regularization I θ ), enlarge the data log likelihood (so as to fit the data) as well as increasing the data smoothness (decreasing N i=1 logS kl (x i , θ)); on the other hand, the imperceptible perturbation ǫ tries to minimize the above objective function. In other words, the proposed novel model tries to find the best parameter θ that is even robust to the worst perturbation as given by the optimal ǫ. Since this adversarial learning is defined in the latent manifold space, we call this model as Manifold Adversarial Training (MAT). In the next section, we will use a strategy similar to [5] and [6] and discuss how to solve this optimization problem practically and efficiently.

H. Practical Algorithm
Similar to [5], we try to solve the above optimization problem in an alternative way. We first solve the inner minimization problem with respect to ǫ, i.e. the worst perturbation, which we denote as ǫ adv . In order to calculate this worst perturbation, the KL-divergence between GMMs needs to be calculated firstly. Since it is difficult to conduct the evaluation directly, we approximate it by matching between the Gaussian elements of the two Gaussian mixture density as described in [21]: where Ψ ij is a K × K stochastic matrix. For simplifying the calculation, we try to optimize the upper bound of KLdivergence with assuming Ψ ij an identity matrix and the mixture weight α = 1/K and β = 1/K. Since it is difficult to evaluate the worst perturbation in a non-convex problem, we relax it to a convex problem with second-order Taylor expansions as in [6] [8]. Since M KL (ǫ, x, θ) reaches the minimum value when ǫ = 0 and P (y|x, θ) is differentiable with respect to x and θ, the first derivative ∇ ǫ M KL (ǫ, x, θ) = 0. We can finally approximate it as follows: where, H(x, θ) is the Hessian matrix calculated by employing the second derivative, H(x, θ) ≡ ∇∇ ǫ M KL (ǫ, x, θ)| ǫ=0 . Before we could obtain the worst perturbation, we first present Lemma III.2 and Lemma III.3 as follows: Lemma III.2. Let B be a real square matrix in R × R and u be the dominant eigenvector of B and e is a vector which is not perpendicular to u. Then, the iterative calculation of e ← Be (15) will make e converge to u. ( . represents the normalization operator) Lemma III.3. Let A be a Hessian matrix function F (x) with respect to x in R × R and r be a vector in R. Then we have: Proof. We give the Taylor expansion of the first derivative of function F (x + ar): then we have: By using Lemma 0.2 and Lemma 0.3 , we could solve the inner minimization problem, i.e., obtain the current worst perturbation. Specifically, we have Theorem 0.4 showing that it can then be written as the production of most sensitive directiond and the scale ǫ, where the most sensitive direction d can be approximated iteratively by the power method. (19) where, . represents the normalization operator and ξ > 0. d is a vector which is not perpendicular to dominant eigenvector of H(x, θ).
Theorem III.4. The steepest directiond in (14) can be approximated by iterative calculation of: where . represents the normalization operator.
Proof. Applying the Lagrange multiplier method on (14), we can easily get that the steepest directiond has the same direction with the dominant eigenvector of H(x, θ). Then we can easily approximate the steepest directiond by iteritive calculation of (refer to Lemma 0.2): where, d is a vector which is not perpendicular to dominant eigenvector of H(x, θ). Then we can calculate Bv using Lemma 0.3: where, ξ is a small value. Then we have: After we solve the inner minimization problem with respect to ǫ, we can then solve the maximization problem with respect to θ. We iterate these two steps until the process converges. The detailed pseudo code algorithm is shown in Algorithm 1.

I. Computational Analysis
We briefly touch the computational analysis here. Compared with previous methods like VAT [6], our proposed method needs to compute the mean and covariance matrix for Gaussian mixtures additionally. In this work, parameters of Gaussian mixtures are evaluated by the latent features of labeled data instead of stochastic gradient. It may be inefficient to compute the inverse of covariance matrix. However, in the lowdimensional latent space, the inverse of covariance matrix can be easily and directly obtained; even in the very highdimensional space, the covariance matrix can be assumed as a diagonal matrix which is also easy to get its inverse. Same as VAT, we just need one iteration of the power method to estimate the worst perturbation ǫ adv . The whole procedure for updating the parameters of deep neural networks consists of Update the parameters of neural network with stochastic gradient: Update the covariance matrix and mean of Gaussian Mixtures using the latent features: 11: for i in 1...m do 12: end for 15: end for two forward and two back propagations. The first forward and backward propagations are used to evaluate the worst perturbation. After calculating the final loss, propagate forward and backward again to update the parameters of neural network and Gaussian mixtures.

IV. EXPERIMENT
For assessing the efficacy of our proposed MAT, we implemented it on several benchmark datasets including MNIST, CIFAR-10, and SVHN. In principle, adversarial methods can be regarded as a robust method which typically have better generalization abilities than the traditional methods. To check if the proposed MAT can indeed improve the classification performance, we first applied MAT on the on benign data of MNIST and CIFAR-10 in both the supervised and semisupervised task. Additionally, to check further the robustness of our proposed adversarial training framework, we perform various experiments to examine how the proposed MAT could defend the attack from various adversarial methods. To illustrate our proposed method, a series of visualizations were made, offering some interesting results which might be used to explain the adversarial examples.

A. Experiments on Benign Data
We first evaluated the performance of our proposed MAT methods against many other comparison methods on the benign data. We first report the performance of various methods in the setting of supervised classification and then perform the comparison in semi-supervised learning.
1) Supervised Learning: We implement the same framework Lenet++ with [13] on MNIST dataset. For Lenet++, there are only two dimensions in the last hidden layer which is convenient for visualization. The based framework for experiment on CIFAR-10 is the same as [6] called Conv-Large. For MNIST dataset, we train the deep model with 60, 000 labeled training samples, and we evaluate it with 10, 000 test samples. For CIFAR-10, we use 50, 000 training samples and 10, 000 test samples. To search good hyper parameters, the training set is divided into 50, 000 training set and 10, 000 validation set. We choose a set of hyper parameters with the best performance on the validation set. For MNIST dataset, the best hyper parameters are obtained as: λ 1 = 0.1, λ 2 = 0.1 and σ = 2. For CIFAR-10, the parameters are obtained as: λ 1 = 0.01, λ 2 = 0.1 and σ = 20.  Table I and Table II list the performance of our proposed method and other competitive methods on MNIST and CIFAR-10. Generally, adversarial methods can be regarded as a robust method which typically have better generalization abilities than the traditional methods. As observed from both the tables, MAT clearly achieves the best performance among all the comparison methods in both the datasets. It is probably unfair to compare with some methods, e.g., ladder network, simply  II  TEST PERFORMANCE ON CIFAR-10 IN SUPERVISED LEARNING   Method  CIFAR-10 Test error rate(%) Network in Network [22] 8.81 All-CNN [23] 7.25 Deeply Supervised Net [24] 7.97 Highway Network [25] 7.72 RPT [6] 6.25 ± 0.04 Baseline 6.76 ± 0.07 VAT 5.81 ± 0.02 MAT 4.40 ± 0.03 because the different based model were used there. However, it is sufficient to show the superiority of our proposed MAT due to its significant lower error rate than most of the other competitive adversarial methods.
2) Semi-supervised Learning: In the methodology part, we have introduced the MAT's final objective function Eq. (13) where the first term is the ordinary soft-max loss function and the last term is the information regularization term. These two terms are exploited on the training data with labels. And the second term is used to smooth the manifold and output distribution which does not need the label information. Therefore, our method can readily be extended to semi-supervised learning. Following the same setting as [6], we implement our MAT both on CIFAR-10 and SVHN with 4, 000 labeled data and 1, 000 labeled data respectively. We use the same base model with [6] called Cov-Large with batch normalization and dropout. We exploit the mini-batch of size 32 for both labeled data and unlabeled data on CIFAR-10. For SVHN, we use the labeled batch with size 32 and unlabeled batch with size 128.

RPT [6]
18.56 ± 0.29 VAT 14.82 ± 0.13 MAT 12.85 ± 0.21 Table III and Table IV show the performance of our proposed model as well as recent state-of-the-art semi-supervised learning methods on CIFAR-10 and SVHN. Once again, our proposed method MAT achieves the highest performance, which further shows its advantages.
3) Visualization: In order to illustrate why the proposed MAT could perform excellent, we take the MNIST as one example to visualize various methods. In particular, Figure 2 shows the illustration of the last embedding space of Lenet++ with center loss, softmax, VAT, and MAT on the MNIST test data. 1 It is obvious that the latent features of our method are represented more compactly and discriminatively with respect to class centers. VAT learns the similar latent space with softmax, but apparently both of their features are not as discriminative as our proposed MAT. This may illustrate why our proposed MAT method could usually generate better classification performance than VAT and the traditional CNN.
Since our proposed MAT seeks a smooth latent space which we believe would benefit the classification, we plot in Figure 3 the learning curve (accuracy rate) and the smoothness for the three different methods: MAT (our proposed method), VAT, and the baseline (traditional CNN) on CIFAR-10. 2 For MAT, the best set of hyperparameters, i.e., λ 1 = 0.01 λ 2 = 0.1 and σ = 20 was used, while for VAT, the best setting reported in [6] was directly applied. Compared with the other two methods, our proposed MAT clearly increases the smoothness R smooth , which is defined by the average of the smoothness of output space and latent space: When we calculate the smoothness, the perturbation needs to be normalized to unit vector. Figure 3(b) shows that our proposed method indeed learns the smoother latent and output space.

B. Performance on Defending Adversarial Examples
We now turn to examining how the proposed MAT method could defend the attacks from various adversarial examples generation approaches in comparison with other competitive methods. Visualization is also presented so as to obtain further understandings on adversarial examples.
1) Robustness to Adversarial Attacks: We implement different methods on FGSM and 2-norm adversarial attack for MNIST, CIFAR-10, and SVHN datasets. Particularly, we generate in the test sets of MNIST and CIFAR-10 10,000 adversarial examples according to FSGM and 2-norm attacks [5] respectively. For SVHN, we generate 26,032 adversarial examples. We increase the level of adversarial noise gradually from 0 to 8 in MNIST and SVHN with the step size as 1 and from 0 to 13 in CIFAR-10 with the step size as 1.6. We then test the performance of various training methods on these adversarial examples. The performance is plotted in Fig. 4. As clearly observed, the proposed MAT shows better robustness against the two types of adversarial examples. Particularly, when the adversarial noises are small, all the adversarial training methods show similar results but perform much better than the CNN (exploiting no adversarial training); when the adversarial attacks are stronger, the proposed method overall demonstrates clearly better performance in almost all the cases, verifying its significant robustness. One exception can be identified in Fig. 4(c) where the FGSM training method performs the best while our proposed MAT performs the second. However, since the FGSM training was specifically designed to defend the FGSM adversarial examples, it will not be very surprising that FGSM can do better against the FGSM adversarial examples. Nonetheless, except Fig. 4(c), our MAT still demonstrates the best robustness in all the other cases.
Again, taking MNIST as one example, we show the TSNE embedding for the second last layer given by VAT and MAT under various levels of 2-norm adversarial attack. As Figure 5 shows, when the degree of adversarial noise increases (from 0 to 6 with step size 2), the clusters of different classes are less discriminative. In contrast to VAT, our proposed MAT obtains much more discriminative features on different levels of adversarial perturbation.
2) Visualization of Adversarial Examples: In this subsection, we visualize some adversarial examples generated by both our proposed method and the other comparison adversar-ial methods. To visualize the adversarial examples, we directly train the deep neural network with original images and then calculate the worst perturbation through back propagation. For better visualization, the image up-sampling was applied.
Specifically, we compare the adversarial examples generated by our method (MAT), VAT, and traditional adversarial learning on MNIST, SVHN, and CIFAR10. Figure 6-8 show several adversarial examples (generated by different methods) from the three datasets. It is interesting to note that after the adversarial training given by MAT, the original images tend to be morphed into some other categories. Particularly, digit "3" was changed into "8", "7" was changed to "9", and "9" was changed to "0" on MNIST (see Figure 6). Similar cases can also be observed in Figure 7 where the worst perturbation generated by MAT tends to change the number "6" to number "8" (the 1st image in the 2nd row of Figure ??), the number "3" to "5" (the 4th image in the 2nd row of Figure 7), and numbers "66" to "00" (the 8th image in the 2nd row of Figure 7). In comparison, the other adversarial methods also appear to generate adversarial examples in a same sense but not as obvious as MAT. This may partly explain why the adversarial examples would be difficult to be recognized by many traditional deep neural networks. Furthermore, generating the most challenging examples (being similar to or even changed to other categories), MAT indeed leads to the most serious adversarial attacks that deep neural networks hardly recognize. Finally, it is observed that the major changed areas occur on the number strokes, which can be considered as the "manifold" of these digit images. This implies that our MAT training method could truly generate adversarial attacks in the manifold.
Last, we also show in Figure 8 adversarial examples on CIFAR-10 generated from various methods. Unfortunately, very different from MNIST and SVHN, it appears that the resulting examples give no obvious clues on how the images would be changed. This might be explained by the fact that general objects are sufficiently more complicated than digits. We would leave the topic of generating explainable adversarial examples in the future.
V. CONCLUSION We present the Manifold Adversarial Training (MAT), a novel method to smooth the distributional manifold in the latent space. Compared with other adversarial training methods, our proposed MAT learns adversarially a robust feature representation in the latent space, making the latent space both informative and discriminative. Specifically, we first represent the latent features with Gaussian mixtures. We then define the smoothness of distributional manifold based on the KL-divergence between Gaussian mixtures for the original and adversarial examples. We implemented MAT on MNIST, CIFAR-10 and SVHN in both supervised and unsupervised tasks. The results showed that our proposed model is much better than those state-of-the-art methods.


Convolutional neural network compression for natural language processing

Convolutional neural networks are modern models that are very efficient in many classification tasks. They were originally created for image processing purposes. Then some trials were performed to use them in different domains like natural language processing. The artificial intelligence systems (like humanoid robots) are very often based on embedded systems with constraints on memory, power consumption etc. Therefore convolutional neural network because of its memory capacity should be reduced to be mapped to given hardware. In this paper, results are presented of compressing the efficient convolutional neural networks for sentiment analysis. The main steps are quantization and pruning processes. The method responsible for mapping compressed network to FPGA and results of this implementation are presented. The described simulations showed that 5-bit width is enough to have no drop in accuracy from floating point version of the network. Additionally, significant memory footprint reduction was achieved (from 85% up to 93%).

Abstract-Convolutional neural networks are modern models that are very efficient in many classification tasks. They were originally created for image processing purposes. Then some trials were performed to use them in different domains like natural language processing. The artificial intelligence systems (like humanoid robots) are very often based on embedded systems with constraints on memory, power consumption etc. Therefore convolutional neural network because of its memory capacity should be reduced to be mapped to given hardware. In this paper, results are presented of compressing the efficient convolutional neural networks for sentiment analysis. The main steps are quantization and pruning processes. The method responsible for mapping compressed network to FPGA and results of this implementation are presented. The described simulations showed that 5-bit width is enough to have no drop in accuracy from floating point version of the network. Additionally, significant memory footprint reduction was achieved (from 85% up to 93%).

I. INTRODUCTION
N ATURAL Language processing (NLP) is considered one of three main pillars of Deep Learning along with image and video processing. Mobile devices are becoming increasingly dominant both in terms of their number as well as a computing load and network traffic. They also become more interactive in terms of NLP algorithms implemented for applications such as translation or voice-typing. Consequently, there is a significant benefit from reducing memory footprint and computational burden for deployment of NLP modules on embedded devices. It is worth noting that despite an abundance of research efforts in Deep Learning architectures compression for image processing [4], [5], [11], there are just a few projects aiming at NLP neural architectures compression [1], [2], [3], [12].
This paper shows a case study of employing common compression techniques for compressing NLP architectures. A series of quantization and pruning techniques were implemented and deployed in FPGA platform. The main goal of the paper was to examine feasibility and efficiency of using FPGAs in a domain of embedded neural computations. A set of common datasets were used to conduct reliable experiments. The experiments reviled a strong relationship between a size and the structure of datasets and performance of quantization and pruning methods used. All the layers of the neural architecture [6] were analyzed separately with respect to their ability to be quantized and pruned.
The structure of the paper is structured as follows. Section II provides an overview of CNN architectures, quantization and pruning processes. Section III covers the dataset used for the experiments. Section IV describes neural model used for the experiments as well as FPGA implementation details of the architecture. Finally, section VI contains the results of the experiments. Section VII summarizes contribution of this work.

II. CONVOLUTIONAL NEURAL NETWORK
CNNs are composed of neurons that have learnable weights and biases. Each neuron receives inputs structured as multidimensional vectors (tensor), perform a convolution operation on this input it with multi-dimensional filters and then optionally follows outputs with pooling and non-linearity functions. Typically, the layers of a CNN have neurons that generate outputs feature maps y i , i = 1...N as follows where F ij are two-dimensional (2D) convolutional kernels of dimensions H × W , * represents the convolution operation and b i are the bias terms. The number of multiply-accumulate (MAC) operations and cycles spent on the execution of in a practical implementation is often used as the metric for a complexity of a CNN. Assuming each output feature map y i has P elements (where P is equal to height multiplied by width of given feature map) the total number of MAC calculations for a convolutional filtering operation is M ACs = P HW M N . Our role in quantization is to reduce the complexity of each one of these operations and in network compression through pruning the goal is to reduce the total number of operations.

A. Quantization process
Quantization is the procedure of constraining values from a continuous set or more dense domain to a relatively discrete set (e.g. integers), more sparse domain in which quantized input values will be represented. In our case, the domain is floating point representation. A floating-point number can be represented as m ∈ Z is the mantissa, b = 2 is the base, and e ∈ Z is the exponent. In case of single precision floating point format according to the IEEE-754 standard the mantissa is assigned 23bits, the exponent is assigned 8bits, and 1bit is assigned for Similarly, a reduced precision IEEE-754 mini-float format assigns 10bits of mantissa, 5bits of exponent and a sign bit. Currently, GPUs with parallel processing being applied to machine learning operate mainly using the single precision format. In contrast, embedded DSPs and the latest GPUs operate with 8-bit integer and 16-bit fixed-point processing that restricts numbers to a range fxp : 2 −frac_bits ·{−2 total_bits−1 , . . . , 2 total_bits−1 −1} (4) where total_bits = 8, 16 are conventional bit-widths and frac_bits is a shift down (or up) that determines fractional length, or number of fractional bits, as well as the integer length int_bits = total_bits − frac_bits − 1 for signed numerical representation. When the fractional length varies over individual coefficients or data samples. This format is also referred to as dynamic fixed-point [10].
It is possible to define a general mapping from a set of floating-point data x ∈ S to fixed-point q as follows (assuming signed representation) In our case µ = 0 and σ = 2 −frac_bits where and frac_bits = total_bits − int_bits − 1. The scaling factor σ is essentially just a shift up or down. A drawback is that a great deal of precision may be lost if the distribution of the dataset S is skewed by a large mean. Yet another approach can define the number of integer and fractional bits to represent regions of a distribution that will represent a large percentage of the range. In these cases, there will be saturation of a small percentage of the data, such as outliers, through the quantization procedure which may or may not affect the accuracy in a significant way. To determine the effects of saturation on can experiment with different saturation levels. Therefore histogram analysis is used to analyze outliers and set best levels of saturation.
Another approach is quantization that maps floating point values to integers: The µ parameter can be set to min(X) (X is a input set of values to be quantized) or can be zero value. In a first case it is known as a asymmetric integer quantization (e.g. used in Tensorflow framework), in the second it is called symmetric. The compression system presented in the paper are based on the symmetric quantization and dynamic fixed point was also implemented.
All presented approaches are examples of linear quantization. It is possible also to use nonlinear version to minimize quantization loss but its hardware implementation is more sophisticated and it is more difficult to achieve gain in used hardware resources.

B. Pruning
After training neural model we acquire a set of weights for each trainable layer. Those weights are not evenly distributed over the range of possible values for a selected data format. As presented in figure 1, most weights are concentrated around 0 and are very close to it. Therefore, their impact on resulting activation value is not significant. The mechanism of pruning removes weights which values are below a certain threshold level. Authors of [9] examined how pruning affects convolutional neural networks for image classification, showing little accuracy drop and high compression ratio. In this work pruning is applied to 1D convolution layer according to algorithm 1, as weights are concentrated around 0, for negative values threshold was calculated for an absolute value. Example from figure 1 presents that threshold level of 0.02 is sufficient to remove around 50% of weights. Depending on specific network implementation storing weights may require significant amounts of memory, removing weights through pruning have a direct impact on lowering storage requirements.
The order in which quantization and pruning are applied should have little effect, especially for higher precisions. If quantization is applied first then weights from a certain range will be brought to the same value and pruning will only be able to remove none or all of them. If pruning will be applied first it will be possible to cut out weights more precisely. This will only have an effect on border quantization buckets and therefore will have any significant impact only for very small precisions.

III. DATASETS
Experiments are performed on the same datasets as used in [6]. Summary statistics of the datasets are provided in table I. Some datasets are divided into training, validation and testing data. If validation data is not specified then random 10% of training is used for it. If testing data is not prepared then 10fold cross validation is performed.

IV. NETWORK ARCHITECTURE
We employed a model CNN-static from [6]. The model uses pre-trained 300-dimensional word embeddings GloVe [7], which are not fine-tuned during training. Unknown words are randomly initialized. Outputs of two parallel CNN layers with The concatenation is connected to dense layer with 128 outputs. CNNs and dense layer have rectified linear units as activation functions. The last layer uses softmax activation function in case of many outputs and sigmoid activation function in case of binary classification. After each layer dropout was applied with probability 50%. Training data was divided into mini-batches of size 128. The training involves Nadam [8] as optimizing algorithm. Training lasts 25 epochs, the best model is saved.

A. FPGA Implementation
Pruning mechanism is especially effective in FPGA implementations. On standard CPUs processing throughput improvement is not significant, multiply and accumulate operation is replaced by conditional, and on most modern processors both of those operations are supported by hardware. Greater improvements can be achieved when pruning is extensive enough to apply algorithms for sparse operations, such extreme pruning, however, will likely has a considerable impact on accuracy. When creating FPGA implementation of a certain neural network there is a higher degree of freedom. Its possible to hard code weights into logic during synthesis, as a result a complete cost of resolving condition if certain calculation should be performed, is moved from inference to compilation. Weights that will be removed by pruning are not included in the synthesized design, saving both logic and routing resources.
FPGA implementation can also give additional freedom for quantization adjustments, as it is not bound to any standard data width or format. On a standard CPU moving from 16 bit precision to 15 bit will not give any performance improvements, a calculation will be executed using the same hardware. For FPGAs however, it will result in synthesizing a smaller logic design with fewer connections, which can be beneficial for latency, throughput or calculations density.
Weights in the used neural network are mostly concentrated in Convolution operations, around 99% of total trained weights, with 1% in fully connected layers, as embeddings  layer is not implemented in FPGA it is excluded from these calculations. Therefore FPGA implementation is concentrated on reducing and accelerating convolutional layer. The design was implemented and synthesized using Xilinx Vivado HLS environment according to algorithm 1. To allow for more effective pruning, all weights are hard coded into logic. Loops from lines 2, 3 and 4 are fully unrolled. Due to limitations of synthesis tools layer from the used neural network could not be used. Therefore experimental results presented in section VI-A are obtained from a smaller version of a similar network, the results, however, can be extrapolated. The process of transforming Keras model to an IP core that can be used in FPGA is presented in the figure 2.
In the first step, Keras network model is created and trained, when accuracy is satisfactory network model and weights need to be exported. Weights from an HDF5 file are processed by a python script to C++ header format. It is important that their values will be known during compilation as Vivado HLS will not be able to create hard coded network otherwise. A separate header is created for each layer. Json file is used to create infer function of a neural network. This is a top level function that performs all the computations on input data. At this, we can decide how the data and weights should be represented. To validate the generation process a simulation was performed. If results are satisfactory, Vivado HLS will synthesize networks IP core. Next, the IP core is put into Vivado design and bitstream is generated.

V. PRECISION REDUCTION
The precision reduction was applied to two types of places: weights of layers and outputs of activation functions. The process is performed by finding minimum and maximum of values to be reduced and uniformly dividing value range into 2 bits buckets. Each bucket has assigned a middle point of represented subrange. In case of layer weights, it is straightforward. Precision reduction of activation functions requires to run computation on training data and remember minimum and maximum of activation functions for chosen layers. During testing, it is possible that computed values will be lower than minimum and larger than maximum -they need to be clipped.
The model used in this work has 5 layers: embeddings, two CNNs, and two dense layers. Activation functions are reduced after the two CNNs and the first dense layer. A number of bits used in each of the 8 places may be different.

VI. EXPERIMENTS
The first set of experiments involved precision reduction to the same number of bits on every combination of the 8 places. Tested number of bits includes 32, 16, 8, 7, 6, 5, 4, 3, 2 and 1. It gives 2 8 ·10 experiments for each dataset/fold. All presented results are based on integer quantization. In all experiments this type of quantization gives little bit better accuracy than in dynamic fixed point (difference up to 2%).
The results show that precision reduction to 16 bits does not change accuracy. Embeddings can be reduced to 5 bits, a majority of datasets weights of convolutional layers can be reduced to 4 bits. It is important to mark that convolutional layers are parallel and reducing one of the layers to 1 bit does not change overall accuracy significantly because the second is operating with full precision.
Tables II and III presents results for SST-2 and MR datasets for each place precision reduced separately. In the case of 32 bits there is no reduction as it is default precision. 16 bits do not change accuracy results. Accuracy increases for some settings. The maximal drop in accuracy less than 1% is observed down to 4-5 bits. Usually, quantization of activa-tions decreases accuracy more than quantizations of weights. Precision reduction to less than 4 bits of activations of first dense layer affects negatively accuracy the most. The tables present also results for precision reduction of all places. a joint reduction can be better than only for one place.
The second experiment was focused on finding the strongest precision reduction with arbitrary chosen maximal accuracy decrease. Here we are not limiting the model reduction to the same number of bits for each place. In order to test every combination including 10 values of a number of bits, it would be necessary to perform 10 8 tests, which is not computable in a reasonable time. To tackle this problem a random-restart hill climbing algorithm was employed. The algorithm starts with a randomly chosen order of places to be reduced. Then iterates the places in a cyclic way until does not increase model compression. For each layer, it tries to find the smallest number of bits that does not decrease accuracy more than a set threshold. The searching method using algorithm 2 was restarted 50 times. In this work, we set maximal decrease as 0.2% of accuracy without precision reduction. Table IV shows results of precision reduction for each  dataset. In comparison to 32 bit representation, the size of models can be reduced by more than 84%. However, the most space is taken by word embeddings (more than 93%). It could be further reduced by limiting vocabulary and reducing dimensions of embeddings. Another approach was to order the places by their influence on model size, in this case starting from optimizing embedding layer. The approach has not achieved better results than the algorithm 2. The same precision bits parameters were obtained for two datasets.
The best found precision settings are shown in table V. All precision values are equal to or less than 8 bits.
The third experiment involved a combination of quantization and pruning. It was performed on MR database, the number of bits changed in the same way as in previous experiments. Prun-ing was applied according to algorithm 1. pruning threshold started from 0, i.e. no pruning, with 0.005 step to 0.15, where all weights were removed and no multiply and accumulate operations were performed in convolutional layers. Figure 3 presents extensiveness of pruning for each data precision. It is visible that the ratio quickly drops, it is a result of Gaussian distribution of weights concentrated around 0. When quantization steps are greater than pruning threshold step ratio decreases more gradually.
Pruning impact on the accuracy of the network is presented in figure 4. Its effect on accuracy is highly separate from quantization. A significant drop in accuracy appears below 5 bit precision, and over 0.03-0.035 pruning threshold. Between these values and no reduction point, a rectangular plateau of only slight variations in precision can be observed. Described quantization and pruning mechanisms were implemented in FPGA to examine their impact on logic utilization. Experiments were performed on a smaller convolutional layer, with an input of size 64 and 35 feature maps, the filter of size 2 and output of size 63 and 16 feature maps. We checked what effects gives pruning and quantization separately and in conjunction. Experiments were performed on Xilinx Virtex-7 FPGA VC707 Evaluation Kit with Virtex 7 series FPGA. Test points were selected so the accuracy remained within 1% of the original. Results are summarized in table VI. When both techniques are used flip flips utilization drops 3.3 times and look up tables utilization drops 4.2 times, while having a little effect on the accuracy.

VII. CONCLUSIONS AND FUTURE WORK
This paper considers several compression and pruning techniques for convolutional neural networks in sentiment analysis. The results of the experiments show that the size of the model may be reduced more than 84% with the little performance degradation of approximately 1%. Word embeddings can be compressed even more reaching 93%. Future work will concentrate on boosting the performance of given convolutional network by changing the architecture using deeper model and modifying filter scheme. The next issue will focus on using memetic algorithm and reinforcement techniques to find a configuration of the network which gives the best compression.