Sparse evolutionary Deep Learning with over one million artificial neurons on commodity hardware

Microarray gene expression has widely attracted the eyes of the public as an efficient tool for cancer diagnosis and classification. However, the very-high dimensionality and the small number of samples make it difficult for traditional machine learning algorithms to address this problem due to the high amount of computations required and overfitting. So far, the existing approaches of processing microarray datasets are still far from satisfactory and they employ two phases, feature selection (or extraction) followed by a machine learning algorithm. In this paper, we show that MultiLayer Perceptrons (MLPs) with adaptive sparse connectivity can directly handle this problem without features selection. Tested on four datasets, our novel results demonstrate that deep learning methods can be applied directly also to high dimensional non-grid like data, while learning from a small amount of labeled examples with imbalanced classes and achieving better accuracy than the traditional two phases approach. Moreover, we have been able to create sparse MLP models with over one million neurons and to train them on a typical laptop without GPU. This is with two orders of magnitude more than the largest MLPs which can run currently on commodity hardware.


Introduction
In the past decades, data have become the indispensable factors of scientific progress, medical development and economic growth. Without the increase in the number of avail-able data, scientific development cannot have such an incredible speed. Especially, gene expression obtained from DNA microarray has emerged as a powerful solution to cure cancers (Simon et al., 2003). However, most of the datasets in DNA microarray are high-dimensional and redundant that would result in unnecessary calculations, huge memory requirements and even the decrease of generalization ability due to the "curse of dimensionality" (Destrero et al., 2009). Moreover, the invisible relationships and nonstandard structures among different features also make it very time-consuming to find the key features from tens of thousands of features. In order to tackle this problem, various methods have been proposed by researchers. Among them, feature selection is undoubtedly a "de facto" standard as it is not only able to remove the redundant features but to improve the classification performance (Destrero et al., 2009). Following the feature detection phase, standard classifiers can be used to perform classification based on the selected features. One of the most used classifier is Multi-Layer Perceptron (MLP), e.g it represents 61% of a typical Google TPU (Tensor Processing Unit) workload for production neural networks applications, while convolutional neural networks represent just 5% (Jouppi et al., 2017).
Motivation. Yet, MLP can not be employed directly on high dimensional data due to the quadratic number of parameters in its fully connected layers with respect to their number of neurons. This limits MLPs size to several thousands neurons and few thousands input features on commodity hardware, and implicitly their representational power. To address this issue, very recently small steps have been made. (Mocanu et al., 2018) proposes Sparse Evolutionary Training (SET), a method to train scalable MLPs with adaptive sparse connectivity (SET-MLPs). However, due to the limitations of typical deep learning libraries (e.g. optimized operations just for fully-connected layers and dense matrices), in (Mocanu et al., 2018) the largest SET-MLP used has just 12,082 neurons -quite a low representational power. Practically, their SET-MLP implementation uses the typical approach from the literature to work with sparsely connected layers, i.e. fully connected layers with sparsity enforced by a binary mask over their weights -this approach, of course, is far from using the full advantage of sparsity. arXiv:1901.09181v1 [cs.NE] 26 Jan 2019 The first contribution of this paper is an efficient implementation framework which can create and train SET-MLP models with over one million neurons on a typical laptop to handle data with tens of thousands of dimensions. This very high representational power is way over the capacity of state-of-the-art SET-MLPs and fully-connected MLPs. Secondly, we show that our proposed approach can be a good replacement for the current methods which employ both, feature reduction and classifiers, to perform classification on high-dimensional non-image datasets such as microarray gene expression data with imbalanced classes. Thirdly, we show that our proposed solution is robust to the "curse of dimensionality", avoiding overfitting and achieving very good performance in terms of classification accuracy on a dataset with over 20,000 dimensions (input features) and less than 100 samples.
The remaining of this paper is organized as follows. Section 2 introduces and discusses our proposed methods. Section 3 presents the experiments performed and analyses the results. Section 4 discusses two extreme SET-MLP models and shows that a SET-MLP model with over one million neurons can be trained on one CPU thread of a typical laptop. Section 5 discusses related work, while Section 6 concludes the paper and presents future research directions.

