BayesNAS: A Bayesian Approach for Neural Architecture Search

Introduction

Neural Architecture Search (NAS), the process of automating architecture engineering, is thus a logical next step in automating machine learning since (Zoph & Le, 2017). There are basically three existing frameworks for neural architecture search. Reinforcement learning based NAS (Baker et al., 2017; Zoph & Le, 2017; Zhong et al., 2018; Zoph et al., 2018; Cai et al., 2018) methods take the generation of a neural architecture as an agent’s action with the action space identical to the search space. More recent neuro-evolutionary approaches (Real et al., 2017; Liu et al., 2018b; Real et al., 2019; Miikkulainen et al., 2019; Xie & Yuille, 2017; Elsken et al., 2019a) use gradient-based methods for optimizing weights and solely use evolutionary algorithms for optimizing the neural architecture itself. However, these two frameworks take enormous computational power when compared to a search using a single GPU. One-Shot based NAS is a promising approach to significantly reduce search time without any separate training, which treats all architectures as different subgraphs of a supergraph (the one-shot model) and shares weights between architectures that have edges of this super-graph in common (Saxena & Verbeek, 2016; Brock et al., 2018; Pham et al., 2018; Bender et al., 2018; Liu et al., 2019b; Cai et al., 2019; Xie et al., 2019; Zhang et al., 2019a, b). A comprehensive survey on Neural Architecture Search can be found in (Elsken et al., 2019b).

Our approach is a one-shot based NAS solution which treats NAS as a Network Compression/pruning problem on the architecture parameters from an over-parameterized network. However, despite it’s remarkable less searching time compared to reinforcement learning and neuro-evolutionary approaches, we can identify a number of significant and practical disadvantages of the current one-shot based NAS. First, dependencies between a node and its predecessors and successors are disregarded in the process of identifying the redundant connections. This is mainly motivated by the improper treatment of zero operations. On one hand, the logit of zero may dominate some of the edges while the child network still has other non-zero edges to keep it connected (Liu et al., 2019b; Xie et al., 2019; Cai et al., 2019; Zhang et al., 2019b), for example, node 2 in Figure1a. Similarly, as shown in Figure 1 of (Xie et al., 2019), the probability of invalid/disconnected graph sampled will be $\frac{511}{1024}$ when there are three non-zero plus one zero operation. Though post-processing to safely remove isolated nodes is possible, e.g., for chain-like structure, it demands extensive extra computations to reconstruct the graph for complex search space with additional layer types and multiple branches and skip connections. This may prevent the use of modern network structure as the backbone such as DenseNet (Huang et al., 2017), newly designed motifs (Liu et al., 2018b) and complex computer vision tasks such as semantic segmentation (Liu et al., 2019a). On the other hand, zero operations should have higher priority to rule out other possible operations, since zero operations equal to all non-zero operations not being selected. Second, most one-shot NAS methods (Liu et al., 2019b; Cai et al., 2019; Xie et al., 2019; Zhang et al., 2019b; Gordon et al., 2018) rely on the magnitude of architecture parameters to prune redundant parts and this is not necessarily true. From the perspective of Network Compression (Lee et al., 2019), magnitude-based metric depends on the scale of weights thus requiring pre-training and is very sensitive to the architectural choices. Also the magnitude does not necessarily imply the optimal edge. Unfortunately, these drawbacks exist not only in Network Compression but also in one-shot NAS.

In this work, we propose a novel, efficient and highly automated framework based on the classic Bayesian learning approach to alleviate these two issues simultaneously. We model architecture parameters by a hierarchical automatic relevance determination (HARD) prior. The dependency can be translated by multiplication and addition of some independent Gaussian distributions. The classic Bayesian learning framework (MacKay, 1992a; Neal, 1995; Tipping, 2001) prevents overfitting and promotes sparsity by specifying sparse priors. The uncertainty of the parameter distribution can be used as a new metric to prune the redundant parts if its associated entropy $\frac{1}{2}\ln(2\pi e{\gamma_{jk}^{o^{\prime}}})$ is nonpositive. The majority of parameters are automatically zeroed out during the learning process.

Bayesian approach: BayesNAS is the first Bayesian approach for one-shot NAS. Therefore, our approach shares the advantages of Bayesian learning, which prevents overfitting and does not require tuning a lot of hyperparameters. Hierarchical sparse priors are used to model the architecture parameters. Priors can not only promote sparsity, but model the dependency between a node and its predecessors and successors ensuring a connected derived graph after pruning. Furthermore, it provides a principled way to prioritize zero operations over other non-zero operations. In our experiment on CIFAR-10, we found that the variance of the prior, as well as that of posterior, is several magnitudes smaller than posterior mean which renders a good metric for architecture parameters pruning.

Network compression: As a byproduct, our approach can be extended directly to Network Compression by enforcing various structural sparsity over network parameters. Extremely sparse models can be obtained at the cost of minimal or no loss in accuracy across all tested architectures. This can be effortlessly integrated into BayesNAS to find sparse architecture along with sparse kernels for resource-limited hardware.

Related Work

