Data Augmentation for Bayesian Deep Learning

TL;DR

This paper introduces data augmentation techniques for Bayesian deep learning that improve uncertainty quantification and training efficiency by leveraging scale mixture models and advanced algorithms.

Contribution

It develops novel data augmentation strategies for deep learning using scale mixture theory, enabling scalable EM and MCMC methods for uncertainty quantification.

Findings

Enhanced training efficiency over traditional SGD

Effective uncertainty quantification in deep learning models

Scalable algorithms for high-dimensional models

Abstract

Deep Learning (DL) methods have emerged as one of the most powerful tools for functional approximation and prediction. While the representation properties of DL have been well studied, uncertainty quantification remains challenging and largely unexplored. Data augmentation techniques are a natural approach to provide uncertainty quantification and to incorporate stochastic Monte Carlo search into stochastic gradient descent (SGD) methods. The purpose of our paper is to show that training DL architectures with data augmentation leads to efficiency gains. We use the theory of scale mixtures of normals to derive data augmentation strategies for deep learning. This allows variants of the expectation-maximization and MCMC algorithms to be brought to bear on these high dimensional nonlinear deep learning models. To demonstrate our methodology, we develop data augmentation algorithms for a variety of commonly used activation functions: logit, ReLU, leaky ReLU and SVM. Our methodology is compared to traditional stochastic gradient descent with back-propagation. Our optimization procedure leads to a version of iteratively re-weighted least squares and can be implemented at scale with accelerated linear algebra methods providing substantial improvement in speed. We illustrate our methodology on a number of standard datasets. Finally, we conclude with directions for future research.

Tables

Table 1. Table 1: Data Augmentation Strategies. Here ρ τ ​ ( x ) = 1 2 ​ | x | + ( τ − 1 2 ) ​ x subscript 𝜌 𝜏 𝑥 1 2 𝑥 𝜏 1 2 𝑥 \rho_{\tau}(x)=\frac{1}{2}|x|+\left(\tau-\frac{1}{2}\right)x is the check function.

$l(W,b)$	$Q(W,b,\omega )$	$p\left(\omega \right)$
ReLU: $\max (1-z_i,0)$	${\int }_0^{\infty }\frac{1}{\sqrt{2\pi c\lambda }}\exp \left\{-\frac{{\left(x+a\lambda \right)}^2}{2c\lambda }\right\}d\lambda =\frac{1}{a}\exp \left(-\frac{2\max (ax,0)}{c}\right)$	$𝒢ℐ𝒢(1,0,0)$
Logit: $\log \left(1+e^{z_i}\right)$	$\frac{1}{2^b}{e}^{\left(a-b/2\right)\psi }{\int }_0^{\infty }{e}^{-\omega {\psi }^2/2}p\left(\omega \right)d\omega =\frac{{\left(e^{\psi }\right)}^a}{{\left(1+e^{\psi }\right)}^b}$	$𝒫𝒢(b,0)$
Lasso: $\left|\frac{z_i}{\sigma }\right|$	${\int }_0^{\infty }\frac{1}{\sqrt{2\pi c\lambda }}\exp \left\{-\frac{x^2}{2c\lambda }\right\}{e}^{-\frac{1}{2}\lambda }d\lambda =\frac{1}{c}\exp \left(-\frac{\left|x\right|}{c}\right)$	$ℰ\left(\frac{1}{2}\right)$
Check: $\left|z_i\right|+\left(2\tau -1\right)z_i$	${\int }_0^{\infty }\frac{1}{\sqrt{2\pi c\lambda }}\exp \left\{-\frac{{\left(x+\left(2\tau -1\right)\lambda \right)}^2}{2c^2\lambda }\right\}{e}^{-2\tau \left(1-\tau \right)\lambda }d\lambda =\frac{1}{c}\exp \left(-\frac{2}{c}{\rho }_{\tau }\left(x\right)\right)$	$𝒢ℐ𝒢(1,0,2\sqrt{\tau -{\tau }^2})$

Table 2. Table 2: Frequencies of Different Wine Ratings

rating	3	4	5	6	7	8	9
frequency	20	163	1457	2198	880	175	5