Methods
This section introduces our proposed methods. First, it discusses the Sparse Evolutionary Training procedure and its current limitations given by the state-of-the-art deep learning libraries and techniques. Secondly, it describes our novel proposed solution to address those limitations.

Sparse evolutionary training
Inspired by the fact that biological neural networks are prone to be sparse, rather than dense (Strogatz, 2001;Pessoa, 2014), and due to obvious computational resources limitations, there is an increasing interest in conceiving neural networks with a sparse topology (Mocanu et al., 2016). In (Mocanu et al., 2018), the authors proposed a novel concept, Artificial Neural Networks (ANNs) with adaptive sparse connectivity. The basic idea is to replace the fully connected layers with sparsely connected layers before training in any type of neural network, and after that during the training process to optimize together the weights values and the network sparse topology to fit the data distribution. They also proposed a scalable training method, i.e. Sparse Evolutionary Training, to train such networks. Different from the conventional methods, e.g. weights pruning (Cun et al., 1990;Han et al., 2015) which creates sparse topologies during or after the training process, the adaptive topology of ANNs trained with SET is designed to be sparse from the begin. This reduces quadratically the amount of connections.
SET algorithm is given in Appendix A, Algorithm 1. For the sake of convenience, we briefly describe the structure of SET, using the same notations with the original paper. The original sparse topology is initialized by Erdős-Rényi random graph topology (Erdős and Rényi, 1959) where a sparse matrix W k ∈ R n k−1 ×n k is defined to determine the probability of each connection between two consecutive layers of neurons(i.e. h k−1 and h k ). The weight matrix is given by whereby the n k and n k−1 represent the numbers of neurons of hidden layers h k and h k−1 , respectively. ∈ R + is a hyperparameter to control the sparsity level.
Like conventional fully-connected MLP (FC-MLP), SET-MLP also employs backpropagation with stochastic gradient descent to learn the best weights for different datasets. However, the initial sparse network may not be suitable for every type of dataset, since it is generated randomly with no information about the data distribution. To overcome this problem, in each epoch, weights pruning based on magnitude is used to eliminate the non-informative connections. More exactly, a certain fraction ζ of the largest negative weights and the smallest positive weights are removed. After weights removal, an equal number of connections are randomly added to the bipartite layers. Roughly speaking, the connections removal in SET represents natural selection, whereas the emergence of new connections corresponds to the mutation phase.
Very recently, several papers proposed various techniques for adaptive sparse connectivity. (Bellec et al., 2018) proposed DEEP R. (Mostafa and Wang, 2019) proposed dynamic sparse reparameterization for deep convolutional neural networks and showed that their proposed method and SET are faster and achieve better accuracy than DEEP R. (Zhu and Jin, 2018) proposed a simplified SET variant (after weights removal, they do not add random connections) to minimize the size of the network as much as possible as needed in federated learning and low-resource devices. Considering the above, in this paper, we focus on the original SET algorithm because it was shown that it is capable to reach very high accuracy performance (Zhu and Jin, 2018;Mocanu et al., 2018), many times even higher than the fullyconnected counterparts (Mocanu et al., 2018), while being very versatile and suitable for many neural network models (e.g. restricted Boltzmann machines, multilayer perceptrons, and convolutional neural networks) and non-grid like data.
However, the authors of SET have used Keras with Theano back-end to implement their SET-MLP models. This implementation choice, while having the big advantage of offering a wide flexibility of architectural choices (e.g. various activation functions, optimizers, GPUs, and so on) which is very welcomed while conceiving new algorithms, does not offer proper support for sparse matrix operations. This limits considerably the practical aspects of SET-MLP with respect to its maximum possible number of neurons and implicitly to its representational power. Due to these reasons, the maximum size of the SET-MLPs from (Mocanu et al., 2018) is 12,082 neurons.