Network Compression. The de facto standard criteria to prune redundant weights depends on their magnitude and is designed to be incorporated with the learning process. These methods are prohibitively slow as they require many iterations of pruning and learning steps. One category is based on the magnitude of weights. The conventional approach to achieve sparsity is by enforcing penalty terms (Chauvin, 1989; Weigend et al., 1991; Ishikawa, 1996). Weights below a certain threshold could be removed. In recent years, impressive results have been achieved using the magnitude of weight as the criterion (Han et al., 2016) as well as other variations (Guo et al., 2016). The other category is based on the magnitude of Hessian of loss with respect to weights, i.e., higher the value of Hessian, greater the importance of the parameters (LeCun et al., 1990; Hassibi et al., 1993). Despite being popular, both of these categories require pretraining and are very sensitive to architectural choices. For instance, different normalization layers affect the magnitude of weights in different ways. This issue has been elaborated in (Lee et al., 2019) where the gradient information at the beginning of training is utilized for ranking the relative importance of weights’ contribution to the training loss.

Bayesian Learning and Compression. Our approach is based on Bayesian learning. In principle, the Bayesian approach to learn neural networks does not have problems of tuning a large amount of hyperparameters or overfitting the training data (MacKay, 1992b, a; Neal, 1995; Hernández-Lobato & Adams, 2015). Recently, Bayesian learning is applied to estimate layer size and network depth in NAS problem (Dikov et al., 2019). By employing sparsity-inducing priors, the obtained model depends only on a subset of kernel functions for linear models (Tipping, 2001) and deep neural networks where the neurons can be pruned as well as all their ingoing and outgoing weights (Louizos et al., 2017). Other Bayesian methods have also been applied to network pruning (Ullrich et al., 2017; Molchanov et al., 2017a) where the former extends the soft weight-sharing to obtain a sparse and compressed network and the latter uses variational inference to learn the dropout rate that can then be used for network pruning.

Search Space Design

The search space defines which neural architectures a NAS approach might discover in principle. Designing a good search space is a challenging problem for NAS. Some works (Zoph & Le, 2017; Zoph et al., 2018; Pham et al., 2018; Cai et al., 2018; Zhang et al., 2019b; Liu et al., 2019b; Cai et al., 2019) have proposed that the search space could be represented by a Directed Acyclic Graph (DAG). We denote $e_{ij}$ as the edge from node $i$ to node $j$ and $o_{ij}$ stands for the operation that is associated with edge $e_{ij}$.

Similar to other one-shot based NAS approaches (Bender et al., 2018; Zhang et al., 2019b; Liu et al., 2019b; Cai et al., 2019; Gordon et al., 2018), we also include (different or same) scaling scalars over all operations of all edges to control the information flow, denoted as $w_{ij}^{o}$ which also represent architecture parameters. The output of a mixed operation $o_{ij},i<j$ is defined based on the outputs of its edge

Then $z_{j}$ can be obtained as $\sum_{i<j}o_{j}(z_{i})$.

To this end, the objective is to learn a simple/sparse subgraph while maintaining/improving the accuracy of the over-parameterized DAG (Bender et al., 2018). Let us formulate the search problem as an optimization problem. Given a dataset $\mathbf{D}=(\mathbf{X},\mathbf{Y})=\{(\mathbf{x}_{n},\mathbf{y}_{n})\}_{n=1}^{N}$ and the desired sparsity level $\kappa$ (i.e., the number of non-zero edges), one-shot NAS problem can be written as an optimization problem with the following constraints:

To alleviate the negative effect induced by the dependency and magnitude-based metric whose issues have been discussed in Introduction, for each $w_{ij}^{o}$, we introduce a switch $s_{ij}^{o}$ that is analogous to the one used in an electric circuit. There are four features associated with these switches. First, the “on-off” status is not solely determined by its magnitude. Second, dependency will be taken into account, i.e., the predecessor has superior control over its successors as illustrated in Figure 1c. Third, $s_{ij}^{o}$ is an auxiliary variable that will not be updated by gradient descent but computed directly to switch on or off the edge. Lastly, $s_{ij}^{o}$ should work for both proxy and proxyless scenarios and can be better embedded into existing algorithmic frameworks (Liu et al., 2019b; Cai et al., 2019; Gordon et al., 2018). The calculation method will be introduced later in Section 4.

Inspired by the hierarchical representation in a DAG (Liu et al., 2019b, 2018b), we abstract a single motif as the building block of DAG, as shown in Figure 1e. Apparently, any derived motif, path, or network can be constructed by such a multi-input-multi-output motif. It shows that a successor can have multiple predecessors and each predecessor can have multiple operations over each of its successors. Since the representation is general, each directed edge can be associated with some primitive operations (e.g., convolution, pooling, etc.) and a node can represent output of motifs, cells, or a network.

Dependency Based One-Shot Performance Estimation Strategy

In the following, we will formally state the criterion to identify the redundant connections in Proposition 1. The idea can be illustrated by Figure 1b in which both the blue and red edges from node 2 to 3 and from node 2 to 4 might be non-zeros but should be removed as a consequence. To enable this, we have the following proposition.

There is information flow from node $j$ to $k$ under operation $o^{\prime}$ as shown in Figure 1e if and only if at least one operation of at least one predecessor of node $j$ is non-zero and $w_{jk}^{o^{\prime}}$ is also non-zero.

The same expression for Proposition 1 is: there is no information flow from node $j$ to $k$ under operation $o^{\prime}$ if and only if all the operation of all the predecessors of node $j$ are zeros or ${w_{jk}^{o^{\prime}}}$ is zero. This explains the incompleteness of the problem 2 as well as the possible phenomenon that non-zero edges become dysfunctional in Figure 1b.

As can be seen in Remark 2, we will construct a probability distribution jointly over $w_{jk}^{o^{\prime}}$, $w_{ij}^{o}$, $\forall i<j$ in the sequel, denoted as