Table 3. Table 3: Percentage of Each Class (#obs = 21 3451)

	AU	CA	DE	ES	FR	UK	IT	NDF	NL	PT	US	other
% obs	0.25	0.67	0.50	1.05	2.35	1.09	1.33	58.35	0.36	0.10	29.22	4.73

Equations

$$\mathbf{y}=f_{\mathbf{\theta }}(\mathbf{x}),$$

$$y_{\star }=f_{\hat{\mathbf{\theta }}}({\mathbf{x}}_{\star }),$$

$$\mathbf{y}=f_{\mathbf{\theta }}(\mathbf{x})=(f_{W_0,b_0}\circ f_{W_1,b_1}\circ \cdots \circ f_{W_L,b_L})(\mathbf{x}),$$

$$\begin{array}{rl}\mathbf{y}& =f_0(Z_1{W}_0+b_0),\\ Z_l& =f_l(Z_{l+1}{W}_l+b_l),l=1,\dots ,L,\\ Z_{L+1}& =\mathbf{x},\end{array}$$

$$\hat{\mathbf{\theta }}:=\arg \underset{\mathbf{\theta }}{\min }\sum _{i=1}^n\ell (y_i,f_{\mathbf{\theta }}({\mathbf{x}}_i))+\lambda \sum _{j=1}^{\#\mathbf{\theta }}\varphi ({\theta }_j),$$

$$\hat{\mathbf{\theta }}=\arg \underset{\mathbf{\theta }}{\min }\frac{1}{n}\sum _{i=1}^n\ell (y_i,f_{\mathbf{\theta }}(x_i))+\sum _{l=0}^L{\lambda }_l{\varphi }_l(W_l,b_l).$$

$$\hat{\mathbf{\theta }}=\arg \underset{\mathbf{\theta }}{\min }\left\{\frac{1}{n}\sum _{i=1}^n\ell (y_i,f_{\mathbf{\theta }}(x_i))+\sum _{l=0}^L{\lambda }_l{\varphi }_l(W_l,b_l)\right\}$$

$$\arg \underset{\mathbf{\theta }}{\max }p(\mathbf{\theta }|\mathbf{y})=\arg \underset{\mathbf{\theta }}{\max }\exp \left\{-\frac{1}{n}\sum _{i=1}^n\ell (y_i,f_{\mathbf{\theta }}(x_i))-\sum _{l=0}^L{\lambda }_l{\varphi }_l(W_l,b_l)\right\},$$

$$p(\mathbf{\theta }\mid \mathbf{y})=p(\mathbf{y}\mid \mathbf{\theta })p(\mathbf{\theta })/p(\mathbf{y}),$$

$$\displaystyle p({\boldsymbol{y}}|\boldsymbol{\theta})\propto\exp\{-\sum_{i=1}^{n}\ell\big{(}y_{i},f_{\boldsymbol{\theta}}(\boldsymbol{x}_{i})\big{)}\},\quad p(\boldsymbol{\theta})\propto\exp\{-\sum_{l=0}^{L}\lambda_{l}\phi_{l}(W_{l},b_{l})\}.$$

$$\begin{array}{rl}\mathbf{y}\mid Z_1& \sim p(\mathbf{y}\mid Z_1),\\ Z_l& \sim N(f_l(W_l{Z}_{l+1}+b_l),{\tau }_l^2),l=1,2,\dots ,L,\\ Z_{L+1}& =\mathbf{x}.\end{array}$$

$$\mathbf{y}=f_0(Z_1{W}_0+b_0),Z_1=f_{\mathbf{B}}(\mathbf{x}),$$

$$p(W_0,b_0\mid Z_1^{(t)},\mathbf{y})$$

$$\propto \exp \left\{-\frac{1}{n}\sum _{i=1}^n\ell (y_i,f_{\mathbf{\theta }}(x_i)\mid {\mathbf{B}}^{(t)})+{\lambda }_0{\varphi }_0(W_0,b_0)\right\}$$

$${\mathbf{B}}^{(t+1)}$$

$$\displaystyle=\arg\min_{\boldsymbol{B}}\frac{1}{n}\sum_{i=1}^{n}\ell\big{(}Z_{1}^{(t)},f_{\boldsymbol{B}}(x_{i})\big{)}+\sum_{l=1}^{L}\lambda_{l}\phi_{l}(W_{l},b_{l}).$$

$$p(\boldsymbol{\theta}\mid{\boldsymbol{y}})=E_{\boldsymbol{\omega}}\Big{[}p(\boldsymbol{\theta},\boldsymbol{\omega}\mid{\boldsymbol{y}})\Big{]},$$

$$\displaystyle\exp\Big{\{}-\ell\big{(}\boldsymbol{y},f_{\boldsymbol{\theta}}(\boldsymbol{x})\big{)}\Big{\}}$$

$$\displaystyle=\int_{0}^{\infty}\exp\Big{(}-Q\big{(}\boldsymbol{y}\mid f_{\boldsymbol{\theta}}(\boldsymbol{x}),\boldsymbol{\omega}\big{)}\Big{)}p(\boldsymbol{\omega})\mathrm{d}\,\boldsymbol{\omega}$$

$$\exp (-\max (1-x,0))$$

$$\exp (-\log (1+e^x))$$

$$\exp (-|x|)$$

$$\hat{\boldsymbol{\theta}}\mathrel{\mathop{:}}=\arg\max_{\boldsymbol{\theta}}E_{\boldsymbol{\omega}}\Big{[}\exp\Big{(}-\frac{1}{n}\sum_{i=1}^{n}Q\big{(}y_{i}\mid f_{\boldsymbol{\theta}}(\boldsymbol{x}_{i}),\boldsymbol{\omega}\big{)}-\sum_{l=0}^{L}\lambda_{l}\phi_{l}(W_{l},b_{l})\Big{)}\Big{]},$$

$$H(\boldsymbol{\theta})=\mathop{E}\nolimits_{\boldsymbol{\omega}\mid\boldsymbol{\theta}}\Big{[}\exp\Big{(}-\frac{1}{n}\sum_{i=1}^{n}Q\big{(}y_{i}\mid f_{\boldsymbol{\theta}}(\boldsymbol{x}_{i}),\boldsymbol{\omega}\big{)}-\sum_{l=0}^{L}\lambda_{l}\phi_{l}(W_{l},b_{l})\Big{)}\Big{]},$$

$$G\big{(}\boldsymbol{\theta}\mid\boldsymbol{\theta}^{(t)}\big{)}=-\mathop{E}\nolimits_{\boldsymbol{\omega}\mid\boldsymbol{\theta}^{(t)}}\left[\frac{1}{n}\sum_{i=1}^{n}Q\big{(}y_{i}\mid f_{\boldsymbol{\theta}}(\boldsymbol{x}_{i}),\boldsymbol{\omega}\big{)}+\sum_{l=0}^{L}\lambda_{l}\phi_{l}(W_{l},b_{l})\right],$$

$$\log H(\boldsymbol{\theta})\geq G\big{(}\boldsymbol{\theta}\mid\boldsymbol{\theta}^{(t)}\big{)}.$$

$$p(\mathbf{\theta }\mid \mathbf{y})$$

$$\displaystyle=E_{\boldsymbol{\omega}}\Big{[}\exp\Big{(}-\frac{1}{n}\sum_{i=1}^{n}Q\big{(}y_{i}\mid f_{\boldsymbol{\theta}}(\boldsymbol{x}_{i}),\boldsymbol{\omega}\big{)}-\sum_{l=0}^{L}\lambda_{l}\phi_{l}(W_{l},b_{l})\Big{)}\Big{]}$$

$$\displaystyle=\int_{0}^{\infty}\exp\Big{(}-\frac{1}{n}\sum_{i=1}^{n}Q(y_{i}\mid f_{\boldsymbol{\theta}}(\boldsymbol{x}_{i}),\boldsymbol{\omega}\big{)}\Big{)}p(\boldsymbol{\omega})p(\boldsymbol{\theta})\mathrm{d}\,\boldsymbol{\omega}$$

$$p(\boldsymbol{\theta},\boldsymbol{\omega}\mid{\boldsymbol{y}})\propto\exp\Big{(}-Q(\boldsymbol{y}\mid f_{\boldsymbol{\theta}}(\boldsymbol{x}),\boldsymbol{\omega})\Big{)}p(\boldsymbol{\theta})p(\boldsymbol{\omega}).$$

$$\displaystyle p\big{(}\boldsymbol{\theta}^{(t)}\mid\boldsymbol{\omega}^{(t)},\boldsymbol{y}\big{)}$$

$$\displaystyle p\big{(}\boldsymbol{\omega}^{(t+1)}\mid\boldsymbol{\theta}^{(t)},\boldsymbol{y}\big{)}$$

$$p\big{(}y_{\star}\mid f(\boldsymbol{x}_{\star})\big{)}=\int p\big{(}y_{\star}\mid\boldsymbol{\theta},f_{\boldsymbol{\theta}}(\boldsymbol{x}_{\star})\big{)}p(\boldsymbol{\theta}\mid\boldsymbol{y})\mathrm{d}\,\boldsymbol{\theta}\approx\frac{1}{T}\sum_{t=1}^{T}p\big{(}y_{\star}\mid\boldsymbol{\theta}^{(t)},f_{\boldsymbol{\theta}^{(t)}}(\boldsymbol{x}_{\star})\big{)}.$$

$${\pi }_J(\mathbf{\theta })=\exp \left\{-Jf(\mathbf{\theta })\right\}/Z_J\mathrm{for}\mathbf{\theta }\in \mathbf{\Theta }$$

$$Z_1^j|\mathbf{\theta },\mathbf{x},\mathbf{y}\stackrel{iid}{\sim }\mathcal{N}({\mu }_z,{\sigma }_z^2),j=1,\dots ,J,$$

Full text

deep learning,

data augmentation,

MCMC,

back-propagation,

SGD,

\arxiv

1903.09668

\startlocaldefs

\endlocaldefs

,

and

1 Introduction

Deep neural networks (DNNs) have become a central tool for Artificial Intelligence (AI) applications such as, image processing (ImageNet, Krizhevsky et al. (2012)), object recognition (ResNet, He et al. (2016)) and game intelligence (AlphaGoZero, Silver et al. (2016)). The approximability (Poggio et al., 2017; Bauer and Kohler, 2019) and rate of convergence of deep learning, either in the frequentist fashion (Schmidt-Hieber, 2020) or from a Bayesian predictive point of view (Polson and Rockova, 2018; Wang and Rockova, 2020), have been well-explored and understood. Fan et al. (2021) provides a selective overview of deep learning. However, training deep learners is challenging due to the high dimensional search space and the non-convex objective function. Deep neural networks have also suffered from issues such as local traps, miscalibration and overfitting. Various efforts have been made to improve the generalization performance and many of their roots lie in Bayesian modeling. For example, Dropout (Wager et al., 2013) is commonly used and can be viewed as a deterministic ridge $\ell_{2}$ regularization. Sparsity structure via spike-and-slab priors (Polson and Rockova, 2018) on weights helps DNNs adapt to smoothness and avoid overfitting. Rezende et al. (2014) propose stochastic back-propagation through the use of latent Gaussian variables.

In this paper, following the spirit of hierarchical Bayesian modeling, we develop data augmentation strategies for deep learning with a complete data likelihood function equivalent to weighted least squares regression. By using the theory of mean-variance mixtures of Gaussians, our latent variable representation brings all of the conditionally linear model theory to deep learning. For example, it allows for the straightforward specification of uncertainty at each layer of deep learning and for a wide range of regularization penalties. Our method applies to commonly used activation functions such as ReLU, leaky ReLU, logit (see also Gan et al. (2015)), and provides a general framework for training and inference in DNNs. It inherits the advantages and disadvantages of data augmentation schemes. For approximation methods like Expectation-Maximization (EM) and Minorize-Maximization (MM), they are stable as they increase the objective but can be slow in the neighborhood of the maximum point even with acceleration methods such as Nesterov acceleration available and the performance is highly dependent on the properties of the objective function. Stochastic exploratory methods like MCMC have the main advantage of addressing uncertainty quantification (UQ) and are stable in the sense they require no tuning. Hyper-parameter estimation is immediately available using traditional Bayesian methods. DA augments the objective function with extra hidden units which allow for efficient step size selection for the gradient descent search. In some of the applications, data augmentation methods can be formulated in terms of complete data sufficient statistics, a considerable advantage when dealing with large datasets where most of the computational expense comes from repeatedly iterating over the data. By combing the MCMC methods with the J-copies trick (Jacquier et al., 2007), we can move faster towards posterior mode and avoid local maxima. Traditional methods for training deep learning models such as stochastic gradient descent (SGD) have none of the above advantages. We also note that we exploit the advantages of SGD and accelerated linear algebra methods when we implement our weighted least squares regression step.

Data augmentation strategies are commonplace in statistical algorithms and accelerated convergence (Nesterov, 1983; Green, 1984) is available. Our goal is to show similar efficiency improvements for deep learning. Our work builds on Deng et al. (2019) who provide adaptive empirical Bayes methods. In particular, we show how to implement standard activation functions, including ReLU (Polson and Rockova, 2018), logistic (Zhou et al., 2012; Hernández-Lobato and Adams, 2015) and SVM (Mallick et al., 2005) activation functions and provide specific data augmentation strategies and algorithms. The core subroutine of the resulting algorithms solves a least squares problem. Scalable linear algebra libraries such as Compute Unified Device Architecture (CUDA) and accelerated linear algebra (XLA) are available for implementation. To illustrate our approach, empirically we experiment with two benchmark datasets using Pólya-Gamma data augmentation for logit activation functions. For the deep architecture embedded in our approach, we adopt deep ReLU networks. Deep networks are able to achieve the same level of approximation accuracy with exponentially fewer parameters for compositional functions (Mhaskar et al., 2017). Poggio et al. (2017) further show how deep networks can avoid the curse of dimensionality. The ReLU function is favored due to its ability to avoid vanishing gradients and its expressibility and inherent sparsity. Approximation properties of deep ReLU networks have been developed in Montufar et al. (2014), Telgarsky (2017), and Liang and Srikant (2017). Yarotsky (2017) and Schmidt-Hieber (2020) show that deep ReLU networks can yield a rate-optimal approximation of smooth functions of an arbitrary order. Polson and Rockova (2018) provide posterior rates of convergence for sparse deep learning.

There is another active area of research that revives traditional statistical models with the computational power of DL (Bhadra et al., 2021). Examples include Gaussian Process models (Higdon et al., 2008; Gramacy and Lee, 2008), Generalized Linear Models (GLM) and Generalized Linear Mixed Models (GLMM) (Tran et al., 2020) and Partial Least Squares (PLS) (Polson et al., 2021). Our method benefits from the computation efficiency and flexibility of expression of the deep neural network. In addition, our work builds on the sampling optimization literature (Pincus, 1968, 1970) which now uses MCMC methods. Other examples include Ma et al. (2019) who study that sampling can be faster than optimization and Neelakantan et al. (2017) showing that gradient noise can improve learning for very deep networks. Gan et al. (2015) implements data augmentation inside learning deep sigmoid belief networks. Neal (2011) and Chen et al. (2014) provide Hamitonian Monte Carlo (HMC) algorithms for MCMC. Duan et al. (2018) proposes a family of calibrated data-augmentation algorithms to increase the effective sample size.

The rest of our paper is outlined as follows. Section 2 provides the general setting of deep neural networks and shows how DA can be integrated into deep learning using the duality between Bayesian simulation and optimization. Section 3 describes our data augmentation (DA) schemes and two approaches to implement them. Section 4 provides applications to Gaussian regression, support vector machines and logistic regression using Pólya-Gamma augmentation (Polson et al., 2013). Section 5 provides the experiments of DA on both regression and classification problems. Section 6 concludes with directions for future research.

2 Bayesian Deep Learning

In deep learning we wish to recover a multivariate predictive map $f_{\boldsymbol{\theta}}(\cdot)$ denoted by

[TABLE]

where ${\boldsymbol{y}}=(y_{1},\ldots,y_{n})^{\prime},y_{i}\in\mathbb{R}$ denotes a univariate output and ${\boldsymbol{x}}=(\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{n})^{\prime}$, $\boldsymbol{x}_{i}\in\mathbb{R}^{p}$ a high-dimensional set of inputs. Using training data of input-output pairs $\{y_{i},{\boldsymbol{x}}_{i}\}_{i=1}^{n}$ that generalizes well out-of-sample, the goal is to provide a predictive rule for a new input variable ${\boldsymbol{x}}_{\star}$

[TABLE]

where $\hat{\boldsymbol{\theta}}$ is estimated from training data typically using SGD. The interest in deep learners lies in their ability to perform better than the additive rule for those interpolation or prediction problems. Other statistical alternatives include Gaussian processes but they often have difficulty in handling higher dimensions.

Deep learners use compositions (Kolmogorov, 1957; Vitushkin, 1964) of ridge functions rather than additive functions that are commonplace in statistical applications. With $L\in\mathbb{N}$ we denote the number of hidden layers and with $p_{l}\in\mathbb{N}$ the number of neurons at the $l^{th}$ layer. Setting $p_{L+1}=p,p_{0}=p_{1}=1$, we denote with $\mathbf{p}=(p_{0},p_{1},\ldots,p_{L+1})\in\mathbb{N}^{L+2}$ the vector of neuron counts for the entire network. Imagine composing $L$ layers, a deep predictor then takes the form

[TABLE]

where $b_{l}\in\mathbb{R}^{p_{l}}$ is a shift vector, $W_{l}\in\mathbb{R}^{p_{l-1}\times p_{l}}$ is a weight matrix that links neurons between $(l-1)^{th}$ and $l^{th}$ layers and $f_{W_{l},b_{l}}(x)=f_{l}(W_{l}x+b_{l})$ is a semi-affine function. We denote with $\boldsymbol{\theta}=\{(W_{0},b_{0}),(W_{1},b_{1}),\ldots,(W_{L},b_{L})\}$ as the stacked parameters. We can rewrite the compositions in (2.1) with a set of latent variables $Z=(Z_{1},Z_{2},$ $\ldots,Z_{L})^{\prime}$ as

[TABLE]

where $Z_{l}\in\mathbb{R}^{n\times p_{l}}$ is the matrix of hidden nodes in $l$-th layer. We only consider the case $p=1$ and $Z_{1}\in\mathbb{R}^{n}$ in our work. We provide discussion on extensions to cases $p>1$ for some of our applications in Section 4.

2.1 Bayesian Simulation and Regularization Duality

The problem of deep learning regularization (Polson and Sokolov, 2017) is to find a set of parameters $\boldsymbol{\theta}$ which minimizes a combination of a negative log-likelihood $\ell(\boldsymbol{y},f_{\boldsymbol{\theta}}(\boldsymbol{x}))$ and a penalty function $\phi(\boldsymbol{\theta})$ defined by

[TABLE]

where $\lambda$ controls regularization and $\#{\boldsymbol{\theta}}$ denotes the number of parameters in $\boldsymbol{\theta}$.

When the function $f_{\boldsymbol{\theta}}(\boldsymbol{x})$ is a deep learner defined as (2.1), we can specify different amount of penalty $\lambda_{l}$ and form of regularization function $\phi_{l}(\cdot)$ for each layer. Then the objective function can be written as

[TABLE]

Commonly used regularization techniques for deep learners include $L^{2}$(weight decay), spike-and-slab regularization (Polson and Rockova, 2018) and dropout (Wager et al., 2013), which can also be viewed as a variant of $L^{2}$-regularization.

As such the optimization problem in (2.4) of training a deep learner $f_{\boldsymbol{\theta}}(\cdot)$ involves a highly nonlinear objective function. Stochastic gradient descent (SGD) is a popular tool based on back-propagation (a.k.a. the chain rule), but it often suffers from local traps and overfitting due to the non-convex nature of the problem. We propose data augmentation techniques which can be seamlessly applied in this context and provide efficiency gains. This is achieved via the hierarchical duality between optimization with regularization and finding the maximum a posteriori (MAP) estimate (Polson and Scott, 2011), as described in the following lemma.

Lemma 2.1**.**

The regularization problem

[TABLE]

is equivalent to finding the the Bayesian MAP estimator defined by

[TABLE]

which corresponds to the mode of a posterior distribution characterized as

[TABLE]

Here $p(\boldsymbol{\theta})$ can be interpreted as a prior probability distribution and the log-prior as the regularization penalty.

2.2 A Stochastic Top Layer

By exploiting the duality from Lemma 2.1, we wish to use a Bayesian framework to add stochastic layers – so as to fully account for the uncertainty in estimating the predictive rule $f_{\boldsymbol{\theta}}(\cdot)$. Thus, we convert the sequence of composite functions in the deep learner specified in (2.2) to a stochastic version given by

[TABLE]

Now the hidden variables $Z=(Z_{1},\ldots,Z_{L})^{\prime}$ can be viewed as data augmentation variables and hence will also allow the contribution of fast scalable algorithms for inference and prediction.

For the ease of computation, we only replace the top layer of the DNN with a stochastic layer. We denote network structure below the top layer with $\boldsymbol{B}=\{(W_{1},b_{1}),\ldots,(W_{L},b_{L})\}$, and the network structure can be rewritten as

[TABLE]

where the function $f_{0}(Z_{1}W_{0}+b_{0})$ is the top layer structure and function $f_{\boldsymbol{B}}({\boldsymbol{x}})$ is the network architecture below the top layer. Considering the objective function in (2.4), we implement the solutions with a two-step iterative search. At iteration $t$, we have

DA-update for the top layer ${W_{0},b_{0}}$ as the MAP estimator of the distribution

[TABLE] 2. 2.

SGD-update for the deep architecture $\boldsymbol{B}$

[TABLE] 3. 3.

Sample $Z_{1}^{(t+1)}$ from a normal distribution $\mathcal{N}(\mu_{z}^{(t)},\sigma_{z}^{(t)})$ where $\mu_{z}^{(t)}$ and $\sigma_{z}^{(t)}$ are determined jointly by $\{\theta^{(t)},\boldsymbol{x},\boldsymbol{y}\}$.

The main contribution of our work comes from two aspects: (1) we update top layer weights $\{W_{0},b_{0}\}$ conditional on $\boldsymbol{B}$ as in (2.6), which is also equivalent to conditioning on $Z_{1}$, with data augmentation techniques as later shown in Section 3; (2) the latent variables $Z_{1}$ is sampled from a normal distribution rather than optimized by gradient descent methods. $Z_{1}$ serves as a bridge that connects a weighted $L^{2}$-norm model $f_{0}$ and a deep learner $f_{\boldsymbol{B}}$. Commonly used activation functions $\{f_{l}\}_{l=1}^{L}$ are linear affine functions, rectified linear units (ReLU), sigmoid, hyperbolic tangent (tanh), and etc. We illustrate our methods with a deep ReLU network, i.e., $\{f_{l}\}_{l=1}^{L}$ are ReLU functions, due to its expressibility and inherent sparsity. In the next section, we introduce our data augmentation strategies and show how the stochastic layers can be achieved via data augmentation.

3 Data Augmentation for Deep Learning

Data augmentation introduces a vector of auxiliary variables, denoted by $\boldsymbol{\omega}=(\omega_{1},\ldots,\omega_{n})^{\prime}$ with $\omega_{i}\in\mathbb{R}$, such that the posterior can be written as

[TABLE]

where the augmented auxiliary distribution, $p(\boldsymbol{\theta},\boldsymbol{\omega}\mid{\boldsymbol{y}})$ factorizes nicely into complete conditionals $p(\boldsymbol{\theta}\mid\boldsymbol{\omega},{\boldsymbol{y}})$ and $p(\boldsymbol{\omega}\mid\boldsymbol{\theta},{\boldsymbol{y}})$. A crucial ingredient is that $p(\boldsymbol{\theta}\mid\boldsymbol{\omega},{\boldsymbol{y}})$ is easily managed typically via conditional Gaussians.

Data augmentation tricks allow us to express the likelihood as an expectation of a weighted $L^{2}$-norm. Specifically, we write

[TABLE]

where $p(\boldsymbol{\omega})$ is the prior on the auxiliary variables $\boldsymbol{\omega}=(\omega_{1},\ldots,\omega_{n})^{\prime}$ and the function Q\big{(}\boldsymbol{y}\mid f_{\boldsymbol{\theta}}(\boldsymbol{x}),\boldsymbol{\omega}\big{)} is designed to be a quadratic form, given the data augmentation variables $\boldsymbol{\omega}$. The function $f_{\boldsymbol{\theta}}(\boldsymbol{x})=(f_{0}\circ\cdots\circ f_{L})(\boldsymbol{x})$ is a deep learner.

Table 1 shows that standard activation functions such as ReLU, logit, lasso and check can be expressed in the form of (3.1). Commonly used activation functions for deep learning, with an appropriate stochastic assumptions for $w$ (for notation of simplicity, we derive the standard form for the single observation case) can be expressed as

[TABLE]

Here $\mathcal{GIG}$ denotes the Generalized Inverse Gaussian distribution, ${\mathcal{PG}}$ represents the Pólya Gamma distribution (Polson et al., 2013), and ${\mathcal{E}}$ represents the exponential distribution.

Using the data augmentation strategies, the objectives are represented as mixtures of Gaussians. DA can perform such an optimization with only the use of a sequence of iteratively re-weighted $L^{2}$-norms. This allows us to use XLA techniques to accelerate the training.

Remark 3.1**.**

The log-posterior is optimized given the training data, $\{y_{i},\boldsymbol{x}_{i}\}_{i=1}^{n}$. Deep learning possesses the key property that $\nabla_{\boldsymbol{\theta}}\log p({\boldsymbol{y}}|\theta,{\boldsymbol{x}})$ is computationally inexpensive to evaluate using tensor methods for very complicated architectures and fast implementation on large datasets. One caveat is that the posterior is highly multi-modal and providing good hyperparameter tuning can be expensive. This is clearly a fruitful area of research for state-of-the-art stochastic MCMC algorithms to provide more efficient algorithms. For shallow architectures, the alternating direction method of multipliers (ADMM) is an efficient solution to the optimization problem.

Similarly we can represent the regularization penalty $\exp(-\lambda\phi(\boldsymbol{\theta}))$ in data augmentation form. Hence, we can then replace the optimization problem in (2.4) with

[TABLE]

using the duality in Lemma 2.1.

There are two approaches to Monte Carlo optimization which could handle our data augmentation (Geyer, 1996), missing data methods like Expectation-Maximization (EM) algorithms or stochastic search methods like Markov Chain Monte Carlo (MCMC). The first approach is based on a probabilistic approximation of the objective function (3.1) and is less concerned with exploring $\boldsymbol{\Theta}$. The second type is more exploratory which aims to optimize the objective function by visiting the entire range of $\boldsymbol{\Theta}$ and is less tied to the properties of the function.

For EM algorithms, we consider constructing a surrogate optimization problem which has the same solution to (3.1) (Lange et al., 2000). Specifically, we define a new objective function as

[TABLE]

which is a concave function to be maximized. A natural choice of the surrogate function can be constructed using Jensen’s inequality as

[TABLE]

where each $\omega_{i}$ is drawn from conditional distribution $p(\omega_{i}\mid\boldsymbol{\theta})\propto p(\omega_{i},\boldsymbol{\theta})$ and the minorization is satisfied as

[TABLE]

Maximizing G\big{(}\boldsymbol{\theta}\mid\boldsymbol{\theta}^{(t)}\big{)} with respect to $\boldsymbol{\theta}$ drives $H(\boldsymbol{\theta})$ uphill. The ascent property of the EM algorithm relies on the nonnegativity of the Kullback-Leibler divergence of two conditional probability densities (Hunter and Lange, 2004; Lange, 2013a). The EM algorithm enjoys the numerical stability as it steadily increases the likelihood without wildly overshooting or undershooting. It simplifies the optimization problem by (1) avoiding large matrix inversion; (2) linearizing the objective function; (3) separating the variables of the optimization problem (Lange, 2013b). In Section 4.3 we show how Pólya-Gamma augmentation leads to an EM algorithm for logistic regression.

The exploratory alternative to solve (3.1) is stochastic search methods such as MCMC. The data augmentation strategies enable us to sample from the joint posterior

[TABLE]

where the prior is related to the regularization penalty, via p(\boldsymbol{\theta})\propto\exp\big{(}-\sum_{l=0}^{L}\lambda_{l}\phi_{l}(W_{l},b_{l})\big{)}.

Hence, we can provide an MCMC algorithm in the augmented space $(\boldsymbol{\theta},\boldsymbol{\omega})$ and simulate from the joint posterior distribution, denoted by $p(\boldsymbol{\theta},\boldsymbol{\omega}\mid\boldsymbol{y})$, namely

[TABLE]

A sequence can be simulated using MCMC Gibbs conditionals,

[TABLE]

Then we recover stochastic draws $\boldsymbol{\theta}^{(t)}\sim p(\boldsymbol{\theta}\mid{\boldsymbol{y}})$ from the marginal posterior. These draws can be used in prediction to account for predictive uncertainty, namely

[TABLE]

As $Q(\boldsymbol{y}\mid f_{\boldsymbol{\theta}}(\boldsymbol{x}),\boldsymbol{\omega})$ is conditionally quadratic, the update step for $\boldsymbol{\theta}\mid\boldsymbol{\omega},\boldsymbol{y}$ can be achieved using SGD or a weighted $L^{2}$-norm – the weights $\boldsymbol{\omega}$ are adaptive and provide an automatic choice of the learning rate, thus avoiding backtracking which can be computationally expensive. And the performance of MCMC search is less tied to the statistical properties (i.e. convexity or concavity) of the objective function. We provide examples of how Gaussian regression and SVMs can be implements in Section 4.1 and Section 4.2.

3.1 MCMC with J-copies

The MCMC methods offer a full description of the objective function (3.1) over the entire space $\boldsymbol{\Theta}$. Inspired by the simulated annealing algorithm (Metropolis et al., 1953), we introduce a scaling factor $J$ to allow for faster moves on the surface of (3.1) to maximize. It also helps avoiding the trapping attraction of local maxima. In addition, the corresponding posterior is connected to the Boltzmann distribution, whose density is prescribed by the energy potential $f(\theta)$ and temperature $J$ as

[TABLE]

where $Z_{J}=\int_{\boldsymbol{\Theta}}\exp\left\{-Jf(\boldsymbol{\theta})\right\}d\boldsymbol{\theta}$ is an appropriate normalizing constant.

To simulate the posterior mode without evaluating the likelihood directly (Jacquier et al., 2007), we sample $J$ independent copies of hidden variable $Z_{1}$ . Denoted the copies with $Z_{1}^{1},\ldots,Z_{1}^{J}$, we sample them simultaneously and independently from the posterior distribution

[TABLE]

where $\mu_{z},\sigma_{z}$ are determined by $\{\boldsymbol{x},\boldsymbol{y},\boldsymbol{\theta}\}$. And we stack the J copies as

[TABLE]

where ${\boldsymbol{y}}^{(S)}$, $Z_{1}^{(S)}$ and ${f_{\boldsymbol{B}}}({\boldsymbol{x}}^{(S)})$ are $(n\times J)$-dimensional vectors. We use $Z_{1}^{(S)}$ to amplify the information in ${\boldsymbol{y}}$, which is especially useful in the finite sample problems. Figure 1 illustrates our network architecture.

With the stacked system, the joint distribution of the parameters $\boldsymbol{\theta}$ and the augmented hidden variables $Z_{1}^{(S)}$ given data ${\boldsymbol{y}},{\boldsymbol{x}}$ can be written as

[TABLE]

Hence, the marginal joint posterior

[TABLE]

concentrates on the density proportional to $p({\boldsymbol{x}},{\boldsymbol{y}}\mid\boldsymbol{\theta})^{J}p(\boldsymbol{\theta})$ and provides us with a simulation solution to finding the MAP estimator (Pincus, 1968, 1970).

Another alternative to simulate from the posterior mode is Hamiltonian Monte Carlo (Neal, 2011), which is a modification of Metropolis-Hastings (MH) sampler. Adding an additional momentum variable $\boldsymbol{\nu}$ to the Boltzmann distribution in (3.3), and generating draws from joint distribution

[TABLE]

where $M$ is a mass matrix. Chen et al. (2014) adopt this approach in a deep learning setting.

3.2 Connection to Diffusion Theory

An alternative to the MCMC algorithm can be derived from diffusion theory (Phillips and Smith, 1996). For example, we can approximate the random walk Metropolis-Hastings algorithm with the Langevin diffusion $L_{t}$ defined by the stochastic differential equation $dL_{t}=dB_{t}+\frac{1}{2}\nabla\log f(L_{t})dt$, where $B_{t}$ is the standard Brownian motion. More specifically, let $d:=\left|\boldsymbol{\theta}\right|$, we write the random walk like transition as

[TABLE]

where ${\boldsymbol{\epsilon}}_{t}\sim\mathcal{N}_{d}(0,I_{d})$ and $\sigma^{2}$ corresponds to the discretization size.

This can also be derived by taking a second-order approximation of $\log(f)$, namely

[TABLE]

where $H(\boldsymbol{\theta}^{(t)})=-\nabla^{2}\log f(\boldsymbol{\theta}^{(t)})$ is the Hessian matrix. By taking exponential transformation on both sides, the random walk type approximation to $f(\boldsymbol{\theta}^{(t+1)})$ is

[TABLE]

where $\widetilde{\boldsymbol{\theta}}^{(t)}=\boldsymbol{\theta}^{(t)}+H^{-1}(\boldsymbol{\theta}^{(t)})\nabla\log f(\boldsymbol{\theta}^{(t)})$. If we simplify this approximation by replacing $H(\boldsymbol{\theta}^{(t)})$ with $\sigma^{-2}I_{p}$, the Taylor approximation leads to updating step as

[TABLE]

Roberts and Rosenthal (1998) give further discussion on the choice of $\sigma$ that would yield an acceptance rate of 0.574 to achieve optimal convergence rate.

Mandt et al. (2017) show that SGD can be interpreted as a multivariate Ornstein-Uhlenbeck process

[TABLE]

here $\eta$ is the constant learning rate, $A$ is the symmetric Hessian matrix at the optimum and $\frac{C}{S}$ is the covariance of the mini-batch (of size $S$) gradient noise, which is assumed to be approximately constant near the local optimum of the loss. They also provide results on discrete-time dynamics on other Stochastic Gradient MCMC algorithms, such as Stochastic Gradient Langevin dynamics (SGLD) by Welling and Teh (2011) and Stochastic Gradient Fisher Scoring by Ahn et al. (2012).

Combing their results and the Langevin dynamics of MCMC algorithms, we can write the approximation of our DA-DL updating scheme as

[TABLE]

Similar adaptive dynamics are also observed in other methods. Geman and Hwang (1986) show the convergence of the annealing process using Langevin equations. Slice sampling (Neal, 2003) adaptively chooses the step size based on the local properties of the density function. By constructing local quadratic approximations, it could adapt to the dependencies between variables. Murray et al. (2010) further propose elliptical slice sampling that operates on the ellipse of states.

4 Applications

To illustrate our methodology, we provide three examples: (1) a standard Gaussian regression model with squared loss; (2) a binary classification model under the support vector machine framework; (3) a logistic regression model paired with a Pólya mixing distribution. For the Gaussian regression and SVM models, we implement with $J$-copies stacking strategy to provide full posterior modes.

Before diving into the examples, we introduce the notations we use throughout this section. We continue to denote the output with $\boldsymbol{y}=(y_{1},\ldots,y_{n})^{\prime},y_{i}\in\mathbb{R}$, the input with $\boldsymbol{x}=(\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{n})^{\prime},\boldsymbol{x}_{i}\in\mathbb{R}^{p}$, the latent variable of the top layer with $Z_{1}=(z_{1,1},\ldots,z_{1,n})^{\prime},z_{1,i}\in\mathbb{R}$ and the stacked version as in (3.4). We introduce stochastic noises ${\boldsymbol{\epsilon}}_{0}=(\epsilon_{0,1},\ldots,\epsilon_{0,n})^{\prime}$ in the top layer and ${\boldsymbol{\epsilon}}_{z}=(\epsilon_{z,1},\ldots,\epsilon_{z,n})^{\prime}$ in the second layer, where $\epsilon_{0,i}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,\tau_{0}^{2})$ and $\epsilon_{z,i}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,\tau_{z}^{2})$. The scale parameters $\tau_{0}$ and $\tau_{z}$ are pre-specified and determine the level of randomness or uncertainty for the DA-update and SGD-update respectively. We use $\eta$ to denote the learning rate used in the SGD updates and $T$ is number of training epochs. We use $\left\|\cdot\right\|$ to denote $\ell_{2}$-norm such that $\left\|\boldsymbol{y}\right\|=\sqrt{\sum_{i=1}^{n}y_{i}^{2}}$ and the matrix-type norm as $\left\|\boldsymbol{y}\right\|_{\Sigma}=\sqrt{y^{T}\Sigma y}$.