Proposed solution
In this paper, we address the above limitations of the SET original implementation and we show how vanilla SET-MLP can be implemented from scratch using just pure Python, SciPy, and Cython. Our approach enables the construction of SET-MLPs with at least two orders of magnitude larger, i.e. over 1,000,000 neurons. What is more, such SET-MLPs do not need GPUs and can run perfectly fine on a standard laptop.

SPARSE MATRICES OPERATIONS
The key element of our very efficient implementation is to use sparse data structures from SciPy. It is important to use the right representation of a sparse matrix for different operations, because different sparse matrix formats have different advantages and disadvantages (as briefly discussed in Appendix A) and one format cannot handle all operations necessary for sparse weights matrices to implement a SET-MLP. Still, the conversions from one format to another are very fast and efficient. Thus in our implementation which was done in pure Python 3, we have used for specific SET-MLP operations, specific sparse matrix formats and their fast conversion capabilities, as follows.
Initialize sparsely connected layers. The sparse matrices which store the sparsely connected layers are creating using the Linked List (LIL) format and then are transformed into Compressed Sparse Row (CSR) format.
Feed-forward phase. During it, the sparse weights matrices are stored and used in the CSR format.
Backpropagation phase -computing gradients. The only operations which can not be implemented with SciPy sparse matrix operations is computing the gradients for backpropagation (Rumelhart et al., 1986) due to the simple fact that by multiplying the vector of backpropagation errors from layer h k with the vector of activation neurons from layer h k−1 will perform a considerable amount of unnecessary multiplications (for nonexistent connections) and will create a dense matrix for updates. This dense matrix, besides being very slow to process, will have a quadratically number of parameters with respect to its number of rows and columns and will fill a 16GB RAM very fast (in practice, for less then 10000 neurons per layer given all the other necessary information which have to be stored in the computer memory).  (Haslinger et al., 2004) 111 11,340 3 5.9 MB SMK-CAN-187 (Zhao et al., 2010) 187 19,993 2 11.9 MB GLI-85 (Freije et al., 2004) 85 22,283 2 8.7 MB To avoid this situation, we have implemented in Cython the computations necessary for the batch weight updates. In this way, we compute in a much faster manner than in pure Python the gradient updates just for the existing connections. For this step, the sparse weight matrices are stored and used in the Coordinate list (COO) format.
Backpropagation phase -weights update. For this, the sparse weights matrices are used in the CSR format.

IMPLEMENTATION OF WEIGHTS EVOLUTION
The key aspect of the SET method that sets it apart from the conventional ANN training is the evolutionary scheme which modifies the connectivity of the layers at the end of every epoch. As the weight evolution routine is executed quite often, the routine needs to be implemented in an efficient manner to ensure that the SET-MLP training can be done as fast as possible. Furthermore, as the layer connections are extremely sparse in the SET scheme, the implementations should ensure that the sparsity level is maintained. Actually, it shall exploit the sparsity while removing and adding new weights. Two implementations of the weight evolution scheme were coded in native Python using Numpy sparse matrix routines. The first implementation is readable and intuitive, but does not exploit the full capabilities of the Numpy library in its various operations. The second implementation is not as readable, but vectorizes most of the operations using Numpy routines and performs the same operations in much lesser time. These implementations are both explained and compared in Appendix A.

Experiments and results
In this section, we evaluate and discuss the performance of our efficient SET-MLP implementation on four benchmark microarray datasets which are publicly available, as detailed in Table 1. For a good understanding of SET-MLP performance, we compare it against another sparse MLP model (implemented by us in the same manner) and in which the bipartite layers are initialized with an Erdős-Rényi topology, but which does not evolve over time and has a fixed sparsity pattern, dubbed MLP F ixP rob as in (Mocanu et al., 2018).

