Approximate Empirical Bayes for Deep Neural Networks

Han Zhao∗, Yao-Hung Hubert Tsai∗, Ruslan Salakhutdinov, Geoff Gordon{hanzhao, yaohungt, rsalakhu, ggordon}@cs.cmu.edu
Machine Learning Department
Carnegie Mellon University

Abstract
We propose an approximate empirical Bayes
framework and an efficient algorithm for learning the weight matrix of deep neural networks.
Empirically, we show the proposed method
works as a regularization approach that helps
generalization when training neural networks
on small datasets.

1

Introduction

The empirical Bayes methods provide us a powerful tool
to obtain Bayesian estimators even if we do not havecomplete information about the prior distribution. The
literature on the empirical Bayes methods and its applications is abundant [2, 6, 8, 9, 10, 15, 21, 23]. Existingstudies on parametric empirical Bayes methods focus on
the setting where the likelihood function and the prior are
assumed to have specific forms, e.g., exponential family
distribution and its conjugate prior, so that the marginal
distribution of data has a closed form from which an estimator of the hyperparameter in the prior distribution canbe obtained. While such assumption helps to simplify thesetting in order to demonstrate the power of the empiricalmethod, it restricts us from using more expressive andrich models.
Motivated by the success of the empirical Bayes methodin the Gaussian settings, in this paper we explore extending it to expressive nonlinear models using deep neuralnetworks. Although deep neural networks have been
widely applied in various domains [14, 16, 18, 19], its
parameters are learned via the principle of maximum likelihood in both classification and regression settings, hence
its success crucially hinges on the availability of large
scale datasets [25]. While different kinds of regularization approaches have been studied and designed for neural networks, e.g., weight decay [17], early stopping [4],∗Equal contribution