Our models differ from standard deep learning models and some newly proposed Bayesian approaches in the adoption of stochastic noises ${\boldsymbol{\epsilon}}_{0}$ and ${\boldsymbol{\epsilon}}_{z}$. It distinguishes our model from other deterministic neural networks. By letting ${\boldsymbol{\epsilon}}_{z}$ follow a spiky distribution that puts most of its mass around zero, we can control the estimation approximating to posterior mode instead of posterior mean. The randomness allows us to adopt a stacked system and make the best use of data especially when the dataset is small.

4.1 Gaussian Regression

We consider the regression model as

[TABLE]

The posterior updates are given by

[TABLE]

where $\bar{{\boldsymbol{y}}}=\frac{1}{n}\sum_{i=1}^{n}y_{i}$ and $C_{z}$ is a normalizing constant. The latent variable $Z_{1}$ is drawn from following normal distribution $Z_{1}\sim\mathcal{N}(\mu_{Z},\sigma_{Z}^{2})$ with the mean and variance specified as

[TABLE]

The $J$ copies of $Z_{1}$ are simulated and stacked as

[TABLE]

The updating scheme for this Gaussian regression is summarized in Algorithm 1.

The model can also be generalized to multivariate $y$. Let $y_{i}$ be a $q$-dimension vector, we denote each dimension as $y_{ik},k=1,\ldots,q$, and the model is written as