Evaluation metrics and experimental setup
To evaluate the performance of the proposed method we have used the accuracy metric and the confusion matrix to get detailed visual information. The rows of the confusion  matrix represent the predicted classes and the columns correspond to the true classes. The diagonal cells represent the numbers of samples that are correctly classified, whereas the off-diagonal cells are the number of incorrectly classified samples. The row at the bottom of the confusion matrix gives the proportion of all examples belonging to each class that are correctly (green) and incorrectly (red) classified. The column on the far right of the confusion matrix represents the proportion of all the samples predicted to belong to each class that are correctly (green) and incorrectly (red) classified.
All the experiments have been executed on a typical laptop using a single thread of the CPU. The laptop configuration is as follows: (1)  For both models, SET-MLP and MLP F ixP rob , we have used two hidden layers with ReLU activation functions, and backpropagation with stochastic gradient descent and momentum for the connection weights optimization. The numbers of neurons of each layer is given by Table 2.

Results on the Leukemia dataset
The Leukemia dataset is obtained from NCBI GEO repository with the accession number GSE13159. It contains 2096 samples with 54,675 features each. The samples are divided in 18 classes. Among these 2096 samples, 1397 samples are selected as training data and 699 as testing data. Table 3 shows the number of test samples in each class. It worths to be highlighted that both, the training and testing sets, are unbalanced. For this dataset, the number of hidden neurons in each layer was set to 27,500, a value which is way above the usual number of neurons in fully-connected MLP models. We performed a small random search, and we set the learning rate equaling 0.005, the batch size equaling 5, momentum to 0.9 and weight decay to 0.0002. As specific SET-MLP settings, we chose sparsity level = 10 and pruning rate ζ = 0.3.
The accuracy of SET-MLP and a MLP F ixP rob with two hidden layers on the Leukemia test set is shown in Fig. 1 over 500 training epochs. We mention that we chose to discuss in details the models with two hidden layer as they offered the best performance. The x-axis shows the training epochs; the y-axis shows the test accuracy. The figure shows that the accuracy of SET-MLP tends to stabilize towards 90% as the training progresses, while MLP F ixP rob stabilizes around 80% accuracy. To understand better the performance of our approach, Fig. 2 shows the confusion matrix on Leukemia dataset for the peak accuracy of SET-MLP (88.10%). We highlight that, to the best of our knowledge, this accuracy is higher than best results (81.11%) reported in the literature (Kumar and Rath, 2015) for this dataset. There in, an ensemble classifier is proposed to deal with microarray data. This classifier connects several feature selection algorithms with MapReduce based proximal support vector machine (mrPSVM) to classify the microarray data. Their experimental results not only show that the ensemble of mrPSVM classifier and feature selection approaches is a state-of-art method to deal with microarray datasets, but provides also concrete information about training data and testing data. In our experiments, we employed the same training and testing data to guarantee the validity of the comparison.

Results on the CLL-SUB-111 dataset
The CLL-SUB-111 is an unbalanced dataset contains gene expressions from high density oligonucleotide arrays, where   The test accuracy for 500 epochs of SET-MLP and MLP F ixP rob on the CLL-SUB-111 dataset is depicted in Fig. 3a. The x-axis shows the training epochs; the y-axis shows the test accuracy. After approximately 200 epochs, SET-MLP reaches an accuracy of 81.08% and MLP F ixP rob an accuracy of about 65%. Among the feature selection based methods to CLL-SUB-111, an accuracy of 78.38% was obtained by using Incremental Wrapper-based Attribute Selection(IWSS) (Bermejo et al., 2012). Although CLL-SUB-111 seriously suffers from an extreme small number of samples, we are still able to obtain outstanding performance with SET-MLP without any overfitting. Fig. 4a depicts the confusion matrix of SET-MLP on the CLL-SUB-111 dataset. We can observe that SET-MLP has excellent recall for class 1 (100.0%), even though there are extremely unfavorable conditions, i.e. very few training samples.