where $c$ is a possible expression like in Remark 2 to encode Proposition 1.

In the following, we will show how the “switches” $s$ can be used to implement Proposition 1. If we assume $s$ has two states $\{\text{ON},\text{OFF}\}$, $w_{jk}^{o^{\prime}}$ is redundant when $s_{jk}^{o^{\prime}}$ is OFF or all $s_{ij}^{o}$ are OFF, $\forall i<j,o\in\mathcal{O}$. How to use $s$ to encode the redundancy of $w_{jk}^{o^{\prime}}$, i.e., ${w_{jk}^{o^{\prime}}\sum_{i<j}|w_{ij}^{o}|=0}$? One possible solution is

If $s$ is a continuous variable with $s=\infty$ for ON and for OFF, set union and intersection can be arithmetically represented by addition and multiplication respectively. $s$ does not directly determine the magnitude of $w$ but plays the role as uncertainty or confidence for zero magnitude.

A straightforward way to encode this logic is to assign a probability distribution, for example Gaussian distribution, over $w_{jk}^{o^{\prime}}$

Since $w_{ij}^{o},\forall i,j,o$ are independent with each other, we construct the following distribution to express equation 3:

Since $s_{ij}^{o}>0$ in equation 5 always holds, regardless of what $s_{ij}^{o}$ is, we can use the following simpler alternative to substitute equation 5 to encode Proposition 1:

Interestingly, equation 7 and 4 are equivalent. This means that we may find an algorithm that is able to find the sparse solution in a probabilistic manner. However, Gaussian distribution, in general, does not promote sparsity. Fortunately, some classic yet powerful techniques in Bayesian learning are applicable, i.e., sparse Bayesian learning (SBL) (Tipping, 2001; Pan, 2017) and automatic relevance determination (ARD) prior (MacKay, 1996; Neal, 1995) in Bayesian neural networks.

2 Zero Operation Ruling All

In our paper, we do not include zero operation as a primitive operation. Instead, between node $i$ and $j$ we compulsively add one more node $i^{\prime}$ and allow only a single identity operation (see Figure 1f). The associated weight $w_{ii^{\prime}}$ is trainable and initialized to $1$ as well as its switch $s_{ii^{\prime}}$. The idea is that if $s_{ii^{\prime}}$ is OFF, all the operations from $i^{\prime}$ to $j$ will be disabled as a consequence. Then $\gamma_{jk}^{o^{\prime}}$ in equation 6 can be substituted by

Bayesian Learning Search Strategy

The likelihood for the network weights $\bm{\mathcal{W}}$ and the noise precision $\sigma^{-2}$ with data ${\mathcal{D}=(\mathbf{X},\mathbf{Y})}$ is

To complete our probabilistic model, we specify a Gaussian prior distribution for each entry in each of the weight matrices in $\bm{\mathcal{W}}$. In particular,

where $\gamma_{jk}^{o^{\prime}}$ is defined in equation 8. $\sigma^{-2}$, $\lambda$ and $\mathbf{s}$ are hyperparameters. Importantly, there is an individual hyperparameter associated independently with every edge weight and a single one with all network weight. Follow Mackay’s evidence framework (MacKay, 1992a), ’hierarchical priors’ are employed on the latent variables using Gamma priors on the inverse variances. The hyper-priors for $\sigma^{-2}$, $\lambda$ and $\mathbf{s}$ are chosen to be a gamma distribution (Berger, 2013), i.e., $p(\lambda)=\text{Gam}(\lambda\,|\,a^{\lambda}$, $b^{\lambda}),p(\beta)=\text{Gam}(\beta\,|\,a^{\beta},b^{\beta})$ with $\beta=\sigma^{-2}$, and $p(s_{ij}^{o})=\text{Gam}(s_{ij}^{o}\,|\,a^{s_{ij}^{o}},b^{s_{ij}^{o}})$. Essentially, the choice of Gamma priors has the effect of making the marginal distribution of the latent variable prior the non-Gaussian Student’s t therefore promoting the sparsity (Tipping, 2001, Section 2 and 5.1). To make these priors non-informative (i.e., flat), we simply fix $a$ and $b$ to zero by assuming uniform scale priors for analysis and implementation. This formulation of prior distributions is a type of hierarchically constructed automatic relevance determination (HARD) prior which is built upon classic ARD prior (Neal, 1995; Tipping, 2001).

The posterior distribution for the parameters $\bm{\mathcal{W}}$, $\gamma$ and $\lambda$ can then be obtained by applying Bayes’ rule:

where $p(\mathbf{Y}\,|\,\mathbf{X})$ is a normalization constant. Given a new input vector $\mathbf{x}_{\star}$, we can make predictions for its output $\mathbf{y}_{\star}$ using the predictive distribution given by

where ${p(\mathbf{y}_{\star}|\mathbf{x}_{\star},\bm{\mathcal{W}},\mathbf{w},\lambda,\mathbf{s},\sigma^{2})=\mathcal{N}(\mathbf{y}_{\star}\,|\,\text{Net}(\mathbf{x}_{\star}),\sigma^{2})}$. However, the exact computation of $p(\bm{\mathcal{W}},\mathbf{w},\lambda,\mathbf{s},\sigma^{2}\,|\,\mathcal{D})$ and $p(\mathbf{y}_{\star}\,|\,\mathbf{x}_{\star})$ is not tractable in most cases. Therefore, in practice, we have to resort to approximate inference methods.