[TABLE]

where $W_{0}=(W_{01},\ldots,W_{0q})^{\prime}$ is now a q-dimensional vector with $W_{0k}$ computed similarly to (4.1), $b_{0}=(b_{01},\ldots,b_{0q})^{\prime}$ is also q-dimensional with $b_{0k}$ calculated as (4.2). The posterior update for $Z_{1}$ becomes

[TABLE]

which is a multivariate normal distribution with the mean and variance as

[TABLE]

4.2 Support Vector Machines (SVMs)

Support vector machines require data augmentation for rectified linear unit (ReLU) activation functions. Polson and Scott (2011) and Mallick et al. (2005) write the support vector machine model as

[TABLE]

where $p({\boldsymbol{\lambda}})$ follows a flat uniform prior. The augmentation variable ${\boldsymbol{\lambda}}=(\lambda_{1},\ldots,\lambda_{n})^{\prime}$ can be regarded as slacks admitting fuzzy boundaries between classes.

By incorporating the augmentation variable ${\boldsymbol{\lambda}}$, the ReLU deep learning model can be written as

[TABLE]

From a probabilistic perspective, the likelihood function for the output ${\boldsymbol{y}}$ is given by

[TABLE]

Derived from this augmented likelihood function, the conditional updates are

[TABLE]