Results on the SMK-CAN-187 dataset
The SMK-CAN-187 is a RNA dataset obtained from normal bronchial epithelium of smokers with and without lung cancer, which is publicly available from (Zhao et al., 2010). It has 19,993 features and 187 samples. Out of these 187 samples, 124 samples are chose to be training data and 63 are testing data. There are 32 samples labeled as '1' and 31 are labeled as '2'. Same as Leukemia and CLL-SUB-111 datasets, SMK-CAN-187 is also unbalanced. The number of neurons in the hidden layers was set to 16,000, learning rate to 0.005, and the batch size to 5. Moreover, momentum was 0.9 and weight decay was set to 0.0002. Additionally, the SET-MLP specific parameters, and ζ were 10 and 0.3, respectively. Fig. 3b shows the accuracy performance of SET-MLP and MLP F ixP rob on the SMK-CAN-187 dataset. The accuracy of SET-MLP is 79.4% which is better than the accuracy reported in (Wang et al., 2016)(74.87±2.32%) in which feature selection was performed by preserving class correlation. It is noteworthy that on this dataset the accuracy of SET-MLP is stable to around 75% after 350 epochs, while MLP F ixP rob seems to suffer from overfitting. The confusion matrix of SET-MLP on SMK-CAN-187 is given in Fig.  4b. On class 2, SET-MLP achieves a remarkable performance of 90.3% recall.

Results on the GLI-85 dataset
The GLI-85 dataset has 22,283 features and 85 samples. We chose this dataset to analyze as it reflects an extreme case for the situation when very little labeled data is available. Out of these 85 samples, 56 samples are training data and 29 are testing data. There are 21 samples labeled as '2' and 8 are labeled as '1'. The batch size was set to 1, while the remaining hyperparameters were the same as for Leukemia. With only 85 samples, any model trained on this dataset is clearly prone to overfitting. Due to the very small number of samples, the best results on this dataset were obtained using just one hidden layer, as discussed next. Fig. 3c shows the accuracy performance of SET-MLP and MLP F ixP rob on the GLI-85 dataset. We can observe that peak accuracy of SET-MLP is 100.0% which is much better than the accuracy reported in (Taheri and Nezamabadi-pour, 2014)(94%) in which an ensemble including three filter methods with a meta-heuristic algorithm is used. It is noteworthy that on this dataset the accuracy of MLP F ixP rob is stable to 100%, while SET-MLP reaches 100% just in few cases, as it suffers some fluctuations. These suggest that MLP F ixP rob is even more stable than SET-MLP when a very small number of samples is available. The confusion matrix of SET-MLP on GLI-85 is given in Fig. 4c for the peak accuracy of 100%.

Results analysis
To understand better the connections reduction made by the SET procedure in a SET-MLP model in comparison with a fully-connected MLP (FC-MLP) which has the same amount of neurons, Fig. 5 and Table 4 provide the number of connections for the SET-MLP models discussed above       For a better understanding of SET computational requirements, Table 5 shows the average training and testing time per epoch of the SET-MLPs used on the datasets. We can observe, as expected, that as the number of features and samples increases also the training time is increasing. Still, it worths to be highlighted that, although, the average training time of Leukemia is relatively long (66.12s), it fulfills an almost impossible mission, that is, running such a large model on a commodity laptop.
In the paper, we have discussed the performance of the SET-MLP models with two hidden layers on the Leukemia, CLL-SUB-111, SMK-CAN-187 datasets and with one hidden layer on the GLI-85 datasets. In Appendix B, we explain our choices on the number of hidden layers by presenting comparatively the performance of SET-MLP models with one, two, and three hidden layers on all datasets and by discussing the beneficial effect of dropout (Hinton et al., 2012) on SET-MLP.