It should be noted that $\lambda$ is the same for all network parameters. However, it can be different for $\bm{\mathcal{W}}$ or constructed to represent the structural sparsity for Convolutional kernels in NN aiming for Network Compression, which is related to Bayesian compression (Louizos et al., 2017) and structural sparsity compression (Wen et al., 2016). We give some examples in Figure 2 and more can be found in the Appendix B.2 where extremely sparse networks on MNIST and CIFAR-10 can be obtained without accuracy deterioration. Since our main focus is on architecture parameters, without breaking the flow, we will fix $\lambda$ which is equivalent to the weight decay coefficient in SGD and $\sigma^{2}=0.01$ that is equivalent to the regularization coefficient for network parameters.

In case of uniform hyperpriors, we only need to maximize the term $p(\mathbf{Y}\,|\,\lambda,\mathbf{s},\sigma^{2})$ (MacKay, 1992a; Berger, 2013)

We assume that the distribution of data likelihood belongs to the exponential family

where $E_{D}(*)$ is the energy function over data.

2 Laplace Approximation and Efficient Hessian Computation

In related Bayesian models, the quantity in equation 14 is known as the marginal likelihood and its maximization is known as the type-II maximum likelihood method (Berger, 2013). And neural networks can also be treated in a Bayesian manner known as Bayesian learning for neural networks (MacKay, 1992b; Neal, 1995). Several approaches have been proposed based on, e.g., the Laplace approximation (MacKay, 1992b), Hamiltonian Monte Carlo (Neal, 1995), expectation propagation (Jylänki et al., 2014; Hernández-Lobato & Adams, 2015), and variational inference (Hinton & Van Camp, 1993; Graves, 2011). Among these methods, we adopt Laplace approximation. However, Laplace approximation requires computation of the inverse Hessian of log-likelihood, which can be infeasible to compute for large networks. Nevertheless, we are motivated by 1) its easy implementation, especially using recent popular deep learning open source software; 2) versatility for modern NN structures such as CNN and RNN as well as their modern variations; 3) close relationship between computation of Hessian and Network Compression using Hessian metric (LeCun et al., 1990; Hassibi et al., 1993); 4) acceleration effect to training convergence by second-order optimization algorithm (Botev et al., 2017) to which it is related. In this paper, we propose the efficient calculation/approximation of Hessian for convolutional layer and architecture parameter. The detailed calculation procedures are explained in Appendix C.2 and C.3 respectively.

3 Optimization Algorithm

As analyzed before, the optimization objective of searching architecture becomes removing redundant edges. The training algorithm is iteratively indexed by $t$. Each iteration may contain several epochs. The pseudo code is summarized in Algorithm 1. The cost function is simply maximum likelihood over the data $D$ with regularization whose intensity is controlled by the re-weighted coefficient $\omega$

The derivation can be found in Appendix A.1 and A.2. The algorithm mainly includes five parts. The first part is to jointly train $\mathcal{W}$ and $\mathbf{w}$. The second part is to freeze the architecture parameters and prepare to compute their Hessian. The third part is to update the variables associated with the architecture parameters. The fourth part is to prune the architecture parameters and the pruned net will be trained in a standard way in the fifth part. As discussed previously on the drawback of magnitude based pruning metric,

we propose a new metric based on maximum entropy of the distribution. Since $p({w_{jk}^{o^{\prime}}})$ in equation 5 is Gaussian with zero mean ${\gamma_{jk}^{o^{\prime}}}$ variance, the maximum entropy is $\frac{1}{2}\ln(2\pi e{\gamma_{jk}^{o^{\prime}}})$. We set the threshold for ${\gamma_{jk}^{o^{\prime}}}$ to prune related edges when $\frac{1}{2}\ln(2\pi e{\gamma_{jk}^{o^{\prime}}})\leq 0$, i.e., ${\gamma_{jk}^{o^{\prime}}}\leq 0.0585$.

The algorithm can be easily transferred to other scenarios. One scenario involves proxy tasks to find the cell. Similar to equation 16, we group same edge/operation in the repeated stacked cells where $g$ is the index. The cost function for proxy tasks is then given as follows in the form of re-weighted group Lasso:

The details are summarized in Algorithm 2 of Appendix A.3. Another scenario is on Network Compression with structural sparsity, which is summarized in Algorithm 3 of Appendix B.

Experiments

The experiments focus on two scenarios in NAS: proxy NAS and proxyless NAS. For proxy NAS, we follow the pipeline in DARTS (Liu et al., 2019b) and SNAS (Xie et al., 2019). First BayesNAS is applied to search for the best convolutional cells in a complete network on CIFAR-10. Then a network constructed by stacking learned cells is retrained for performance comparison. For proxyless NAS, we follow the pipeline in ProxylessNAS (Cai et al., 2019). First, the tree-like cell from (Cai et al., 2018) with multiple paths is integrated into the PyramidNet (Han et al., 2017). Then we search for the optimal path(s) within each cell by BayesNAS. Finally, the network is reconstructed by retaining only the optimal path(s) and retrained on CIFAR-10 for performance comparison. Detailed experiments setting is in Appendix D.1.

Unlike DARTS and SNAS that rely on validation accuracy during or after search, we use $\gamma$ in BayesNAS as performance evaluation criterion which enables us to achieve it in an one-shot manner.