where $\Lambda=\text{diag}(\lambda_{1},\ldots,\lambda_{n})$ is the diagonal matrix of the augmentation variables.

In order to generate the latent variables, we use conditional Gibbs sampling as

[TABLE]

with the means and variances given by

[TABLE]

where $\mathcal{IG}$ denotes the Inverse Gaussian distribution and $\mathbf{1}=(1,\ldots,1)^{\prime}$ is a n-dimensional unit vector.

The $J$-copies strategy can also be adopted here. $Z_{1}^{j}$ and ${\boldsymbol{\lambda}}^{j}$ needs to be sampled independently for $j=1,\ldots,J$. Algorithm 2 summarizes the updating scheme with $J$-copies for SVMs.

4.3 Logistic Regression

The aim of this example is to show how EM algorithm can be implemented via a weighted $L^{2}$-norm in deep learning. Adopting the logistic regression model from Polson and Scott (2013), we focus on the penalization of $W_{0}$, with parameter optimization given by

[TABLE]

The outcomes $y_{i}$ are coded as $\pm 1$, and $\tau$ is assumed fixed.

For likelihood function $\ell$ and regularization penalty $\phi$, we assume

[TABLE]

where $\mu_{W},\kappa_{W}$ are pre-specified terms controlling the prior of the penalty term and $\lambda$ is endowed with a Pólya distribution prior $P(\lambda)$. Let $\omega_{i}^{-1}$ have a Pólya distribution with $\alpha=1,\kappa=1/2$, the following three updates will generate a sequence of estimates that converges to a stationary point of posterior