Discussion: extreme SET-MLP models
While in the previous section we have analyzed the qualitative performance of our proposed approach, in this section we briefly discuss two extreme SET-MLP models on the largest dataset used in this paper, i.e. Leukemia. The goal is to assess how fast SET-MLP can achieve a good accuracy and to see how large can be a trainable SET-MLP model on a typical laptop. For each model, we used a SET-MLP with two hidden layers, and a Softmax layer as output. For the first model, the number of hidden neurons per layer was set to 1,000, while for the second model the number of hidden neurons per layer was set to 500,000. In both cases, we have used a very eager learning rate (0.05) and we trained the models for 5 epochs. On each hidden layer we applied a dropout rate of 0.4. The other hyperparameters were set as in the previous section for Leukemia and we have used the same training/testing data splitting. Table 6 presents SET-MLP performance in comparison with the best state-of-the-art results of mrPSVM from (Kumar and Rath, 2015). We clarify that the goal of this experiment is not to obtain the best accuracy possible with SET-MLP.  (Goodfellow et al., 2016) the size of the largest neural networks which run currently in the cloud is about 10 to 20 million neurons. Therefore, our results emphasize even more the capabilities of SET-MLPs and open the path for new research directions.

Related work
Traditional feature selection methods can be roughly divided into three categories: filter methods, wrapper methods and embedded methods. Independent of classifiers, filter methods employ a certain criterion to rank different features and use the features with the highest scores to fulfill classification tasks. According to whether features are evaluated in an individual way or in a batch way, filter methods can be divided into univariate filters (Liu and Setiono, 1995) and multivariate filters (Hall, 1999). Considering the evaluation role of classifiers, wrapper methods employ classification performance as feedback to assess the selected features repeatedly, which, in turn, helps to find the best feature subset. The WrapperSubsetEval (Hall et al., 2009) is a general wrapper method which can be connected with various learning algorithms. Considering the advantages of filter methods and wrapper methods, embedded methods are capable of utilizing the biases of classifiers, while reducing the computation cost. Lasso regularization (Tibshirani, 1996) was employed to objective functions to eliminate the nonimportant features, by reducing the non-important features whose coefficients are close to zero. However, as the explosive growth of the internet, more and more datasets with ultrahigh dimensions are available, which are much out of the capacity of commodity computers. To address this problem, distributed computing has been proposed. Ibrahim et al. (Ibrahim et al., 2014) generalized Principal Components Analysis (ensemblePCA) to Deep Belief Network (DBN) and to transfer high-dimensional data to low-dimensional nodes. A distributed decentralized algorithm for k-Nearest Neighbor (kNN) graph (Plaku and Kavraki, 2007) was proposed to distribute the computation of kNN graph with ultra big datasets by utilizing the sequential structure of kNN data. MapReduce (Dean and Ghemawat, 2008) is an efficient programming model used by Google to compute different types of data and to process large raw data. Chu et al. (Chu et al., 2007) applied MapReduce to ten learning algorithms and obtained a progressive speed-up. Developed by Doug Cut- Table 6. Two extreme SET-MLP models on Leukemia against state-of-the-art (mrPSVM with ANOVA (Kumar and Rath, 2015) for feature selection). The numbers in brackets for SET-MLP reflect the number of neurons per layer from input to output. The accuracy of SET-MLP is reported as the mean and standard deviation of 5 runs. The density level represents the percentage of the number of existing connections in the SET-MLP model from the total number of connections in its corresponding FC-MLP.  (Kumar and Rath, 2015) conventional n/a 1265 81.1 mrPSVM with ANOVA (Kumar and Rath, 2015) Hadoop cluster n/a 291 81.1 ting, Apache Hadoop is an open-source software framework based on MapReduce that can accomplish data distributed storage and processing (White, 2012). It partitions big data into different blocks, in a distributed fashion, and then employs parallel tools to process data in different blocks on one big cluster. Kumar et al. (Kumar and Rath, 2015) proposed a classifier framework combining MapReduce and proximal Support Vector Machine(mrPSVM). The experimental results on several high-dimensional, low-sample benchmark datasets demonstrated that the ensemble of mrPSVM classifier and other feature selection methods outperforms the classical approaches.
As a promising method that has been widely used in image recognition, speech recognition and language translation, deep learning has also been employed in classification tasks with high-dimensional data. In (Fakoor et al., 2013), an autoencoder approach was combined with a softmax regression classifier to address cancer data classification problems. Danaee et la. (Ibrahim et al., 2014) introduced a multi-level feature selection method using a deep learning and active learning approach is proposed. The experiments verify the superiority of this approach in terms of classification accuracy. Additionally, a Stacked Denoising Autoencoder (SDAE) was employed to represent high-dimensional features with lower dimensional and more meaningful features (Danaee et al., 2017).
Although the above-mentioned algorithms can have a good performance on some datasets, one common disadvantage of all these algorithms is that the systems are hierarchical, meaning that at least a dimensionality reduction technique and an efficient classifier are needed. There is a serious dependency between the two subsystems. When one subsystem is not good enough, the overall performance can not be guaranteed due to information loss caused by the feature selection phase. That's why the aim of this paper was to introduce a method capable to process and classify highdimensional data with one unitary model, i.e. SET-MLP.