Our setup follows DARTS and SNAS, where convolutional cells of 7 nodes are stacked for multiple times to form a network. The input nodes, i.e., the first and second nodes, of cell $k$ are set equal to the outputs of cell $k-1$ and cell $k-2$ respectively, with $1\times 1$ convolutions inserted as necessary, and the output node is the depthwise concatenation of all the intermediate nodes. Reduction cells are located at the 1/3 and 2/3 of the total depth of the network to reduce the spatial resolution of feature maps. Details about all operations included are shown in Appendix D.1. Unlike DARTS and SNAS, we exclude zero operations.

In the searching stage, we train a small network stacked by 8 cells using BayesNAS with different $\lambda_{w}$. This network size is determined to fit into a single GPU. Since we cache the feature maps in memory, we can only set batch size as 18. The optimizer we use is SGD optimizer with momentum 0.9 and fixed learning rate 0.1. Other training setups follow DARTS and SNAS (Appendix D.1). The search takes about $3$ hours on a single GPUAll the experiments were performed using NVIDIA TITAN V GPUs.

The normal and reduction cells learned on CIFAR-10 using BayesNAS are shown in Figure 3a and 3b. A large network of 20 cells where cells at 1/3 and 2/3 are reduction cells is trained from scratch with the batch size of 128. The validation accuracy is presented in Table 1. The test error rate of BayesNAS is competitive against state-of-the-art techniques and BayesNAS is able to find convolutional cells with fewer parameters when compared to DARTS and SNAS.

2 Proxyless Search

Using existing tree-like cell, we apply BayesNAS to search for the optimal path(s) within each cell. Varying from proxy search, cells do not share architecture in proxyless search.

The backbone used is PyramidNet with three layers each consisting of $18$ bottleneck blocks and $\alpha=84$. All $3\times 3$ convolution in bottleneck blocks are replaced by the tree-cell that has in total $9$ possible paths within. The groups for grouped convolution is set to $2$. For the detailed structure of the tree-cell, we refer to (Cai et al., 2018).

In the searching stage, we set batch size to 32 and learning rate to 0.1. We use the same optimizer as for proxy search. The $\lambda$ of BayesNAS for each possible path is set to $1\times 10^{-2}$.

Because each cell can have a different structure in proxyless setting, we demonstrate only two typical types of cell structure among all of them in Figure 4a and Figure 4b. The first type is a chain-like structure where only one path exists in the cell connecting the input of the cell to its output. The second type is an inception structure where divergence and convergence both exist in the cell. Our further observation reveals that some cells are dispensable with respect to the entire network. After the architecture is determined, the network is trained from scratch with the batch size of 64, learning rate as 0.1 and cosine annealing learning rate decay schedule (Loshchilov & Hutter, 2017). The validation accuracy is also presented in Table 1. Although test error increases slightly compared to (Cai et al., 2019), there is a significant drop in the number of model parameters to be learned which is beneficial for both training and inference.

Transferability to ImageNet

For ImageNet mobile setting, the input images are of size 224$\times$224. A network of 14 cells is trained for 250 epochs with batch size 128, weight decay $3\times 10^{-5}$ and initial SGD learning rate 0.1 (decayed by a factor of 0.97 after each epoch). Results in Table 2 show that the cell learned on CIFAR-10 can be transfered to ImageNet and is capable of achieving competitive performance.

Conclusion and Future Work

We introduce BayesNAS that can directly learn a sparse neural network architecture. We significantly reduce the search time by using only one epoch to get the candidate architecture. Our current implementation is inefficient by caching all the feature maps in memory to compute the Hessian. However, Hessian computation can be done along with backpropagation which will potentially further reduce the searching time and scale our approach to larger search spaces.

Acknowledgements

The work of Hongpeng Zhou is sponsored by the program of China Scholarships Council (No.201706120017).

References

Appendix A BayesNAS Algorithm Derivation

In this subsection, we explain the detailed algorithm of updating hyper-parameters for the abstracted single motif as shown in Figure 1e. The proposition about optimization objective will be illustrated firstly.

Suppose the likelihood of the architecture parameters of a neural network $\mathbf{w}$ could be formulated as one exponential family distribution $p(\mathbf{Y}\,|\,\mathbf{w},\mathbf{X},{\mathbf{s}})\sim\exp\left(-E_{D}(\mathbf{Y};\text{Net}(\mathbf{X};\mathbf{w});\mathbf{s})\right)$, where $\mathcal{D}=(\mathbf{X},\mathbf{Y})$ is the given dataset, $\mathbf{s}$ stands for the uncertainty and $E_{D}(*)$ represents the energy function over data. The sparse prior with super Gaussian distribution for each architecture parameter has been defined in equation 11. The unknown architecture parameter of the network $\mathbf{w}$ and hyperparameter ${\mathbf{s}}$ can be approximately obtained by solving the following optimization problem

specially, for the architecture parameter $\mathbf{w_{jk}^{o^{\prime}}}$ which is associated with one operation of the edge $e_{jk}$ ($j<k$), the optimization problem could be reformulated as:

where ${\mathbf{w}_{jk}^{o^{\prime}}}^{*}$ is arbitrary, and

It should also be noted that ${\mathbf{s}}_{jk}^{o^{\prime}}$ represents the uncertainty of ${\mathbf{w}_{jk}^{o^{\prime}}}$ without considering the dependency between edge $e_{jk}^{o^{\prime}}$ and $\sum_{i<j}e_{ij}^{o}$, where $o^{\prime}$ and $o$ stands for one possible operation in corresponding edges.

Given the likelihood with exponential family distribution

as explained in equation 5, we define the prior of $w_{jk}^{o^{\prime}}$ with Gaussian distribution