[TABLE]

where $z_{i}^{(t)}=y_{i}z_{1,i}^{T}{W_{0}^{(t)}}=y_{i}\text{logit}(\hat{y}_{i}^{t})$, ${\boldsymbol{x}}_{*}$ is a matrix with rows $x_{i}^{*}=y_{i}z_{1,i}$, $\Omega=\text{diag}(\omega_{1},\ldots,\omega_{n})$ and $\Lambda=\text{diag}(\lambda_{1},\ldots,\lambda_{p})$ are diagonal matrices. ${\boldsymbol{x}}_{*}$ can be written as ${\boldsymbol{x}}_{*}=\text{diag}({\boldsymbol{y}})Z_{1}$, $\phi^{\prime}(\cdot)$ denotes the derivative of standard normal density function.

In the non-penalized case, with $\lambda_{i}=0$ for every $i$, the updates can be simplified as weighted least squares

[TABLE]

We focus on the non-penalized binary classification case and Algorithm 3 summarizes our approach. Further generalizations are available. For example, a ridge-regression penalty, along with the generalized double-pareto prior (Armagan et al., 2013) can be implemented by adding a sample-wise $L^{2}$-regularizer. A multinomial generalization of this model can be found in Polson and Scott (2013).

5 Experiments

We illustrate the performance of our methods on both synthetic and real datasets, compared to the deep ReLU networks without the data augmentation layer. We refer to the latter as DL in our results. We denote the data augmented gaussian regression in Algorithm 1 as DA-GR, the SVM implementation in Algorithm 2 as DA-SVM and the logistic regression in Algorithm 3 as DA-logit. For appropriate comparison, we adopt the same network structures, such as the number of layers, the number of hidden nodes, and regularizations like dropout rates, for DL and our methods. The differences between our methods and DL are that (1) the top layer weights ${W_{0},b_{0}}$ of DL are updated via SGD optimization, while the weights ${W_{0},b_{0}}$ of our methods are updated via MCMC or EM; (2) for binary classification, DA-logit and DL adopt a sigmoid activation function in the top layer to produce a binary output, while DA-SVM uses a linear function in the top layer and the augmented sampling layer transforms the continuous value into a binary output. For all experiments, the datasets are partitioned into 70% training and 30% testing randomly. For the optimization we use a modification of the SGD algorithm, the Adaptive moment estimation (Adam, Kingma and Ba (2015)) algorithm. The Adam algorithm combines the estimate of the stochastic gradient with the earlier estimate of the gradient, and scales this using an estimate of the second moment of the unit-level gradient. We have also explored RMSprop (Tieleman and Hinton, 2012) optimizer and we observe similar decreases in regression or classification errors.