Conclusion
Microarray data have been treated in the literature as a difficult task due to their very high number of features but little number of examples. Besides that, this type of data suffer from imbalance and data shift problems.
In this paper, an efficient implementation of SET-MLP, a sparse multilayer perceptron trained with the sparse evolutionary training procedure, is proposed to deal with high dimensional microarray datasets. This implementation makes use just of Python 3, sparse data structures from SciPy, and Cython. With this implementation, we have created for the first time in literature sparse MLP models with over one million neurons which can be trained on a standard laptop using a single CPU thread and without GPU. This is with two orders of magnitude more than state-of-the-art fully connected MLPs and SET-MLPs trained on commodity hardware.
Besides that, we demonstrated on four microarray datasets with tens of thousands of input features and with up to just two thousands of samples that our approach reduces the number of connections quadratically in large MLPs (about 99.9 % connections reduction), while outperforming the state-of-the-art methods on these datasets for the classification task. Moreover, our proposed SET-MLP models showed to be robust to overfitting, imbalanced and data shift problems, which is not so usual for fully connected MLPs. The last but not the least, the results suggest that our proposed approach can cope efficiently with the "curse of dimensionality", being capable of learning from small amounts of labeled data, and outperforming the state-of-theart methods (ensembles of classifiers and feature selection methods) which are currently employed on high dimensional non-grid like data (or tabular data).
In the future, we intend to focus also on other types of neural layers, such as convolutional layers in CNN which have been widely-used to deal with graphic data with gridlike topology. Furthermore, we intend to extend this work to address problems from other fields which suffer from the "curse of dimensionality" and which have ultra high dimensional data (e.g. social networks, financial networks, semantic networks). The last but not the least future research direction, would be to parallelize our implementation to use efficiently all CPU threads of a typical workstation and to incorporate it into usual Deep Learning frameworks, such as TensorFlow or PyTorch. This probably would allow us to scale with one order of magnitude more the SET-MLP models (up to the level of few tens of millions of neurons), while still using commodity hardware.

A.1. Sparse data structures
Below the SciPy sparse data structures used to implement SET-MLPs are briefly discussed, while the interested reader is referred to 3 for detailed information.
• Compressed Sparse Row (CSR) sparse matrix: The data is stored in three vectors. The first vector contains nonzero values, the second one stores the extents of rows, and the third one contains the column indices of the nonzero values. This format is very fast for many arithmetic operations, but slow for changes to the sparsity pattern.
• Linked List (LIL) sparse matrix: This format saves nonzero values in row-based linked lists. Items in the rows are also sorted. The format is fast and flexible in changing the sparsity patterns, but inefficient for arithmetic matrix operations.
• Coordinate list (COO) sparse matrix: This format saves the nonzero elements and their coordinates (i.e. row and column). It is very fast in constructing new sparse matrices, but it does not support arithmetic matrix operations and slicing.
• Dictionary Of Keys (DOK) sparse matrix: This format has a dictionary that maps row and column pairs to the value of nonzero elements. It is very fast in incrementally constructing new sparse matrices, but can not handle arithmetic matrix operations.