Dropout [22] and the more recent DeCov [5] method, in
this paper we propose a regularization approach for theweight matrix in neural networks through the lens of theempirical Bayes method. We aim to address the problem of overfitting when training large networks on small
dataset. Our key insight stems from the famous argument
by Efron [6]: It is beneficial to learn from the experienceof others. Specifically, from an algorithmic perspective,
we argue that the connection weights of neurons in thesame layer (row/column vectors of the weight matrix) will
be correlated with each other through the backpropagation learning, and hence by learning from other neurons inthe same layer, a neuron can “borrow statistical strength”
from other neurons to reduce overfitting.
As an illustrating example, consider the simple setting
where the input x ∈ Rd is fully connected to a hidden
layer h ∈ Rp , which is further fully connected to thesingle output ŷ ∈ R. Let σ(·) be the nonlinear activation
function, e.g., ReLU [20], W ∈ Rp×d be the connection
matrix between the input layer and the hidden layer, anda ∈ Rp be the vector connecting the output and the hidden
layer. Without loss of generality ignore the bias term ineach layer, we have: ŷ = aT h, h = σ(W x). Consider
using the usual `2 loss function `(ŷ, y) = 12 |ŷ − y|2 andtake the derivative of `(ŷ, y) w.r.t. W . We obtain the
update formula in backpropagation as follows:
W ← W − α(ŷ − y)(ah0 ) xT ,where h0 is the componentwise derivative of h w.r.t. its
input argument, and α > 0 is the learning rate. For any
fixed a > 0 (or a < 0), realize that (a h0 ) xT is a rank 1
matrix, and furthermore the component of h0 is either 0 or
1. Hence the update for each row vector of W is linearly
proportional to x. Note that the observation holds for any
input pair (x, y), so the update formula implies that therow vectors of W are correlated with each other. Similar
analysis also holds for the column vectors of W . The
above observation leads us to the following question: canwe define a prior distribution over W that captures thecorrelations between its row vectors and column vectors?

2Preliminary

In this section, we first introduce the notation usedthroughout the paper and then give a brief discussion
on the empirical Bayes method [1, 7, 11] and its corresponding estimator [8, 9].
2.1

Notation and Setup

We use lowercase letter, such as y, to represent scalar andlowercase bold letter, such as x, to denote vector. Capitalletter, e.g., X, is reserved for matrix. Calligraphic letter,such as D, is used to denote set. We write tr(A) as thetrace of a matrix A and vec(A) as A’s vectorization (bystacking its column vectors on top of one another). [n]is used to represent the set {1, . . . , n} for any integer n.
a1 a2 stands for the Hadamard product (elementwise
product) between two vectors a1 and a2 .
Suppose we have access to a training set D of n pairs ofdata instances (xi , yi ), i ∈ [n]. We consider the supervised learning setting where xi ∈ X ⊆ Rd and yi ∈ Y.
For a regression problem, Y = R; for a binary classification problem, Y = {1, −1}. Let p(y | x, w) be theconditional distribution of y given x with parameter w.
The parametric form of the conditional distribution is assumed be known. In this paper, we consider a Bayesiansetting that the model parameter w is sampled from aprior distribution p(w | θ) with hyperparameter θ. On
the other hand, given D, the posterior distribution of w
is denoted by p(w | D, θ). From a Bayesian perspective, given an unseen instance x, the goal is to infer thepredictive distribution:Zp(y | x, D, θ) = p(y | x, w) · p(w | D, θ) dw, (1)from which we can compute the mean, or the median,or other statistic (depending on the choice of the loss
function) as our estimator of the unseen target variable y.
2.2

The Empirical Bayes Method

To solve the Bayesian inference in Equation (1), we needaccess to the value of the hyperparameter θ. However,complete information about the hyperparameter θ is usually not available in practice. Considering this issue, empirical Bayes method [8, 21] proposes to estimate θ fromthe data directly using the marginal distribution:θ̂ = arg max p(D | θ)θZ= arg max p(D | w) · p(w | θ) dw.

(2)θ

Under specific choice of the likelihood function p(x, y |w) and the prior distribution p(w | θ), e.g., exponential

family distribution and its corresponding conjugate prior,
we can solve the integral in (2) in closed form to obtain ananalytic solution of θ̂, which can be subsequently pluggedin (1) to obtain a Bayesian estimator for the model parameter w or to compute the predictive distribution of an
unseen instance x.
At a high level, by learning the hyperparameter θ in theprior distribution directly from data, the empirical Bayes
method provides us a principled and convenient way toobtain an estimator of the model parameter w. In fact,
when both the prior and the likelihood functions are normal, it has been formally shown that the empirical Bayes
estimators, e.g., the James-Stein estimator [15] and theEfron-Morris estimator [9], dominate the classic maximum likelihood estimator (MLE) in terms of quadratic
loss for every choice of the model parameter w. At acolloquial level, the success of empirical Bayes estimators can be attributed to the effect of “learning from theexperience of others” [6], which also makes it a powerful tool in recent studies of multitask learning [26] andmeta-learning [12].

3Empirical Bayes Regularization

3.1

Approximate Empirical Bayes

When the likelihood function p(D | w) is implemented
as a neural network, the marginalization in (2) over model
parameter w cannot be computed exactly. Nevertheless,instead of performing expensive Monte-Carlo simulation,
we can use point estimate of w to approximate the integral
as follows:Zp(D | w) · p(w | θ) dw ≈ p(D | ŵ) · p(ŵ | θ). (3)Note that the above approximation will be accurate if thelikelihood under ŵ dominates the likelihoods under othermodel parameters.
Given an estimate ŵ, by maximizing the R.H.S. of (3),we can obtain θ̂ as an approximator of the maximum
marginal likelihood estimator of θ. As a result, we canuse θ̂ to further refine the estimate ŵ by maximizing theposterior distribution as follows:ŵ ← max p(w | D) = max p(D | w) · p(w | θ̂). (4)ww

The maximizer of (4) can in turn be used to better approximate the integral in (3). Formally, we can define
the following optimization problem that characterizes ourframework of the approximate empirical Bayes (AEB)
method:max max log p(D | w) + log p(w | θ)
w

θ

(5)

It is worth to connect the optimization problem (5) to theclassic maximum a posteriori (MAP) inference and also
discuss their difference: if we drop the inner optimization
over the hyperparameter θ in the prior distribution, thenfor any fixed value θ̂, (5) reduces to MAP with the priordefined by the specific choice of θ̂, and the maximizer ŵ
corresponds to the mode of the posterior distribution givenby θ̂. From this perspective, the optimization problemin (5) actually defines a series of MAP inference problems, and the sequence {ŵj (θ̂j )}j defines a solution path
towards the final approximate empirical Bayes estimator.
On the algorithmic side, the optimization problem (5)also suggests a natural block coordinate ascent algorithm
where we alternatively optimize over w and θ until theconvergence of the objective function. We will detail the
discussion of the algorithm in next section where we givea specific prior distribution over the parameters of neuralnetworks.
3.2

The Model and Algorithms

Inspired by the observation from Sec. 1, we propose todefine a matrix-variate normal distribution [13] over theconnection weight matrix W :W ∼ MN (0p×d , Σr , Σc ),(6)

where Σr ∈ Sp++ and Σc ∈ Sd++ are the row and column
covariance matrices, respectively. The probability density
function p(W | Σr , Σc ) is given by:−1T
exp −tr(Σ−1
r W Σc W )/2p(W | Σr , Σc ) =(2π)pd/2 det(Σr )d/2 det(Σc )p/2
Equivalently, one can understand the matrix-variate normal distribution over W as a multivariate normal distribution with a Kronecker product covariance structure overvec(W ): vec(W ) ∼ N (0p×d , Σc ⊗ Σr ). It is then easy
to check that the marginal prior distributions over the row
and column vectors of the weight matrix W are given by:
Wi: ∼ N (0d , [Σr ]ii · Σc ),W:j ∼ N (0p , [Σc ]jj · Σr )We point out that the Kronecker product structure of thecovariance matrix exactly captures our prior about theconnection matrix W : the fan-in/fan-out weights of neurons in the same layer (row/column vectors of W ) arecorrelated with the same correlation matrix, and they only
differ at the scales (variances). Alternatively, if we treat
learning each row/column vector as a separate task, thenthe distribution (6) corresponds to the common prior overdifferent task vectors [26].
−1Define Ωr := Σ−1r and Ωc := Σc to be the corresponding precision matrices and plug in the prior distribution (6)into the general approximate empirical Bayes framework(5). After routine algebraic simplifications, we reach thefollowing concrete optimization problem that corresponds
to the model we described in Sec. 1:
1
1/2 2
min min|ŷ(W, a) − y|2 + λ||Ω1/2
r W Ωc ||FW,a Ωr ,Ωc2− λ (d log det(Ωr ) + p log det(Ωc ))subject to uIp  Ωr  vIp , uId  Ωc  vId (7)where λ is a constant that only depends on p and d, and
0 < u ≤ v are two constants that guarantee the optimization problem (7) is well-posed.
It is not hard to show that in general the optimization problem (7) is not jointly convex in terms of {a, W, Ωr , Ωc },and this holds even if the activation function is linear
and we do not have the hidden layer. However, for any
fixed a, W , the partial optimization over Ωr and Ωc isbi-convex, and more importantly, there is an efficientalgorithm that finds the optimal Ωr (Ωc ) for any fixed
a, W, Ωc (Ωr ) in O(p2 max{d, p}) time [26]. Given apair of constants u, v, we define the following thresholding function T[u,v] (x):u, x < uT[u,v] (x) = x, u ≤ x ≤ v .
(8)v, x > v
We summarize our block coordinate descent algorithm tosolve (7) in Alg. 1. In each iterative loop, Alg. 1 takes afirst-order algorithm A, e.g., the stochastic gradient descent, to optimize the parameters of the neural network bybackpropagation. It then proceeds to compute the optimal
solutions for Ωr and Ωc using the InverseThresholding asa sub-procedure. Alg. 1 terminates when some convergence conditions on model parameters w are satisfied.

4

Experiments

In Sec. 3 we develop our model and algorithms basedon a simple neural network with one hidden layer and asingle output. However, it is straightforward to extend
our AEB framework to more sophisticated and practicalmodels with various structures. In this section we demonstrate the effect of our AEB method on regularizing deepneural networks. Specifically, we consolidate our analysis on two different tasks. The first one is a multi-classclassification problem where we define a matrix-variate
prior over the weight matrix in the last softmax layer. In
the second experiment, we consider a multitask learning
setting where neural networks are used to learn the shared
representations among all the tasks [3], and we apply our
AEB method on the last layer weight matrix, where eachrow corresponds to a separate task vector.
Multi-class Classification. We consider the MNIST
dataset where the target is to classify an input 28 × 28

Algorithm 1 Block Coordinate Descent for Approximate Empirical Bayes

We compare our AEB method with classic regularization methods in the literature of deep learning, including
weight decay and dropout. To show the effect of regularization, we gradually increase the training set size from 60
to 60,000 for 11 different experiments. For each training
set size, we repeat the experiments for 10 times, and weshow the mean along with its standard deviation as error
bars in Fig. 1a.
It is clear from Fig. 1a that AEB helps generalization by
reducing overfitting, and its regularization effect is particularly beneficial when the training set size is small. Fromthe experiments we can also observe that the variance of
the CNN-AEB model is significantly smaller than all theother methods, even if training set size is small.
Multitask Learning. In the second experiment we
consider a dataset for multitask learning. This data
relates to an inverse dynamics problem for a sevendegree-of-freedom (DOF) SARCOS anthropomorphic
robot arm [24]. The goal of this task is to map froma 21-dimensional input space (7 joint positions, 7 joint
velocities, 7 joint accelerations) to the corresponding 7
joint torques. Hence there are 7 tasks and the inputs areshared among all the tasks. The training set and test set
contain 44,484 and 4,449 examples, respectively.
Inspired by the admissible property of the empirical Bayes

100

0.6

95

0.5

90

0.4

Explained Variancedigit image into {0, . . . , 9} categories. In this experiment,
we would like to show that AEB provides an effective
regularization on the network parameters when the training set size is small. To this end, we use a convolutional
neural network as our baseline model with the following
structure: CONV5×5×1×10 -CONV5×5×10×20 -FC320×50 FC50×10 . The notation CONV5×5×1×10 denotes the convolutional layer with kernel size 5 × 5 from depth 1 to
10; the notation FC320×50 denotes the fully connected
layer with size 320 × 50. We ignore bias term for brevity.
As a comparison, in our AEB model we place a matrixvariate normal prior over the weight matrix of the lastsoftmax layer, and we use Alg. 1 to optimize both themodel weights of the convolutional network and the twocovariance matrices of the weight matrix in the last layer.

Accuracy(0)(0)
Input: Initial value w0 := {a(0) , W (0) }, Ωr and Ωc , first-order optimization algorithm A, constants 0 < u ≤ v.
1: for t = 1, . . . , ∞ until convergence do
(t−1)
(t−1)
2:
Fix Ωr, Ωc, optimize w(t) by backpropagation and algorithm A
(t)
(t−1)
3:
Ωr ← InvThresholding(W (t) Ωc
W (t)T , d, u, v)(t)(t)T (t)4:Ωc ← InvThresholding(WΩr W (t) , p, u, v)
5: end for
6:
7: procedure I NV T HRESHOLDING(∆, m, u, v)
8:Compute SVD: Qdiag(r)QT = SVD(∆)
9:
Threshold r0 ← T[u,v] (m/r) (see (8) for the definition of T[u,v] (·))
10:
return Qdiag(r0 )QT
11: end procedure

85
80
CNN
CNN­WeightDecay
CNN­Dropout
CNN­AEB

75
102

103

Train Size

104

MTL
MTL-AEB

0.3
0.2
0.1
0.0

1

2

3

4
Task index

5

6

7

(a) Different regularization (b) Explained variance on thesarcos dataset.
methods on MNIST.

Figure 1: Best viewed zoomed in and in color.
estimator, we are interested in investigating whether our
AEB estimator can lead to better generalization for multiple regression problems, as empirical Bayes estimators do. To do so, we use a baseline fully-connected
network (MTL) with the following structure: FC21×256 FC256×100 -FC100×7 , where the 7 output units correspondto 7 different tasks. Again, in MTL-AEB we place a priorover the task weight matrix and we optimize for bothmodel parameters and two covariance matrices. We plot
and show the explained variance (between 0 and 1, thelarger the better) of MTL and MTL-AEB in Fig. 1b. The
results in Fig. 1b confirm our hypothesis: by making useof the effect “learning from the experience of ours”, AEB
estimator consistently reduces the test set mean squared
error and improves generalization.

5

Conclusion

We propose an approximate empirical Bayes method forlearning the model parameters of deep neural networks.
On two datasets we demonstrate that our AEB method
helps to reduce overfitting and hence improves generalization when the training set is small. One future direction
is to develop a more efficient approximate solution to
optimize the two covariance matrices so that our method
can scale to larger networks. It is also interesting to see
whether applying our AEB method to all the layers canfurther reduce overfitting or not.

References[1] J. M. Bernardo and A. F. Smith. Bayesian theory,
2001.
[2] B. P. Carlin and T. A. Louis. Bayes and empirical
Bayes methods for data analysis. Chapman andHall/CRC, 2010.
[3] R. Caruana. Multitask learning. In Learning tolearn, pages 95–133. Springer, 1998.
[4] R. Caruana, S. Lawrence, and C. L. Giles. Overfitting in neural nets: Backpropagation, conjugate
gradient, and early stopping. In Advances in neural information processing systems, pages 402–408,
2001.
[5] M. Cogswell, F. Ahmed, R. Girshick, L. Zitnick,and D. Batra. Reducing overfitting in deep networks
by decorrelating representations. arXiv preprint
arXiv:1511.06068, 2015.
[6] B. Efron. Large-scale inference: empirical Bayes
methods for estimation, testing, and prediction, volume 1. Cambridge University Press, 2012.
[7] B. Efron and T. Hastie. Computer age statisticalinference, volume 5. Cambridge University Press,
2016.
[8] B. Efron and C. Morris. Stein’s estimation ruleand its competitors—an empirical Bayes approach.
Journal of the American Statistical Association, 68(341):117–130, 1973.
[9] B. Efron and C. Morris. Stein’s paradox in statistics.
Scientific American, 236(5):119–127, 1977.
[10] B. Efron, R. Tibshirani, J. D. Storey, and V. Tusher.
Empirical bayes analysis of a microarray experiment.
Journal of the American statistical association, 96(456):1151–1160, 2001.
[11] A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson,
A. Vehtari, and D. B. Rubin. Bayesian data analysis.
CRC press, 2013.
[12] E. Grant, C. Finn, S. Levine, T. Darrell, and T. Griffiths. Recasting gradient-based meta-learning as hierarchical bayes. arXiv preprint arXiv:1801.08930,
2018.
[13] A. K. Gupta and D. K. Nagar. Matrix variate distributions. Chapman and Hall/CRC, 2018.
[14] K. He, X. Zhang, S. Ren, and J. Sun. Deep residual
learning for image recognition. In Proceedings ofthe IEEE conference on computer vision and patternrecognition, pages 770–778, 2016.
[15] W. James and C. Stein. Estimation with quadratic
loss. In Proceedings of the fourth Berkeley symposium on mathematical statistics and probability,
volume 1, pages 361–379, 1961.

[16] A. Krizhevsky and G. Hinton. Learning multiple
layers of features from tiny images. 2009.
[17] A. Krogh and J. A. Hertz. A simple weight decay
can improve generalization. In Advances in neural information processing systems, pages 950–957,
1992.
[18] Y. LeCun, B. E. Boser, J. S. Denker, D. Henderson, R. E. Howard, W. E. Hubbard, and L. D.
Jackel. Handwritten digit recognition with a backpropagation network. In Advances in neural information processing systems, pages 396–404, 1990.
[19] Y. LeCun, Y. Bengio, and G. Hinton. Deep learning.
nature, 521(7553):436, 2015.
[20] V. Nair and G. E. Hinton. Rectified linear units
improve restricted boltzmann machines. In Proceedings of the 27th international conference on machinelearning (ICML-10), pages 807–814, 2010.
[21] H. Robbins. An empirical bayes approach to statistics. Technical report, Columbia University, NewYork City, United States, 1956.
[22] N. Srivastava, G. Hinton, A. Krizhevsky,
I. Sutskever, and R. Salakhutdinov. Dropout:
A simple way to prevent neural networks fromoverfitting. The Journal of Machine Learning
Research, 15(1):1929–1958, 2014.
[23] C. Stein. Inadmissibility of the usual estimator forthe mean of a multivariate normal distribution. Technical report, Stanford University, Stanford, UnitedStates, 1956.
[24] S. Vijayakumar and S. Schaal. Locally weighted
projection regression: Incremental real time learning in high dimensional space. In Proceedings of
the Seventeenth International Conference on Machine Learning, pages 1079–1086. Morgan Kaufmann Publishers Inc., 2000.
[25] C. Zhang, S. Bengio, M. Hardt, B. Recht, and
O. Vinyals. Understanding deep learning requires rethinking generalization. arXiv preprint
arXiv:1611.03530, 2016.
[26] H. Zhao, O. Stretcu, R. Negrinho, A. Smola, and
G. Gordon. Efficient multi-task feature and relationship learning. arXiv preprint arXiv:1702.04423,
2017.