The marginal likelihood could be calculated as:

Typically, this integral is intractable or has no analytical solution.

The mean and covariance can be fixed if the family is Gaussian. Performing a Taylor series expansion around some point ${\mathbf{w}_{jk}^{o^{\prime}}}^{*}$, $E_{D}({\mathbf{w}_{jk}^{o^{\prime}}})$ can be approximated as

where $\mathbf{g}(\cdot)$ is the gradient and $\mathbf{H}(\cdot)$ is the Hessian of the energy function $E$

To derive the cost function in equation A.1.2, we introduce the posterior mean and covariance:

Now the approximated likelihood $p(\mathbf{Y}|{\mathbf{w}_{jk}^{o^{\prime}}})$ is a exponential of quadratic, then Gaussian,

We can write the approximate marginal likelihood as

where $\mathbf{m}_{jk}^{o^{\prime}}$ and $C_{jk}^{o^{\prime}}$ are given in equation A.1.6. From equation A.1.6a and equation A.1.6b, the data-dependent term can be re-expressed as

Using equation A.1.11, we can evaluate the integral in equation A.1.9 to obtain

Applying a $-2\log(\cdot)$ transformation to equation A.1.9, we have

Therefore we get the following cost function to be minimised in equation A.1.2 over ${\mathbf{w}_{jk}^{o^{\prime}}},{{\mathbf{s}}_{jk}^{o^{\prime}}},$

Once the estimation on ${\mathbf{w}_{jk}^{o^{\prime}}}$ and ${{\mathbf{s}}_{jk}^{o^{\prime}}}$ are obtained, the cost function is alternatively optimised. The new estimated ${\mathbf{w}_{jk}^{o^{\prime}}}$ can substitute ${\mathbf{w}_{jk}^{o^{\prime}}}^{*}$ and repeat the estimation iteratively.

We note that in equation A.1.5, ${\mathbf{w}_{jk}^{o^{\prime}}}^{*}$ may not be the mode (i.e., the lowest energy state), which means the gradient term $\mathbf{g}$ may not be zero. Therefore the selection of ${{\mathbf{w}_{jk}^{o^{\prime}}}(1)}^{*}$ remains to be problematic. We give the following Corollary to address this issue, which is more general.

Instead of minimising $\mathcal{L}({\mathbf{w}_{jk}^{o^{\prime}}},{{\mathbf{s}}_{jk}^{o^{\prime}}})$, we can solve the following optimization problem to get ${\mathbf{w}_{jk}^{o^{\prime}}},{{\mathbf{s}}_{jk}^{o^{\prime}}},$

We look at the first part of $\mathcal{L}({\mathbf{w}_{jk}^{o^{\prime}}},{{\mathbf{s}}_{jk}^{o^{\prime}}})$ in equation A.1.2, and define them as

Such quadratic approximation to $E_{D}({\mathbf{w}_{jk}^{o^{\prime}}})+{\mathbf{w}_{jk}^{o^{\prime}}}{{\mathbf{s}}_{jk}^{o^{\prime}}}^{-1}{\mathbf{w}_{jk}^{o^{\prime}}}$ is actually the same as the approximation procedure in Trust-Region Methods where a region is defined around the current iterate within which they trust the model to be an adequate representation of the objective function (Nocedal & Wright, 2006, pp.65).

To obtain each step, we seek a solution of the subproblem at iteration $t$

then inject ${\mathbf{w}_{jk}^{o^{\prime}}}^{*}$ into $\min_{{\mathbf{w}_{jk}^{o^{\prime}}},{{\mathbf{s}}_{jk}^{o^{\prime}}},}\mathcal{L}({\mathbf{w}_{jk}^{o^{\prime}}},{{\mathbf{s}}_{jk}^{o^{\prime}}})$, we can optimise equation A.1.16 instead of equation A.1.2, i.e., $\min_{{\mathbf{w}_{jk}^{o^{\prime}}},{{\mathbf{s}}_{jk}^{o^{\prime}}},}\hat{\mathcal{L}}({\mathbf{w}_{jk}^{o^{\prime}}},{{\mathbf{s}}_{jk}^{o^{\prime}}})$.

A.2 Algorithm for Proxyless Tasks

In this Section, we propose iterative optimization algorithms to estimate ${\mathbf{w}_{jk}^{o^{\prime}}}$ and ${{\mathbf{s}}_{jk}^{o^{\prime}}}$ alternatively.

We first target for the estimation of unknown parameter ${\mathbf{w}_{jk}^{o^{\prime}}}$ and hyperparameter ${\mathbf{s}}_{jk}^{o^{\prime}}$. In the sequel, we show that the stated program can be formulated as a convex-concave procedure (CCCP) for ${\mathbf{w}_{jk}^{o^{\prime}}}$ and ${\mathbf{s}}_{jk}^{o^{\prime}}$.

can be formulated as a convex-concave procedure (CCCP), where ${\mathbf{w}_{jk}^{o^{\prime}}}^{*}$ can be arbitrary real vector.

is convex jointly in ${\mathbf{w}_{jk}^{o^{\prime}}}$, ${{\mathbf{s}}_{jk}^{o^{\prime}}}$ due to the fact that $f(\mathbf{x},Y)=\mathbf{x}\mathbf{Y}^{-1}\mathbf{x}$ is jointly convex in $\mathbf{x}$, $\mathbf{Y}$ (see, (Boyd & Vandenberghe, 2004, p.76)). Hence $u$ as a sum of convex functions is convex.