To illustrate how the choice of $J$ could affect the speed of convergence, we include different implementations of DA-GR and DA-SVM with $J=2,5,10$. We have explored different sampling noise variance $\tau_{0},\tau_{Z}$, but the choices, in general, do not affect the results significantly.

5.1 Friedman Data

The benchmark (Friedman, 1991) setup uses a regression of the form

[TABLE]

where $\boldsymbol{x}_{i}=(x_{i1},x_{i1},\ldots,x_{ip})$ and only the first 5 covariates are predictive of $y_{i}$. We run the experiments with $n=100,1\,000$ and $p=10,50,100,1\,000$ to explore the performance in both low dimensional and high dimensional scenarios. We implement both one-layer ($L=1$) and two-layer $(L=2)$ ReLU networks with 64 hidden units in each layer. For DA-GR model, we let $\tau_{0}=0.1,\tau_{z}=1$. The experiments are repeated 50 times with different random seeds.

Figure 2 reports the three quartiles of the out-of-sample squared errors (MSEs). The top row is the performance of the one-layer networks and the bottom row is the performance of the two-layer networks. The two-layer networks perform better and converge faster. For DA-GR, when $J=5$ or $J=10$, it converges significantly faster and the prediction errors are also smaller. When $J=2$, the performance of DA-GR is relatively similar to the deep learning model with only SGD updates. This is due to the fact that DA-GR with J-copies learns the posterior mode which is equivalent to the minimization point of the objective function, and it concentrates on the mode faster when $J$ becomes larger.

The computation costs of DA are higher as shown in Figure 3. This is not entirely unexpected since we introduce sampling steps. When $J$ increases, the computation costs also increase slightly. Given the improvement in convergence speed and prediction errors, our data augmentation strategies are still worthwhile even with some extra computation costs. In addition, for each epoch, we can draw the sample-wise posteriors in parallel and the gap between the computation time can be further mitigated.

5.2 Boston Housing Data

Another classical regression benchmark dataset is the Boston Housing dataset111https://archive.ics.uci.edu/ml/machine-learning-databases/housing/, see, for example, Hernández-Lobato and Adams (2015). The data contains $n=506$ observations with 13 features. To show the robustness of DA, we repeat the experiment 20 times with different training subsets. We adopt the ReLU networks with one hidden layer of 64 units and set the dropout rate to be 0.5. For the DA-GR model, we let $\tau_{0}=0.1,\tau_{Z}=1$.