A.2. Weights evolution -Implementation I
In this implementation, the sparse weight matrices in the CSR format are converted to three vectors representing the indices of the rows, columns of the non-zero elements along with the element values (either using the COO or LIL format). The values are then compared in a for-loop to the threshold to keep the weights or discard them, as per the user specified ζ values. To ensure that the total number of non-zeros in the weight matrix remains the same, random connections between neurons need to be created. Again a for-loop is used to create new random connections in an incremental manner and ensure that the total number of non-zeros are equal to the original number of non-zeros.
Most of the processing time in the code occurs in the forloops and the while loops and this is confirmed by a code profiling tool in python 4 . Furthermore, as we are constantly accessing the weights by the row and column index, this method does not exploit the sparsity of the weight matrix. The code profile of the processing time demonstrated that the removal of weights of the weight matrix takes about 15% of total time in an epoch and adding new random connections takes about 50% of the total time during an epoch (of course, these percentages depend also on the size of the datasets and they become smaller when the dataset gets larger). The detailed algorithm is given in Algorithm 2.

A.3. Weights evolution -Implementation II
In order to make full use of advantages of different sparse matrix formats, we also propose Fast Weights Evolution (FWE). In FWE, the sparse weight matrices in the CSR format are also converted to three vectors representing the indices of the rows, columns of the non-zero elements along with the element values using the COO format. The value vector is compared a single time with the minimum and maximum threshold values using the vectorized operations in numpy. This enables the identification of the indices of small weights for fast deletion of the weights. Next, the remaining row and column indices are stored together into an array and a list of all the arrays of the non-zero elements is created. This is used directly to determine the random row and column indices of the additional weights to ensure that the number of connections between the neurons are constant. As the weights are sparse, the size of the list is much smaller than the full size of the weight matrix and performing all the computations with the list will be faster. The detailed algorithm is given in Algorithm 3. The comparison of running time of these two implementations is given in Table 7, which shows that Implementation II is more efficient than Implementation I.
We know that it is hard for one to reproduce an efficient implementation of an algorithm given just the above details, and we mention that our SET-MLP proof-of-concept implementation is available online 5 B. Comparative study of SET-MLPs with one, two, and three hidden layers on the four benchmark datasets.
The amount of neurons per hidden layer and the other hyperparameters are set to be the same with the models from the  Figure 6. Experiments with SET-MLPs on all four datasets to understand the effect of the number of hidden layers(n h ). For each dataset, three cases for the number of hidden layers are considered, i.e. n h = {1, 2, 3}. Each row represents the test classification accuracy of SET-MLPs with one, two, or three hidden layers on the same dataset. Every model from each row has been trained with same hyperparameters as in the paper, except for number of hidden layers. paper. Fig. 6 summarizes these experiments. From the first row, it can be inferred that SET-MLP with two hidden layers reaches the highest peak accuracy (88.12%) and has relatively the most robust performance on the Leukemia dataset. Similarly, SET-MLP with two hidden layer reaches outstanding accuracy (81.11%) on the CLL-SUB-111 dataset, while the accuracy can not reach 80% with one or three hidden layers.
As expected, but at the same time having the most interesting results, due to the very small number of samples of GLI-85 (Fig. 6, third row), SET-MLP with one hidden layer avoids overfitting in exchange to quite an oscillating behavior. At the same time, SET-MLP with two or three hidden layers even if they are capable to reach also perfect accuracy of 100%, after about 200 epochs they have a dramatical drop in accuracy to about 80%. We hypothesis that this situation happens due to overfitting as the number of training samples is extremely insufficient. If this is the case, adding dropout regularization to SET-MLP is able to figure out this problem. We applied dropout with 0.5 dropout rate to both hidden layers. The performance is shown in Fig. 7. It is clear that the accuracy of SET-MLP with dropout keeps the same trend as before, without any drop in accuracy after 200 epochs. replace FC layer with Sparse Connected(SC) layer with a Erdős-Rényi topology given by and Eq.1; 6: end for 7: initialize training algorithm parameters; 8: %Training; 9: for each training epoch i do