is jointly concave in ${{\mathbf{s}}_{jk}^{o^{\prime}}}$, $\bm{\Pi}$. We exploit the properties of the determinant of a matrix

which is a $\log$-determinant of an affine function of semidefinite matrices $\bm{\Pi}$, ${{\mathbf{s}}_{jk}^{o^{\prime}}}$ and hence concave.

Therefore, we can derive the iterative algorithm solving the CCCP. We have the following iterative convex optimization program by calculating the gradient of concave part.

Using basic principles in convex analysis, we then obtain the following analytic form for the negative gradient of $v({{\mathbf{s}}_{jk}^{o^{\prime}}})$ at ${{\mathbf{s}}_{jk}^{o^{\prime}}}$ is (using chain rule):

Combined with equation A.1.6b, we denote a new hyper-parameter $w_{ij}^{o^{\prime}}$ as following:

Therefore, the iterative procedures equation A.2.5 and equation A.2.6 for ${\mathbf{w}_{jk}^{o^{\prime}}}(t)$ and ${{\mathbf{s}}_{jk}^{o^{\prime}}}(t)$ can be formulated as

the optimal ${{\mathbf{s}}_{jk}^{o^{\prime}}}(t)$ can be obtained as:

${\mathbf{w}_{jk}^{o^{\prime}}}(t)$ can be obtained as follows

We can then inject this into equation A.2.11, which yields

The update rules for ${\mathbf{s}}_{jk}^{o^{\prime}}$ without considering the dependency has been explained above. However, as illustrated in Sec .4, the dependency between a node and its predecessors should not be disregarded. It means the dependency between edge $e_{jk}^{o^{\prime}}$ and $\sum_{i<j}e_{ij}^{o}$ should be taken into consideration, then Gaussian prior could be defined as equation 5 and equation 11:

based on this prior, the uncertainty of ${\mathbf{w}_{jk}^{o^{\prime}}}(t)$ should be computed as:

A.3 Algorithm for Proxy Tasks

Our algorithm can be easily transferred to the scenario of proxy tasks to find the cell. Suppose a network is assembled by stacking $O$ different kinds of cells together, such as $\aleph_{1}$ normal cells and $\aleph_{O}$ reduction cells in (Liu et al., 2019b). Then optimal $O$ cells are required to be designed in a NAS task. As explained before, we design a switch ${\mathbf{s}}$ for each architecture parameter $w$ to determine the “on-off” of the corresponding edge in our method. In order to find such optimal cells, we propose that switches on the same position of the identical kind of cells should also be same. Based on this, the architecture parameters could be divided into different groups. The general grouped architecture parameters are given as follows:

Similar to equation A.2.11, if the group $o$ is consist of $\aleph_{o}$ elements, where $o=1,\ldots,O$, the optimal $s_{jk,o}^{o^{\prime}}$ can be obtained as:

The calculation of $\omega_{jk,g}^{o^{\prime}}$ for group $o$ is:

and both ${\mathbf{s}}$ and $\omega$ for the different elements in identity group should keep the same:

It should be noted that the detailed derivation procedures can be referred to A.1 and A.2. The pseudo code is summarised in Algorithm 2.

Appendix B Structural Bayesian Deep Compression

In addition to applying the proposed Bayesian approach to address NAS problem, we also explore the possibility of our method on network structural compression problem. In this section, we extend to compress deep neural networks by proposing a series of generic and easily implemented reweighted group Lasso algorithms to solve maximization of marginal likelihood $\int p(\mathbf{Y}\,|\,\bm{\mathcal{W}})p(\bm{\mathcal{W}})d\bm{\mathcal{W}}$ where $p(\bm{\mathcal{W}})$ can be specified as various sparse structured priors over network weights as shown in Table S3. The proposed Algorithm is generic for the weights in fully connected and convolutional neural networks. The training algorithm is iteratively indexed by $t$. Each iteration contains several epochs. Within each iteration $t$, there are three parts, The first part is simply a reweighted group Lasso type optimization. In the regularization terms $R(\omega^{l}\circ{\mathcal{W}^{l}})$, each weight of layer $l$ is scaled by a factor $(\omega^{l})$. The update of $(\omega^{l})$ needs the Hessian of each $\bm{\mathcal{W}}^{l}$. The Hessian of each layer can be computed recursively given the Hessian in the next layer through backward passing. The hyper-parameters $\gamma^{l}$, $C^{l}$ and $\omega^{l}$ will be updated every iteration $t$ with $T_{\max}$ being the maximal iterations. $\alpha^{l}$ is an introduced intermediate variable during the update process. The pseudo code is summarized in Algorithm 3.

B.2 Experiments

We first perform LeNet-300-100 and LeNet-5 on MNIST dataset (LeCun, 1998). For LeNet-300-100, we apply shape-wise, row-wise and column-wise regularization as shown in Fig. S5(a), S5(b) and S5(c), for the 2D weight matrices. The hyper-parameters $\gamma$, $\omega$ and $\alpha$ are updated every ten epochs for a total of $T_{\max}=10$ loops. The learned structure is $465-37-90$ with $1.54\%$ test error and $0.04$ FLOPS (Molchanov et al., 2017b). Comparison with other methods can be found in Table S4. For LeNet-5, we apply shape-wise and filter-wise regularization for the conv layer as shown in Fig. S5(a) and S5(j); row-wise and column-wise regularization for fc layer as shown in Fig. S5(b) and S5(c). The learned structure is $5-10-65-25$ with $1.00\%$ test error and $0.57$ FLOPS. Comparison with other methods can be found in Table S5.