Figure 4 shows the prediction errors of all methods. DA-GR with $J=10$ performs significantly better than the others, in terms of both prediction errors and convergence rates. Meanwhile, DA-GR with $J=2$ behaves similarly to SGD at the beginning, but it converges significantly faster than SGD after a few epochs. This again, shows that with the J-copies strategy, our method helps the optimization converge at a faster speed, and injecting the noise helps the model generalize well out-of-sample.

5.3 Wine Quality Data Set

The Wine Quality Data Set 222P. Cortez, A. Cerdeira, F. Almeida, T. Matos and J. Reis, ‘Wine Quality Data Set’, UCI Machine Learning Repository. contains 4 898 observations with 11 features. The output wine rating is an integer variable ranging from 0 to 10 (the observed range in the data is from 3 to 9). The frequency of each rating is reported in Table 2.

The most frequent ratings are 5 and 6. Since we focus on binary classification problems, we provide two types of classifications, both of which have relatively balanced categories: (1) wine with a rating of 5 or 6 (Test 1); (2) wine with a rating of $\leq 5$ or $>5$ (Test 2). We use the same network architectures adopted in Friedman’s example with $\tau_{0}=\tau_{z}=0.1$.

Figure 5 provides results for the two types of binary classifications. In both cases, DA-SVM performs better than DA-logit and DL. The advantage of large $J$ is still significant and helps converge especially in the early phase. DA-logit outperforms DL in Test 1 when the network is shallow (L=1), while in other cases performs similarly to DL.

5.4 Airbnb Data Set

The Airbnb Kaggle competition333https://www.kaggle.com/c/airbnb-recruiting-new-user-bookings provides a more challenging application with 21 3451 observations in total, and classified by destination into 12 classes: 10 most popular countries, other and no destination found (NDF), where other corresponds to any other country which is not among the top 10 and NDF corresponds to situations that no booking was made. The countries are denoted with their standard codes, as ‘AU’ for Australia, ‘CA’ for Canada, ‘DE’ for Germany, ‘ES’ for Spain, ‘FR’ for France, ‘UK’ for United Kingdom, ‘IT’ for Italy, ‘NL’ for Netherlands, ‘PT’ for Portugal, ‘US’ for United States. Table 3 reports the percentage of each class. We follow the preprocessing steps in Polson and Sokolov (2017). The list of variables contains information from the sessions records (number of sessions, summary statistics of action types, device types and session duration), and user tables such as gender, language, affiliate provider etc. All categorical variables are converted to binary dummies, which leads to 661 features in total. For the neural network architecture, we use a two-layer ReLU network with 64 hidden units on each layer and set the dropout rate to be 0.3. For the SVM model, we let $\tau_{0}=\tau_{z}=0.1$.

Our goal is to test the binary classification models on this dataset. We consider two types of binary responses, both of which have relatively balanced amounts of observations in each category.

Spain (1.05%) vs United Kingdom(1.09%) 2. 2.

United Kingdom (1.09%) vs Italy (1.33%)

Figure 6 demonstrates the binary classifications for Spain versus UK and UK versus Italy. For both cases, the out-of-sample misclassification rates are not small and the fluctuations over epochs are big, suggesting that a better model structure may be needed. However, we still observe that DA-SVM with $J=5$ or $J=10$ has smaller classification errors over epochs and the out-of-sample errors decrease faster during earlier phase of training.

5.5 Summary of Experiment Results

From the above examples, we observe that DA-logit which is implemented under the EM principle does not show an obvious advantage over the vanilla neural network. It shows some improvements on the convergence speed when the network is shallow in the Wine Quality dataset case as in Figure 5. This could be partially due to the fact that we did not apply regularization on the DA layer for our logit implementation. More importantly, the performance of the EM algorithm is contingent on the statistical properties of the objective function. Although the surrogate function is constructed via only the top layer whose quadratic form ensures concavity, the property of the objective function as a whole becomes complicated when the deep network architecture is more complex. Since our method also inherits the negative side of EM and MM algorithms, convergence to the global maximum is not guaranteed in the absence of concavity. However, this observation could open the possibility of future research where we can combine the EM algorithms with shape-constrained neural networks (Gupta et al., 2020).

On the contrary, the MCMC methods with the J-copies strategy significantly improve the prediction errors and convergence speed of the neural networks for both regression and classification problems. And the advantages become more outstanding when $J$ is larger. The phenomenon suggests that the stochastic exploratory methods are preferable when the statistical property of the objective function is unknown or too complex. And the $J$-copies scheme largely relieves the problem of being trapped into local modes.

One concern of using MCMC methods is the extra computation costs induced by the sampling steps. In our current version where $p_{1}=1$, the sample-wise sampling steps can be computed in parallel. If one wishes to introduce a higher dimension latent variable $Z_{1}$ such that $p_{1}>1$, the computation costs will increase as it may involve sampling from multivariate distributions. In that case, fast sampling implementation such as Bhattacharya et al. (2016) is recommended to speed up the process.

6 Discussion

Various regularization methods have been deployed in neural networks to prevent overfitting, such as early stopping, weight decay, dropout (Hinton et al., 2012), gradient noise (Neelakantan et al., 2017). Bayesian strategies tackle the regularization problem by proposing probability structures on the weights. We show that data augmentation strategies are available for many standard activation functions (ReLU, SVM, logit) used in deep learning.

Using MCMC provides a natural stochastic search mechanism that avoids procedures such as back-tracking and provides full descriptions of the objective function over the entire range $\boldsymbol{\Theta}$. Training deep neural networks thus benefits from additional hidden stochastic augmentation units (a.k.a. data augmentation). Uncertainty can be injected into the network through the probabilistic distributions on only one or two layers, permitting more variability of the network. When more data are observed, the level of uncertainty decreases as more information is learned and the network becomes more deterministic. We also exploit the duality between maximum a posteriori estimation and optimization. We provide a $J$-copies stacking scheme to speed up the convergence to posterior mode and avoid trapping attraction of the local modes. Concerning efficiency, DA provides a natural framework to convert the objective function into weighted least squares and is straightforward to implement with the current deep learning training process.

Our three motivational examples illustrated the advantages of data augmentation. Our work has the potential to be generalized to many other data augmentation schemes and different regularization priors. Probabilistic structures on more units and layers are also possible to allow for more uncertainty.

Our DA-DL methods enjoy the benefits of both worlds. On one hand, with the data augmentation on top, it is robust to random weight initialization. Although we still need to specify the learning rates for the deep architecture, the top layer can learn adaptively and the entire network becomes less sensitive to the choice of learning rate. On the other hand, the fast SGD updates from the deep architecture largely alleviate the computation concerns compared to a fully Bayesian hierarchical model.

There are many directions to future research, including adding more sampling layers so the model could accommodate more randomness and flexibility, and using weighted Bayesian bootstrap (Newton et al., 2021) to approximate the unweighted posteriors by assigning random weight to each observation and penalty. Uncertainty quantification for prediction is also possible. Although we focus on the training aspect of deep learning, one can collect posterior draws $\boldsymbol{\theta}^{(t)}$ from the MCMC procedure when the training process converges. Using (3.2), we can construct predictive intervals and conduct inference.