B.2.2 ResNet-18 on CIFAR-10

We also evaluate our algorithm on Cifar10 dataset using ResNet-18 as initialized backbone (He et al., 2016). In addition to the input conv layer and output fc layer, the other 16 conv layers are separated into 8 blocks with 2 layers each. We apply shape-wise and filter-wise regularization to the conv layer as shown in Fig. S5(a) and S5(j); row-wise and column-wise regularization to the fc layer as shown in Fig. S5(b) and S5(c). The result is given in Table S6. It can be found that two Conv layers in block 4 are pruned away which shows the potential of our method to reduce the number of layers.

Appendix C Efficient Hessian Computation

The mathematical operation in a fully-connected layer could be formulated as:

where $h_{i}$ is the pre-activation value for node $i$ and $a_{i}$ is the activation value. $\sigma()$ is the element-wise activation function. ${\mathcal{W}_{ij}^{o}}$ stands for the weight matrix associated with operation $o$ in edge $e_{ij}^{o}$. In (Botev et al., 2017), a recursive method is proposed to compute the Hessian ${\mathbf{H}}$ for $\mathcal{W}_{ij}^{o}$:

where $\otimes$ stands for Kronecker product; The pre-activation Hessian $H_{j}^{o}$ is known and could be used to compute the pre-activation Hessian recursively for the previous layer:

In order to reduce computation complexity, the original pre-activation Hessian $H$ and Hessian ${\mathbf{H}}$ in Eq C.1.2-C.1.3 are replaced with their diagonal values for recursive computation. Thus the matrix multiplication could be reduced to vector multiplication. The hessian calculation process could be reformulated as:

C.2 Compute the Hessian of Conv Layer

Although the Hessian of weight matrix has been widely used in second-order optimization techniques to speed up the training process (LeCun et al., 1990; Amari, 1998), it still remains infeasible to calculate explicit Hessian directly due to the intensive computation burden (Martens & Grosse, 2015; Botev et al., 2017). Moreover, as most of current deep neural networks include plenty of Convolutional (Conv) layers, it further increases the difficulty of calculation due to the indirect convolution operation. Inspired by the Hessian calculation methods for Fully Connected (FC) layers as shown in (Botev et al., 2017), we propose a recursive and efficient method to compute the Hessian of Conv layers by converting Conv layers to FC layers (Ma & Lu, 2017). Therefore Hessian of the resulting equivalent FC layer is ready to be obtained. The detailed calculation procedures are explained in the following:

$({H^{oM}_{j}})^{n}$ is the pre-activation Hessian which could be computed recursively. With $({H^{oM}_{j}})^{n}$ known, the pre-activation Hessian for $(M_{i})^{n}$ could be calculated as:

where $(h_{i})^{n}$ is the pre-activation value for FC layer and $L$ means the loss function. The pre-activation Hessian $H^{M}_{i}$ could be obtained after concatenating all $(H^{M}_{i})^{n}$ as

the Hessian $\mathbf{H}^{oM}_{ij}$ for $\mathcal{W}^{oM}_{ij}$ can be obtained as:

As analyzed in Sec C.1, the Hessian calculation may cost a lot of time and resource. In order to address this problem, we propose the following approximate method:

C.3 Compute the Hessian of Architecture Parameter

After we have the computation method for the Hessian of a Convolutional layer, we need to consider the Hessian of an architecture parameter. Now the output from node $i$ to $j$ under operation $o$ becomes $w_{ij}^{o}B_{i}$, where $w_{ij}^{o}$ is the architecture parameter and $B_{i}$ stands for the input vector.

Since $B_{i}$ and $H_{j}$ are independent of each other, the Hessian $\mathbf{H}_{ij}^{o}$ could also be calculated more efficiently:

Appendix D Detailed Settings of Experiments

We employ the following techniques in our experiments: centrally padding the training images to $40\times 40$ and then randomly cropping them back to $32\times 32$; randomly flipping the training images horizontally; normalizing the training and validation images by subtracting the channel mean and dividing by the channel standard deviation.

The operations include: 3 $\times$ 3 and 5 $\times$ 5 separable convolutions, 3 $\times$ 3 and 5 $\times$ 5 dilated convolutions, 3 $\times$ 3 max pooling, 3 $\times$ 3 average pooling, and skip connection. All operations are of stride one (excluded the ones adjacent to the input nodes in the reduction cell, which are of stride two) and the convolved feature maps are padded to preserve their spatial resolution. Convolutions are applied in the order of BN-ReLU-Conv and the depthwise separable convolution is always applied twice (Zoph et al., 2018; Real et al., 2019; Liu et al., 2018a, 2019b).

The network parameters are optimized using momentum SGD, with initial learning rate $\eta_{\bm{\theta}}=0.1$, momentum 0.9, and weight decay $1\times 10^{-4}$. The batch size employed is 16 and the initial number of channels is 16.

D.2 Architecture evaluation on CIFAR-10

Following existing works (Zoph et al., 2018; Liu et al., 2018a; Pham et al., 2018; Real et al., 2019; Liu et al., 2019b), we employ the following additional enhancements: cutout (DeVries & Taylor, 2017).

D.3 Architecture transferability evaluation on CIFAR-10

The network is trained with batch size 128, SGD optimizer with weight decay $3\times 10^{-4}$, momentum 0.9 and initial learning rate 0.1, which is decayed using cosine annealing.