# HesScale: Scalable Computation of Hessian
Diagonals

Mohamed Elsayed & A. Rupam Mahmood
Department of Computing Science, Alberta Machine Intelligence Institute (Amii)
University of Alberta
Edmonton, Alberta, Canada
<EMAIL_ADDRESS>
CIFAR AI Chair

###### Abstract

Second-order optimization uses curvature information about the objective
function, which can help in faster convergence. However, such methods
typically require expensive computation of the Hessian matrix, preventing
their usage in a scalable way. The absence of efficient ways of computation
drove the most widely used methods to focus on first-order approximations that
do not capture the curvature information. In this paper, we develop HesScale,
a scalable approach to approximating the diagonal of the Hessian matrix, to
incorporate second-order information in a computationally efficient manner. We
show that HesScale has the same computational complexity as backpropagation.
Our results on supervised classification show that HesScale achieves high
approximation accuracy, allowing for scalable and efficient second-order
optimization.111Code is available at https://github.com/mohmdelsayed/HesScale.

##  1 Introduction

First-order optimization offers a cheap and efficient way of performing local
progress in optimization problems by using gradient information. However,
their performance suffers from instability or slow progress when used in ill-
conditioned landscapes. Such a problem is present because first-order methods
do not capture curvature information which causes two interrelated issues.
First, the updates in first-order have incorrect units (Duchi et al. 2011),
which creates a scaling issue. Second, first-order methods lack
parameterization invariance (Martens 2020) in contrast to second-order methods
such as natural gradient (Amari 1998) or Newton-Raphson methods. Therefore,
some first-order normalization methods were developed to address the
invariance problem (Ba et al. 2016, Ioffe & Szegedy 2015, Salimans & Kingma
2016). On the other hand, some recent adaptive step-size methods try to
alleviate the scaling issue by using gradient information for first-order
curvature approximation (Luo et al. 2019, Duchi et al. 2011, Zeiler 2012,
Reddi et al. 2018, Kingma & Ba 2015, Tran & Phong 2019, Tieleman et al. 2012).
Specifically, such methods use the empirical Fisher diagonals heuristic by
maintaining a moving average of the squared gradients to approximate the
diagonal of the Fisher information matrix. Despite the huge adoption of such
methods due to their scalability, they use inaccurate approximations. Kunstner
et al. (2019) showed that the empirical Fisher does not generally capture
curvature information and might have undesirable effects. They argued that the
empirical Fisher approximates the Fisher or the Hessian matrices only under
strong assumptions that are unlikely to be met in practice. Moreover, Wilson
et al. (2017) presented a counterexample where the adaptive step-size methods
are unable to reduce the error compared to non-adaptive counterparts such as
stochastic gradient descent.

Although second-order optimization can speed up the training process by using
the geometry of the landscape, its adoption is minimal compared to first-order
methods. The exact natural gradient or Newton-Raphson methods require the
computation, storage, and inversion of the Fisher information or the Hessian
matrices, making them computationally prohibitive in large-scale tasks.
Accordingly, many popular second-order methods attempt to approximate less
expensively. For example, a type of truncated-Newton method called Hessian-
free methods (Martens 2010) exploits the fact that the Hessian-vector product
is cheap (Bekas et al. 2007) and uses the iterative conjugate gradient method
to perform an update. However, such methods might require many iterations per
update or some tricks to achieve stability, adding computational overhead
(Martens & Sutskever 2011). Some variations try to approximate only the
diagonals of the Hessian matrix using stochastic estimation with matrix-free
computations (Chapelle & Erhan 2011, Martens et al. 2012, Yao et al. 2021).
Other methods impose probabilistic modeling assumptions and estimate a block
diagonal Fisher information matrix (Martens & Grosse 2015, Botev et al. 2017).
Such methods are invariant to reparametrization but are computationally
expensive since they need to perform matrix inversion for each block.

Deterministic diagonal approximations to the Hessian (LeCun et al. 1990,
Becker & Lecun 1989) provide some curvature information and are efficient to
compute. Specifically, they can be implemented to be as efficient as first-
order methods. We view this category of approximation as scalable second-order
methods. In neural networks, curvature backpropagation (Becker & Lecun 1989,
Mizutani & Dreyfus 2008) can be used to backpropagate the curvature vector.
Although these methods show a promising direction for scalable second-order
optimization, the approximation quality is sometimes poor with objectives such
as cross-entropy (Martens et al. 2012). A scalable second-order method with
high quality approximation is still needed.

In this paper, we present HesScale, a high-quality approximation method for
the Hessian diagonals. Our method is also scalable and has little memory
requirement with linear computational complexity while maintaining high
approximation accuracy.

##  2 Background

In this section, we describe the Hessian matrix for neural networks and some
existing methods for estimating it. Generally, Hessian matrices can be
computed for any scalar-valued function that are twice differentiable. If
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is such a function, then for its
argument ${\bm{\psi}}\in\mathbb{R}^{n}$, the Hessian matrix
${\bm{H}}\in\mathbb{R}^{n\times n}$ of $f$ with respect to ${\bm{\psi}}$ is
given by
$H_{i,j}=\nicefrac{{\partial^{2}f({\bm{\psi}})}}{{\partial\psi_{i}\partial\psi_{j}}}$.
Here, the $i$th element of a vector ${\bm{v}}$ is denoted by $v_{i}$, and the
element at the $i$th row and $j$th column of a matrix ${\bm{M}}$ is denoted by
$M_{i,j}$. When the need for computing the Hessian matrix arises for
optimization in deep learning, the function $f$ is typically the objective
function, and the vector ${\bm{\psi}}$ is commonly the weight vector of a
neural network. Computing and storing an $n\times n$ matrix, where $n$ is the
number of weights in a neural network, is expensive. Therefore, many methods
exist for approximating the Hessian matrix or parts of it with less memory
footprint, computational requirement, or both. A common technique is to
utilize the structure of the function to reduce the computations needed. For
example, assuming that connections from a certain layer do not affect other
layers in a neural network allows one to approximate a block diagonal Hessian.
The computation further simplifies when we have piece-wise linear activation
functions (e.g., ReLU), which result in a _Generalized Gauss-Newton_ (GGN)
(Schraudolph 2002) approximation that is equivalent to the block diagonal
Hessian matrix with linear activation functions. The GGN matrix is more
favored in second-order optimization since it is positive semi-definite.
However, computing a block diagonal matrix is still demanding.

Many approximation methods were developed to reduce the storage and
computation requirements of the GGN matrix. For example, under probabilistic
modeling assumptions, the _Kronecker-factored Approximate Curvature_ (KFAC)
method (Martens & Grosse 2015) writes the GGN matrix ${\bm{G}}$ as a Kronecker
product of two matrices of smaller sizes as:
${\bm{G}}={\bm{A}}\otimes{\bm{B}}$, where
${\bm{A}}=\mathbb{E}[{\bm{h}}{\bm{h}}^{\top}]$,
${\bm{B}}=\mathbb{E}[{\bm{g}}{\bm{g}}^{\top}]$, ${\bm{h}}$ is the activation
output vector, and ${\bm{g}}$ is the gradient of the loss with respect to the
activation input vector. The ${\bm{A}}$ and ${\bm{B}}$ matrices can be
estimated by Monte Carlo sampling and an exponential moving average. KFAC is
more efficient when used in optimization since it requires inverting only the
small matrices using the Kronecker-product property
$({\bm{A}}\otimes{\bm{B}})^{-1}={\bm{A}}^{-1}\otimes{\bm{B}}^{-1}$. However,
KFAC is still expensive due to the storage of the block diagonal matrices and
computation of Kronecker product, which prevent it from being used as a
scalable method.

Computing the Hessian diagonals can provide some curvature information with
relatively less computation. However, it has been shown that the exact
computation for diagonals of the Hessian typically has quadratic complexity
with the unlikely existence of algorithms that can compute the exact diagonals
with less than quadratic complexity (Martens et al. 2012). Some stochastic
methods provide a way to compute unbiased estimates of the exact Hessian
diagonals. For example, the AdaHessian (Yao et al. 2021) algorithm uses the
Hutchinson’s estimator
$\text{diag}({\bm{H}})=E[{\bm{z}}\circ({\bm{H}}{\bm{z}})]$, where ${\bm{z}}$
is a multivariate random variable with a Rademacher distribution and the
expectation can be estimated using Monte Carlo sampling with an exponential
moving average. Similarly, the GGN-MC method (Dangel et al. 2020) uses the
relationship between the Fisher information matrix and the Hessian matrix
under probabilistic modeling assumptions to have an MC approximation of the
diagonal of the GGN matrix. Although these stochastic approximation methods
are scalable due to linear or $O(n)$ computational and memory complexity, they
suffer from low approximation quality, improving which requires many sampling
and factors of additional computations.

##  3 The Proposed HesScale Method

In this section, we present our method for approximating the diagonal of the
Hessian at each layer in feed-forward networks, where a backpropagation rule
is used to utilize the Hessian of previous layers. We present the derivation
of the backpropagation rule for fully connected and convolutional neural
networks in supervised learning. Similar derivation for fully connected
networks with mean squared error is presented before (LeCun et al. 1990,
Becker & Lecun 1989). However, we use the exact diagonals of the Hessian
matrix at the last layer with some non-linear and non-element-wise output
activations such as softmax and show that it can still be computed in linear
computational complexity. We show the derivation for Hessian diagonals for
fully connected networks in the following and provide the derivation for the
convolutional neural networks in Appendix B.

We use the supervised classification setting where there is a collection of
data examples. These data examples are generated from some _target function_
${f}^{*}$ mapping the input ${\bm{x}}$ to the output $y$, where the $k$-th
input-output pair is $({\bm{x}}_{k},y_{k})$. In this task, the _learner_ is
required to predict the output class $y\in\\{1,2,...,m\\}$ given the input
vector ${\bm{x}}\in\mathbb{R}^{d}$ by estimating the target function
${f}^{*}$. The performance is measured with the cross-entropy loss,
$\mathcal{L}({\bm{p}},{\bm{q}})=-\sum_{i=1}^{m}p_{i}\log q_{i}$, where
${\bm{p}}\in\mathbb{R}^{m}$ is the vector of the target one-hot encoded class
and ${\bm{q}}\in\mathbb{R}^{m}$ is the predicted output. The learner is
required to reduce the cross-entropy by matching the target class.

Consider a neural network with $L$ layers that outputs the predicted output
${\bm{q}}$. The neural network is parametrized by the set of weights
$\\{{\bm{W}}_{1},...,{\bm{W}}_{L}\\}$, where ${\bm{W}}_{l}$ is the weight
matrix at the $l$-th layer, and its element at the $i$th row and the $j$th
column is denoted by $W_{l,i,j}$. During learning, the parameters of the
neural network are changed to reduce the loss. At each layer $l$, we get the
activation output ${\bm{h}}_{l}$ by applying the activation function
$\bm{\sigma}$ to the activation input ${\bm{a}}_{l}$:
${\bm{h}}_{l}=\bm{\sigma}({\bm{a}}_{l})$. We simplify notations by defining
${\bm{h}}_{0}\doteq{\bm{x}}$. The activation output ${\bm{h}}_{l}$ is then
multiplied by the weight matrix ${\bm{W}}_{l+1}$ of layer $l+1$ to produce the
next activation input:
${a}_{l+1,i}=\sum_{j=1}^{|{\bm{h}}_{l}|}{{W}_{l+1,i,j}}{h}_{l,j}$. We assume
here that the activation function is element-wise activation for all layers
except for the final layer $L$, where it becomes the softmax function. The
backpropagation equations for the described network are given as follows
Rumelhart et al. (1986):

$\displaystyle\frac{\partial\mathcal{L}}{\partial a_{l,i}}$
$\displaystyle=\sum_{k=1}^{|{\bm{a}}_{l+1}|}\frac{\partial\mathcal{L}}{\partial
a_{l+1,k}}\frac{\partial a_{l+1,k}}{\partial h_{l,i}}\frac{\partial
h_{l,i}}{\partial
a_{l,i}}=\sigma^{\prime}(a_{l,i})\sum_{k=1}^{|{\bm{a}}_{l+1}|}\frac{\partial\mathcal{L}}{\partial
a_{l+1,k}}W_{l+1,k,i},$ (1) $\displaystyle\frac{\partial\mathcal{L}}{\partial
W_{l,i,j}}$ $\displaystyle=\frac{\partial\mathcal{L}}{\partial
a_{l,i}}\frac{\partial a_{l,i}}{\partial
W_{l,i,j}}=\frac{\partial\mathcal{L}}{\partial a_{l,i}}h_{l-1,j}.$ (2)

In the following, we write the equations for the exact Hessian diagonals with
respect to weights
$\nicefrac{{\partial^{2}\mathcal{L}}}{{\partial{W^{2}_{l,i,j}}}}$, which
requires the calculation of
$\nicefrac{{\partial^{2}\mathcal{L}}}{{\partial{a^{2}_{l,i}}}}$ first:

$\displaystyle\frac{\partial^{2}\mathcal{L}}{\partial{a^{2}_{l,i}}}$
$\displaystyle=\frac{\partial}{\partial
a_{l,i}}\left(\sigma^{\prime}(a_{l,i})\sum_{k=1}^{|{\bm{a}}_{l+1}|}\frac{\partial\mathcal{L}}{\partial
a_{l+1,k}}W_{l+1,k,i}\right)$
$\displaystyle={\sigma^{\prime}}(a_{l,i})\sum_{k=1}^{|{\bm{a}}_{l+1}|}\sum_{p=1}^{|{\bm{a}}_{l+1}|}\frac{\partial^{2}\mathcal{L}}{\partial
a_{l+1,k}\partial a_{l+1,p}}\frac{\partial a_{l+1,p}}{\partial
a_{l,i}}W_{l+1,k,i}+\sigma^{\prime\prime}(a_{l,i})\sum_{k=1}^{|{\bm{a}}_{l+1}|}\frac{\partial\mathcal{L}}{\partial
a_{l+1,k}}W_{l+1,k,i}$
$\displaystyle={\sigma^{\prime}}(a_{l,i})^{2}\sum_{k=1}^{|{\bm{a}}_{l+1}|}\sum_{p=1}^{|{\bm{a}}_{l+1}|}\frac{\partial^{2}\mathcal{L}}{\partial
a_{l+1,k}\partial
a_{l+1,p}}W_{l+1,p,i}W_{l+1,k,i}+\sigma^{\prime\prime}(a_{l,i})\sum_{k=1}^{|{\bm{a}}_{l+1}|}\frac{\partial\mathcal{L}}{\partial
a_{l+1,k}}W_{l+1,k,i},$
$\displaystyle\frac{\partial^{2}\mathcal{L}}{\partial{W^{2}_{l,i,j}}}$
$\displaystyle=\frac{\partial}{\partial
W_{l,i,j}}\left(\frac{\partial\mathcal{L}}{\partial
a_{l,i}}h_{l-1,j}\right)=\frac{\partial}{\partial
a_{l,i}}\left(\frac{\partial\mathcal{L}}{\partial
a_{l,i}}\right)\frac{\partial a_{l,i}}{\partial
W_{l,i,j}}h_{l-1,j}=\frac{\partial^{2}\mathcal{L}}{\partial
a^{2}_{l,i}}h^{2}_{l-1,j}.$ (3)

Since, the calculation of
$\nicefrac{{\partial^{2}\mathcal{L}}}{{\partial{a^{2}_{l,i}}}}$ depends on the
off-diagonal terms, the computation complexity becomes quadratic. Following
Becker and Lecun (1989), we approximate the Hessian diagonals by ignoring the
off-diagonal terms, which leads to a backpropagation rule with linear
computational complexity for our estimates
$\widehat{\frac{\partial^{2}\mathcal{L}}{\partial{W^{2}_{l,i,j}}}}$ and
$\widehat{\frac{\partial^{2}\mathcal{L}}{\partial{a^{2}_{l,i}}}}$:

$\displaystyle\widehat{\frac{\partial^{2}\mathcal{L}}{\partial{a^{2}_{l,i}}}}$
$\displaystyle\doteq{\sigma^{\prime}}(a_{l,i})^{2}\sum_{k=1}^{|{\bm{a}}_{l+1}|}\widehat{\frac{\partial^{2}\mathcal{L}}{\partial{a^{2}_{l+1,k}}}}W^{2}_{l+1,k,i}+\sigma^{\prime\prime}(a_{l,i})\sum_{k=1}^{|{\bm{a}}_{l+1}|}\frac{\partial\mathcal{L}}{\partial
a_{l+1,k}}W_{l+1,k,i},$ (4)
$\displaystyle\widehat{\frac{\partial^{2}\mathcal{L}}{\partial{W^{2}_{l,i,j}}}}$
$\displaystyle\doteq\widehat{\frac{\partial^{2}\mathcal{L}}{\partial{a^{2}_{l,i}}}}h^{2}_{l-1,j}.$
(5)

However, for the last layer, we use the exact Hessian diagonals
$\widehat{\frac{\partial^{2}\mathcal{L}}{\partial{a^{2}_{L,i}}}}\doteq\frac{\partial^{2}\mathcal{L}}{\partial{a^{2}_{L,i}}}$
since it can be computed in $O(n)$ for the softmax activation function and the
cross-entropy loss. More precisely, the exact Hessian diagonals for cross-
entropy loss with softmax is simply ${\bm{q}}-{\bm{q}}\circ{\bm{q}}$, where
${\bm{q}}$ is the predicted probability vector and $\circ$ denotes element-
wise multiplication. We found empirically that this small change makes a large
difference in the approximation quality, as shown in Fig. 1(a). Hence, unlike
Becker and Lecun (1989) who use a Hessian diagonal approximation of the last
layer by Eq. 4, we use the exact values directly to achieve more approximation
accuracy. We call this method for Hessian diagonal approximation _HesScale_
and provide its pseudocode for supervised classification in Algorithm 1.

Algorithm 1 HesScale: Computing Hessian diagonals of a neural network layer in
classification

Neural network $f$ and a layer number $l$

First and second order information
$\widehat{\frac{\partial\mathcal{L}}{\partial a_{l+1,i}}}$ and
$\widehat{\frac{\partial^{2}\mathcal{L}}{\partial a_{l+1,i,j}^{2}}}$, unless
$l=L$

Input-output pair $({\bm{x}},y)$

Set loss function $\mathcal{L}$ to cross-entropy loss

Compute preference vector ${\bm{a}}_{L}\leftarrow f({\bm{x}})$ and target one-
hot-encoded vector ${\bm{p}}\leftarrow\texttt{onehot}(y)$

Compute the predicted probability vector
${\bm{q}}\leftarrow{\bm{\sigma}}({\bm{a}}_{L})$ using softmax function
${\bm{\sigma}}$

Compute the error $\mathcal{L}({\bm{p}},{\bm{q}})$

if $l=L$ then $\triangleright$ Computing Hessian diagonals exactly for the
last layer

Compute
$\frac{\partial\mathcal{L}}{\partial{\bm{a}}_{L}}\leftarrow{\bm{q}}-{\bm{p}}$$\triangleright$
$\frac{\partial\mathcal{L}}{\partial{\bm{a}}_{L}}$ consists of elements
$\frac{\partial\mathcal{L}}{\partial a_{L,i}}$

Compute $\frac{\partial\mathcal{L}}{\partial{\bm{W}}_{L}}$ using Eq. 2
$\triangleright$ $\frac{\partial\mathcal{L}}{\partial{\bm{W}}_{L}}$ consists
of elements $\frac{\partial\mathcal{L}}{\partial W_{L,i,j}}$

$\widehat{\frac{\partial^{2}\mathcal{L}}{\partial{\bm{a}}_{L}^{2}}}\leftarrow{\bm{q}}-{\bm{q}}\circ{\bm{q}}$
$\triangleright$
$\widehat{\frac{\partial^{2}\mathcal{L}}{\partial{\bm{a}}_{L}^{2}}}$ consists
of elements $\widehat{\frac{\partial^{2}\mathcal{L}}{\partial a_{L,i}^{2}}}$

Compute $\widehat{\frac{\partial^{2}\mathcal{L}}{\partial{\bm{W}}_{L}^{2}}}$
using Eq. 5$\triangleright$
$\widehat{\frac{\partial^{2}\mathcal{L}}{\partial{\bm{W}}_{L}^{2}}}$ consists
of elements $\widehat{\frac{\partial^{2}\mathcal{L}}{\partial W_{L,i,j}^{2}}}$

else if $l\neq L$ then

Compute $\frac{\partial\mathcal{L}}{\partial{\bm{a}}_{l}}$ and
$\nicefrac{{\partial\mathcal{L}}}{{\partial{\bm{W}}_{l}}}$ using Eq. 1 and Eq.
2

Compute $\widehat{\frac{\partial^{2}\mathcal{L}}{\partial{\bm{a}}_{l}^{2}}}$
and $\widehat{\frac{\partial^{2}\mathcal{L}}{\partial{\bm{W}}_{l}^{2}}}$ using
Eq. 4 and Eq. 5

end if

return $\frac{\partial\mathcal{L}}{\partial{\bm{W}}_{l}}$,
$\widehat{\frac{\partial^{2}\mathcal{L}}{\partial{\bm{W}}_{l}^{2}}}$,
$\frac{\partial\mathcal{L}}{\partial{\bm{a}}_{l}}$, and
$\widehat{\frac{\partial^{2}\mathcal{L}}{\partial{\bm{a}}_{l}^{2}}}$

HesScale is not specific to cross-entropy loss as the exact Hessian diagonals
can be calculated in $O(n)$ for some other widely used loss functions as well.
We show this property for negative log-likelihood function with Gaussian and
softmax distributions in Appendix A. The computations can be reduced further
using a linear approximation for the activation functions (by dropping the
second term in Eq. 4), which corresponds to an approximation of the GGN
matrix. We call this variation of our method _HesScaleGN_.

Based on HesScale, we make a stable optimizer, which we call AdaHesScale,
given in Algorithm 2. We use the same style introduced in Adam (Kingma & Ba
2015), using the squared diagonal approximation instead of the squared
gradients to update the moving average. Moreover, we introduce another
optimizer based on HesScaleGN, which we call AdaHesScaleGN.

Algorithm 2 AdaHesScale for optimization

Neural network $f$ with weights $\\{{\bm{W}}_{1},...,{\bm{W}}_{L}\\}$ and a
dataset $\mathcal{D}$

Small number $\epsilon\leftarrow 10^{-8}$

Exponential decay rates $\beta_{1},\beta_{2}\in[0,1)$

step size $\alpha$

Initialize $\\{{\bm{W}}_{1},...,{\bm{W}}_{L}\\}$

Initialize time step $t\leftarrow 0$.

for $l$ in $\\{L,L-1,...,1\\}$ do $\triangleright$ Set exponential moving
averages at time step 0 to zero

${\bm{M}}_{l}\leftarrow 0;\quad{\bm{V}}_{l}\leftarrow 0$ $\triangleright$ Same
size as ${\bm{W}}_{l}$

end for

for $({\bm{x}},y)$ in $\mathcal{D}$ do

$t\leftarrow t+1$

${\bm{r}}_{L+1}\leftarrow{\bm{s}}_{L+1}\leftarrow\bm{\emptyset}$
$\triangleright$ ${\bm{r}}_{l}$ and ${\bm{s}}_{l}$ stand for
$\frac{\partial\mathcal{L}}{\partial{\bm{a}}_{l}}$ and
$\widehat{\frac{\partial^{2}\mathcal{L}}{\partial{\bm{a}}_{l}^{2}}}$,
respectively

for $l$ in $\\{L,L-1,...,1\\}$ do

${\bm{F}}_{l},{\bm{S}}_{l},{\bm{r}}_{l},{\bm{s}}_{l}\;\leftarrow\;$
HesScale($f,{\bm{x}},y,l,{\bm{r}}_{l+1},{\bm{s}}_{l+1}$). $\triangleright$
Check Algorithm 1

${\bm{M}}_{l}\leftarrow\beta_{1}{\bm{M}}_{l}+(1-\beta_{1}){\bm{F}}_{l}$
$\triangleright$ ${\bm{F}}_{l}$ stands for
$\frac{\partial\mathcal{L}}{\partial{\bm{W}}_{l}}$

${\bm{V}}_{l}\leftarrow\beta_{2}{\bm{V}}_{l}+(1-\beta_{2}){\bm{S}}_{l}^{2}$
$\triangleright$ ${\bm{S}}_{l}$ stands for
$\widehat{\frac{\partial^{2}\mathcal{L}}{\partial{\bm{W}}_{l}^{2}}}$

$\hat{{\bm{M}}}_{l}\leftarrow{\bm{M}}_{l}/(1-\beta_{1}^{t})$ $\triangleright$
Bias-corrected estimate for ${\bm{F}}_{l}$

$\hat{{\bm{V}}}_{l}\leftarrow{\bm{V}}_{l}/(1-\beta_{2}^{t})$ $\triangleright$
Bias-corrected estimate for ${\bm{S}}_{l}$

${\bm{W}}_{l}\leftarrow{\bm{W}}_{l}-\alpha\hat{{\bm{M}}}_{l}\oslash(\hat{{\bm{V}}}_{l}+\epsilon)^{\circ\frac{1}{2}}$
$\triangleright$ $\oslash$ is element-wise division

$\triangleright$ ${\bm{A}}^{\circ\frac{1}{2}}$ is element-wise square root of
${\bm{A}}$

end for

end for

##  4 Approximation quality & scalability of HesScale

In this section, we evaluate HesScale for its approximation quality and
computational cost and compare it with other methods. These measures
constitute the criteria we look for in scalable and efficient methods. For our
experiments, we implemented HesScale using the _BackPack_ framework (Dangel et
al. 2020), which allows easy implementation of backpropagation of statistics
other than the gradient.

We start by studying the approximation quality of Hessian diagonals compared
to the true values. To measure the approximation quality of the Hessian
diagonals for different methods, we use the $L^{1}$ distance between the exact
Hessian diagonals and their approximations. Our task here is supervised
classification, and data examples are randomly generated. We used a network of
three hidden layers with _tanh_ activations, each containing 16 units. The
network weights and biases are initialized randomly. The network has six
inputs and ten outputs. For each example pair, we compute the exact Hessian
diagonals for each layer and their approximations from each method. All
layers’ errors are summed and averaged over 1000 data examples for each
method. In this experiment, we used 40 different initializations for the
network weights, shown as colored dots in Fig. 1(a). Each point represents the
summed error over network layers, averaged over 1000 examples for each
different initialization. In this figure, we show the average error incurred
by each method normalized by the average error incurred by HesScale. Any
approximation that incurs an averaged error above 1 has a worse approximation
than HesScale, and any approximation with an error less than 1 has a better
approximation than HesScale. Moreover, we show the layer-wise error for each
method in Fig. 1(b).

(a) Normalized $L^{1}$ error with respect to HesScale

(b) Layer-wise $L^{1}$ error

Figure 1: The averaged error for each method is normalized by the averaged
error incurred by HesScale. We show 40 initialization points with the same
colors across all methods. The norm of the vector of Hessian diagonals
$|\text{diag}({\bm{H}})|$ is shown as a reference.

Different Hessian diagonal approximations are considered for comparison with
HesScale. We included several deterministic and stochastic approximations for
the Hessian diagonals. We also include the approximation of the Fisher
Information Matrix done by squaring the gradients and denoted by $g^{2}$,
which is highly adopted by many first-order methods (e.g., Kingma and Ba,
2015). We compare HesScale with three stochastic approximation methods:
AdaHessian (Yao et al. 2021), Kronecker-factored approximate curvature (KFAC)
(Martens & Grosse 2015), and the Monte-Carlo (MC) estimate of the GGN matrix
(GGN-MC) (Dangel et al. 2020). We also compare HesScale with two deterministic
approximation methods: the diagonals of the exact GGN matrix (Schraudolph
2002) ($\text{diag}({\bm{G}})$) and the diagonal approximation by Becker and
Lecun (1989) (BL89). In KFAC, we extract the diagonals from the block diagonal
matrix and show the approximation error averaged over 1 MC sample (KFAC-MC1)
and over 50 MC samples (KFAC-MC50), both per each data example. Since
AdaHessian and GGN-MC are already diagonal approximations, we use them
directly and show the error with 1 MC sample (GGN-MC1 & AdaHessian-MC1) and
with 50 MC samples (GGN-MC50 & AdaHessian-MC50).

HesScale provides a better approximation than the other deterministic and
stochastic methods. For stochastic methods, we use many MC samples to improve
their approximation. However, their approximation quality is still poor.
Methods approximating the GGN diagonals do not capture the complete Hessian
information since the GGN and Hessian matrices are different when the
activation functions are not piece-wise linear. Although these methods
approximate the GGN diagonals, their approximation is significantly better
than the AdaHessian approximation. And among the methods for approximating the
GGN diagonals, HesScaleGN performs the best and close to the exact GGN
diagonals. This experiment clearly shows that HesScale achieves the best
approximation quality compared to other stochastic and deterministic
approximation methods.

Next, we perform another experiment to evaluate the computational cost of our
optimizers. Our Hessian approximation methods and corresponding optimizers
have linear computational complexity, which can be seen from Eq. 4 and Eq. 5.
However, computing second-order information in optimizers still incurs extra
computations compared to first-order optimizers, which may impact how the
total computations scale with the number of parameters. Hence, we compare the
computational cost of our optimizers with others for various number of
parameters. More specifically, we measure the update time of each optimizer,
which is the time needed to backpropagate first-order and second-order
information and update the parameters.

We designed two experiments to study the computational cost of first-order and
second-order optimizers. In the first experiment, we used a neural network
with a single hidden layer. The network has 64 inputs and 512 hidden units
with tanh activations. We study the increase in computational time when
increasing the number of outputs exponentially, which roughly doubles the
number of parameters. The set of values we use for the number of outputs is
$\\{2^{4},2^{5},2^{6},2^{7},2^{8},2^{9}\\}$. The results of this experiment
are shown in Fig. 2(a). In the second experiment, we used a neural network
with multi-layers, each containing 512 hidden units with tanh activations. The
network has 64 inputs and 100 outputs. We study the increase in computational
time when increasing the number of layers exponentially, which also roughly
doubles the number of parameters. The set of values we use for the number of
layers is $\\{1,2,4,8,16,32,64,128\\}$. The results of this experiment are
shown in Fig. 2(b). The points in Fig. 2(a) and Fig. 2(b) are averaged over 30
updates. The standard errors of the means of these points are smaller than the
width of each line. On average, we notice that the cost of AdaHessian,
AdaHesScale, and AdaHesScaleGN are three, two, and 1.25 times the cost of
Adam, respectively.

(a) Increasing number of outputs in a neural network

(b) Increasing number of layers in a neural network

Figure 2: The average computation time for each step of an update is shown for
different optimizers. The computed update time is the time needed by each
optimizer to backpropagate gradients or second-order information and to update
the parameters of the network. GGN overlaps with H in (a).

It is clear that our methods are among the most computationally efficient
approximation method for Hessian diagonals.

##  5 Empirical performance of HesScale in Optimization

In this section, we compare the performance of our optimizers—AdaHesScale and
AdaHesScaleGN—with three second-order optimizers: BL89, GGNMC, and AdaHessian.
We also include comparisons to two first-order methods: Adam and SGD. We
exclude KFAC and the exact diagonals of the GGN matrix from our comparisons
due to their prohibitive computations.

Our optimizers are evaluated in the supervised classification problem with a
series of experiments using different architectures and three datasets: MNIST,
CIFAR-10, and CIFAR-100. Instead of attempting to achieve state-of-the-art
performance with specialized techniques and architectures, we follow the
DeepOBS benchmarking work (Schneider et al. 2019) and compare the optimizers
in their generic and pristine form using relatively simple networks. It allows
us to perform a more fair comparison without extensively utilizing specialized
knowledge for a particular task. In the first experiment, we use the MNIST-MLP
task from DeepOBS. The images are flattened and used as inputs to a network of
three fully connected layers (1000, 500, and 100 units) with tanh activations.
We train each method for 100 epochs with a batch size of 128. We show the
training plots in Fig. 6(a) with their corresponding sensitivity plots in
Appendix C, Fig. 8(a). In the second experiment, we use the CIFAR10-3C3D task
from the DeepOBS benchmarking tasks. The network consists of three
convolutional layers with _tanh_ activations, each followed by max pooling.
After that, two fully connected layers (512 and 256 units) with _tanh_
activations are used. We train each method for 100 epochs with a batch size of
128. We show the training plots in Fig. 6(b) with their corresponding
sensitivity plots in Fig. 8(b). In the third experiment, we use the
CIFAR100-3C-3D task from DeepOBS. The network is the same as the one used in
the second task except for the activations are _ELU_. We train each method for
200 epochs with a batch size of 128. We show the training plots in Fig. 7(b)
with their corresponding sensitivity plots in Fig. 9(b). In the fourth
experiment, we use the CIFAR100-ALL-CNN task from DeepOBS with the ALL-CNN-C
network, which consists of 9 convolutional layers (Springenberg et al. 2015).
We use ELU activations insead of ReLU to differentiate between the performance
of AdaHesScale and AdaHesScaleGN. We show the training plots in Fig. 7(a) with
their corresponding sensitivity plots in Fig. 9(a).

In the MNIST-MLP and CIFAR-10-3C3D experiments, we performed a hyperparameter
search for each method to determine the best set of $\beta_{1}$, $\beta_{2}$,
and $\alpha$. The range of $\beta_{2}$ is $\\{0.99,0.999,0.9999\\}$ and the
range of $\beta_{1}$ is $\\{0.0,0.9\\}$. The range of step size is selected
for each method to create a convex curve. Our criterion was to find the best
hyperparameter configuration for each method in the search space that
minimizes the area under the validation loss curve. The performance of each
method was averaged over 30 independent runs. Each independent run had the
same initial representation for the algorithms used in an experiment. Using
each method’s best hyperparameter configuration on the validation set, we show
the performance of each method against the time in seconds needed to complete
the required number of epochs, which better depicts the computational
efficiency of the methods. Fig. 3(a) and Fig. 3(b) show these results on
MNIST-MLP and CIFAR-10 tasks. Moreover, we show the sensitivity of each method
to the step size in Fig. 5(a) and Fig. 5(b).

(a) MNIST-MLP

(b) CIFAR-10 3C3D

Figure 3: MNIST-MLP and CIFAR-10 3C3D classification tasks. Each method is
trained for 100 epochs. We show the time taken by each algorithm in seconds
(left) and we show the learning curves in the number of epochs (right). The
performance of each method is averaged over 30 independent runs. The shaded
area represents the standard error.

(a) CIFAR-100 3C3D

(b) CIFAR-100 All-CNN-C

Figure 4: CIFAR-100 3C3D and CIFAR-100 ALL-CNN classification tasks. Each
method from the first task is trained for 200 epochs and each method from the
second task is trained for 350 epochs. We show the time taken by each
algorithm in seconds (left) and we show the learning curves in the number of
epochs (right). The performance of each method is averaged over 30 independent
runs. The shaded area represents the standard error.

(a) MNIST

(b) CIFAR-10

Figure 5: Sensitivity of the step size for each method on MNIST-MLP and
CIFAR-10 3C3D tasks. We select the best values of $\beta_{1}$ and $\beta_{2}$
for each step size $\alpha$.

In the CIFAR-100-ALL-CNN and CIFAR-100-3C3D experiments, we used the set of
$\beta_{1}$ and $\beta_{2}$ that achieved the best robustness in the previous
two tasks. The values used for $\beta_{1}$ and $\beta_{2}$ were $0.9$ and
$0.999$ respectively. We did a hyperparameter search for each method to
determine the best step size using the specified $\beta_{1}$ and $\beta_{2}$.
The range of step size is selected for each method to create a convex curve,
while ensuring that each method uses a similar search budget and the same
interval between search points. Our criterion was to find the best
hyperparameter configuration for each method in the search space that
minimizes the area under the validation loss curve. The performance of each
method was averaged over 30 independent runs. Each independent run had the
same initial representation for the algorithms used in an experiment. Using
each method’s best hyperparameter configuration on the validation set, we show
the performance of each method against the time in seconds needed to complete
the required number of epochs. Fig. 4(a) and Fig. 4(b) show these results on
CIFAR-100-ALL-CNN and CIFAR-100-3C3D tasks.

Our results show that all optimizers except for BL89 performed well on the
MNIST-MLP task. However, in CIFAR-10, CIFAR-100 3c3d, and CIFAR-100 ALL-CNN,
we notice that AdaHessian performed worse than all methods except BL89. This
result is aligned with AdaHessian’s inability to accurately approximate the
Hessian diagonals, as shown in Fig. 1. Moreover, AdaHessian required more
computational time compared to all methods, which is also reflected in Fig. 2.
While being time-efficient, AdaHesScaleGN consistently outperformed all
methods in CIFAR-10-3C3D and CIFAR-100-3C3D, and it outperformed all methods
except AdaHesScale in CIFAR-100 ALL-CNN. This result is aligned with our
methods’ accurate approximation of Hessian diagonals. Our experiments indicate
that incorporating HesScale and HesScaleGN approximations in optimization
methods can be of significant performance advantage in both computation and
accuracy. AdaHesScale and AdaHesScaleGN outperformed other optimizers likely
due to their accurate approximation of the diagonals of the Hessian and GGN,
respectively.

##  6 Conclusion

HesScale is a scalable and efficient second-order method for approximating the
diagonals of the Hessian at every network layer. Our work is based on the
previous work done by Becker and Lecun (1989). We performed a series of
experiments to evaluate HesScale against other scalable algorithms in terms of
computational cost and approximation accuracy. Moreover, we demonstrated how
HesScale can be used to build efficient second-order optimization methods. Our
results showed that our methods provide a more accurate approximation and
require small additional computations.

##  7 Broader Impact

Second-order information is used in domains other than optimization. For
example, some works alleviating catastrophic forgetting use a utility measure
for the network’s connections to protect them. Typically, an auxiliary loss is
used between such connections, and their old values are weighted by their
corresponding importance. Such methods (LeCun et al. 1990, Hassibi & Stork
1993, Dong et al. 2017, Kirkpatrick et al. 2017, Schwarz et al. 2018, Ritter
et al. 2018) use the diagonal of the Fisher information matrix or the Hessian
matrix as a utility measure for each weight. The quality of these algorithms
depends heavily on the approximation quality of the second-order
approximation. Second-order information can also be used in neural network
pruning. Molchanov et al. (2019) showed that second-order approximation with
the exact Hessian diagonals could closely represent the true measure of the
utility of each weight.

The accurate and efficient approximation for the diagonals of the Hessian at
each layer enables HesScale to be used in many important lines of research.
Using this second-order information provides a reliable measure of connection
utility. Therefore, using HesScale in these types of problems can potentially
improve the performance of neural network pruning methods and regularization-
based catastrophic forgetting prevention methods.

#### Acknowledgments

We gratefully acknowledge funding from the Canada CIFAR AI Chairs program, the
Reinforcement Learning and Artificial Intelligence (RLAI) laboratory, the
Alberta Machine Intelligence Institute (Amii), and the Natural Sciences and
Engineering Research Council (NSERC) of Canada. We would also like to thank
Compute Canada for providing the computational resources needed.

## References

Amari, S. (1998). Natural Gradient Works Efficiently in Learning. Neural
Computation, 10(2), 251–276.

Ba, J. L., Kiros, J. R., & Hinton, G. E. (2016). Layer normalization. arXiv
preprint arXiv:1607.06450.

Becker, S. & Lecun, Y. (1989). Improving the convergence of back-propagation
learning with second-order methods. Proceedings of the 1988 Connectionist
Models Summer School (pp. 29–37).

Bekas, C., Kokiopoulou, E., & Saad, Y. (2007). An estimator for the diagonal
of a matrix. Applied Numerical Mathematics, 57(11), 1214–1229.

Botev, A., Ritter, H., & Barber, D. (2017). Practical Gauss-Newton
optimisation for deep learning. Proceedings of the 34th International
Conference on Machine Learning, 70, 557–565.

Chan, A., Silva, H., Lim, S., Kozuno, T., Mahmood, A. R., & White, M. (2022).
Greedification operators for policy optimization: Investigating forward and
reverse KL divergences. Journal of Machine Learning Research, 23(253), 1-79.

Chapelle, O. & Erhan, D. (2011). Improved preconditioner for hessian free
optimization. NIPS Workshop on Deep Learning and Unsupervised Feature
Learning.

Dangel, F., Kunstner, F., & Hennig, P. (2020). BackPACK: Packing more into
backprop. International Conference on Learning Representations.

Dong, X., Chen, S., & Pan, S. J. (2017). Learning to prune deep neural
networks via layer-wise optimal brain surgeon. Proceedings of the 31st
International Conference on Neural Information Processing Systems (pp.
4860-4874).

Duchi, J., Hazan, E., & Singer, Y. (2011). Adaptive subgradient methods for
online learning and stochastic optimization. Journal of Machine Learning
Research, 12(61), 2121–2159.

Hassibi, B. & Stork, D. (1993). Second order derivatives for network pruning:
Optimal brain surgeon. Advances in Neural Information Processing Systems, 5.

Ioffe, S. and Szegedy, C. (2015). Batch normalization: Accelerating deep
network training by reducing internal covariate shift. International
conference on machine learning (pp. 448-456).

Kingma, D. P. & Ba, J. (2015). Adam: A method for stochastic optimization.
International Conference on Learning Representations.

Kirkpatrick, J., Pascanu, R., Rabinowitz, N., Veness, J., Desjardins, G.,
Rusu, A. A., Milan, K., Quan, J., Ramalho, T., Grabska-Barwinska, A., et al.
(2017). Overcoming catastrophic forgetting in neural networks. Proceedings of
the national academy of sciences, 114(13), 3521–3526.

Kunstner, F., Hennig, P., & Balles, L. (2019). Limitations of the empirical
fisher approximation for natural gradient descent. Advances in Neural
Information Processing Systems, 32.

LeCun, Y., Denker, J., & Solla, S. (1990). Optimal brain damage. Advances in
Neural Information Processing Systems, 2.

Luo, L., Xiong, Y., & Liu, Y. (2019). Adaptive gradient methods with dynamic
bound of learning rate. International Conference on Learning Representations.

Martens, J. (2010). Deep learning via hessian-free optimization. International
Conference on Machine Learning (pp. 735–742).

Martens, J. (2020). New insights & perspectives on the natural gradient
method. Journal of Machine Learning Research, 21(146), 1-76.

Martens, J. & Grosse, R. (2015). Optimizing neural networks with kronecker-
factored approximate curvature. International Conference on Machine Learning
(pp. 2408–2417).

Martens, J. & Sutskever, I. (2011). Learning recurrent neural networks with
hessian-free optimization. International Conference on Machine Learning (pp.
1033-1040).

Martens, J., Sutskever, I., & Swersky, K. (2012). Estimating the hessian by
back-propagating curvature. International Conference on International
Conference on Machine Learning(pp. 963–970).

Mizutani, E. & Dreyfus, S. E. (2008). Second-order stagewise backpropagation
for hessian-matrix analyses & investigation of negative curvature. Neural
Networks, 21(2), 193–203.

Molchanov, P., Mallya, A., Tyree, S., Frosio, I., & Kautz, J. (2019).
Importance estimation for neural network pruning. IEEE/CVF Conference on
Computer Vision and Pattern Recognition (pp. 11264–11272).

Reddi, S. J., Kale, S., & Kumar, S. (2018). On the convergence of Adam and
beyond. International Conference on Learning Representations.

Ritter, H., Botev, A., and Barber, D. (2018). Online structured Laplace
approximations for overcoming catastrophic forgetting. International
Conference on Neural Information Processing Systems (pp. 3742–3752).

Rumelhart, D. E., Hinton, G. E., and Williams, R. J. (1986). Learning
representations by back-propagating errors. Nature, 323(6088), 533–536.

Salimans, T. and Kingma, D. P. (2016). Weight normalization: A simple
reparameterization to accelerate training of deep neural networks. Advances in
neural information processing systems, 29, 901–909.

Schneider, F., Balles, L., and Hennig, P. (2019). DeepOBS: A deep learning
optimizer benchmark suite. International Conference on Learning
Representations.

Schraudolph, N. N. (2002). Fast curvature matrix-vector products for second-
order gradient descent. Neural Computation, 14(7), 1723–1738.

Schwarz, J., Czarnecki, W., Luketina, J., Grabska-Barwinska, A., Teh, Y. W.,
Pascanu, R., and Hadsell, R. (2018). Progress & compress: A scalable framework
for continual learning. International Conference on Machine Learning (pp.
4528–4537).

Springenberg, J., Dosovitskiy, A., Brox, T., and Riedmiller, M. (2015).
Striving for simplicity: The all convolutional net. International Conference
on Learning Representations [Workshop].

Tieleman, T., Hinton, G., et al. (2012). Lecture 6.5-RMSProp: Divide the
gradient by a running average of its recent magnitude. COURSERA: Neural
networks for machine learning, 4(2), 26–31.

Tran, P. T. and Phong, L. T. (2019). On the convergence proof of AMSGrad and a
new version. IEEE Access, 7, 61706–61716.

Wilson, A. C., Roelofs, R., Stern, M., Srebro, N., and Recht, B. (2017). The
marginal value of adaptive gradient methods in machine learning. International
Conference on Neural Information Processing Systems (pp. 4151–4161).

Yao, Z., Gholami, A., Shen, S., Keutzer, K., and Mahoney, M. W. (2021).
Adahessian: An adaptive second order optimizer for machine learning. AAAI
Conference on Artificial Intelligence, 35(12), 10665-10673.

Zeiler, M. D. (2012). Adadelta: An adaptive learning rate method. arXiv
preprint arXiv:1212.5701.

##  Appendix A Hessian diagonals of the log-likelihood function for two common
distributions

Here, we provide the diagonals of the Hessian matrix of functions involving
the log-likelihood of two common distributions: a normal distribution and a
categorical distribution with probabilities represented by a softmax function,
which we refer to as a _softmax distribution_. We show that the exact
computations of the diagonal can be computed with linear complexity since
computing the diagonal elements does not depend on off-diagonals in these
cases. In the following, we consider the softmax and normal distributions, and
we write the exact Hessian diagonals in both cases.

###  A.1 Softmax distribution

Consider a cross-entropy function for a discrete probability distribution as
$f\doteq-\sum_{i=1}^{|{\bm{q}}|}p_{i}\log q_{i}(\bm{\theta})$, where
${\bm{q}}$ is a probability vector that depends on a parameter vector
$\bm{\theta}$, and ${\bm{p}}$ is a one-hot vector for the target class. For
softmax distributions, ${\bm{q}}$ is parametrized by a softmax function
${\bm{q}}\doteq{e^{\bm{\theta}}/\sum_{i=1}^{|{\bm{q}}|}e^{\theta_{i}}}$. In
this case, we can write the gradient of the cross-entropy function with
respect to $\bm{\theta}$ as

$\nabla_{\bm{\theta}}f(\bm{\theta})={\bm{q}}-{\bm{p}}.$

Next, we write the exact diagonal elements of the Hessian matrix as follows:

$\text{diag}({\bm{H}}_{\bm{\theta}})=\text{diag}(\nabla_{\bm{\theta}}({\bm{q}}-{\bm{p}}))={\bm{q}}-{\bm{q}}^{2},$

where ${\bm{q}}^{2}$ denotes element-wise squaring of ${\bm{q}}$, and $\nabla$
operator applied to a vector denotes Jacobian. Computing the exact diagonals
of the Hessian matrix depends only on vector operations, which means that we
can compute it in $O(n)$. The cross-entropy loss is used with softmax
distribution in many important tasks, such as supervised classification and
discrete reinforcement learning control with parameterized policies (Chan et
al. 2022).

###  A.2 Multivariate normal distribution with diagonal covariance

For a multivariate normal distribution with diagonal covariance, the parameter
vector $\bm{\theta}$ is determined by the mean-variance vector pair:
$\bm{\theta}\doteq(\bm{\mu},{\bm{\sigma}}^{2})$. The log-likelihood of a
random vector ${\bm{x}}$ drawn from this distribution can be written as

$\displaystyle\log q({\bm{x}};\bm{\mu},{\bm{\sigma}}^{2})$
$\displaystyle=-\frac{1}{2}({\bm{x}}-\bm{\mu})^{\top}\bm{D}({\bm{\sigma}}^{2})^{-1}({\bm{x}}-\bm{\mu})-\frac{1}{2}\log(|\bm{D}({\bm{\sigma}}^{2})|)+c$
$\displaystyle=-\frac{1}{2}({\bm{x}}-\bm{\mu})^{\top}\bm{D}({\bm{\sigma}}^{2})^{-1}({\bm{x}}-\bm{\mu})-\frac{1}{2}\log(\sum_{i=1}^{|{\bm{\sigma}}|}\sigma^{2}_{i})+c,$

where $\bm{D}({\bm{\sigma}}^{2})$ gives a diagonal matrix with
${\bm{\sigma}}^{2}$ in its diagonal, $|{\bm{M}}|$ is the determinant of a
matrix ${\bm{M}}$ and $c$ is some constant. We can write the gradients of the
log-likelihood function with respect to $\bm{\mu}$ and ${\bm{\sigma}}^{2}$ as
follows:

$\displaystyle\nabla_{\bm{\mu}}\log q({\bm{x}};\bm{\mu},{\bm{\sigma}}^{2})$
$\displaystyle=\bm{D}({\bm{\sigma}}^{2})^{-1}({\bm{x}}-\bm{\mu})=({\bm{x}}-\bm{\mu})\oslash{\bm{\sigma}}^{2},$
$\displaystyle\nabla_{{\bm{\sigma}}^{2}}\log
q({\bm{x}};\bm{\mu},{\bm{\sigma}}^{2})$
$\displaystyle=\frac{1}{2}\big{[}({\bm{x}}-\bm{\mu})^{2}\oslash{\bm{\sigma}}^{2}-\mathbf{1}\big{]}\oslash{\bm{\sigma}}^{2},$

where $\mathbf{1}$ is an all-ones vector, and $\oslash$ denotes element-wise
division. Finally, we write the exact diagonals of the Hessian matrix as

$\text{diag}({\bm{H}}_{\bm{\mu}})=\text{diag}(\nabla_{{\bm{\mu}}}({\bm{x}}-\bm{\mu})\oslash{\bm{\sigma}}^{2})=-\mathbf{1}\oslash{\bm{\sigma}}^{2},$
$\text{diag}({\bm{H}}_{{\bm{\sigma}}^{2}})=\text{diag}\Big{(}\nabla_{{\bm{\sigma}}^{2}}\big{[}\frac{1}{2}[({\bm{x}}-\bm{\mu})^{2}\oslash{\bm{\sigma}}^{2}-\mathbf{1}]\oslash{\bm{\sigma}}^{2}\big{]}\Big{)}=\big{[}0.5\mathbf{1}-({\bm{x}}-\bm{\mu})^{2}\oslash{\bm{\sigma}}^{2}\big{]}\oslash{\bm{\sigma}}^{4}.$

Clearly, the gradient and the exact Hessian diagonals can be computed in
$O(n)$. Log-likelihood functions for normal distributions are used in many
important problems, such as variational inference and continuous reinforcement
learning control.

##  Appendix B HesScale with Convolutional Neural Networks

Here, we derive the Hessian propagation for convolutional neural networks
(CNNs). Consider a CNN with $L-1$ layers followed by a fully connected layer
that outputs the predicted output ${\bm{q}}$. The CNN filters are
parameterized by $\\{{\bm{W}}_{1},...,{\bm{W}}_{L}\\}$, where ${\bm{W}}_{l}$
is the filter matrix at the $l$-th layer with the dimensions $k_{l,1}\times
k_{l,2}$, and its element at the $i$th row and the $j$th column is denoted by
$W_{l,i,j}$. For the simplicity of this proof, we assume that the number of
filters at each layer is one; the proof can be extended easily to the general
case. The learning algorithm learns the target function $f^{*}$ by optimizing
the loss $\mathcal{L}$. During learning, the parameters of the neural network
are changed to reduce the loss. At the layer $l$, we get the activation output
matrix ${\bm{H}}_{l}$ by applying the activation function $\bm{\sigma}$ to the
activation input ${\bm{A}}_{l}$: ${\bm{H}}_{l}=\bm{\sigma}({\bm{A}}_{l})$. We
assume here that the activation function is element-wise activation for all
layers except for the final layer $L$, where it becomes the softmax function.
We simplify notations by defining ${\bm{H}}_{0}\doteq{\bm{X}}$, where
${\bm{X}}$ is the input sample. The activation output ${\bm{H}}_{l}$ is then
convoluted by the weight matrix ${\bm{W}}_{l+1}$ of layer $l+1$ to produce the
next activation input:
${A}_{l+1,i,j}=\sum_{m=0}^{k_{l,1}-1}\sum_{n=0}^{k_{l,2}-1}{{W}_{l+1,m,n}}{H}_{l,(i+m),(j+n)}$.
We denote the size of the activation output at the $l$-th layer by
$h_{l}\times w_{l}$. The backpropagation equations for the described network
are given following Rumelhart et al. (1986):

$\displaystyle\frac{\partial\mathcal{L}}{\partial A_{l,i,j}}$
$\displaystyle=\sum_{m=0}^{k_{l+1,1}-1}\sum_{n=0}^{k_{l+1,2}-1}\frac{\partial\mathcal{L}}{\partial
A_{l+1,(i-m),(j-n)}}\frac{\partial A_{l+1,(i-m),(j-n)}}{\partial A_{l,i,j}}$
$\displaystyle=\sum_{m=0}^{k_{l+1,1}-1}\sum_{n=0}^{k_{l+1,2}-1}\frac{\partial\mathcal{L}}{\partial
A_{l+1,(i-m),(j-n)}}\sum_{m^{\prime}=0}^{k_{l+1,1}-1}\sum_{n^{\prime}=0}^{k_{l+1,2}-1}W_{l+1,m^{\prime},n^{\prime}}\frac{\partial
H_{l,(i-m+m^{\prime}),(j-n+n^{\prime})}}{\partial A_{l,i,j}}$
$\displaystyle=\sum_{m=0}^{k_{l+1,1}-1}\sum_{n=0}^{k_{l+1,2}-1}\frac{\partial\mathcal{L}}{\partial
A_{l+1,(i-m),(j-n)}}W_{l+1,m,n}\sigma^{\prime}(A_{l,i,j})$
$\displaystyle=\sigma^{\prime}(A_{l,i,j})\sum_{m=0}^{k_{l+1,1}-1}\sum_{n=0}^{k_{l+1,2}-1}\frac{\partial\mathcal{L}}{\partial
A_{l+1,(i-m),(j-n)}}W_{l+1,m,n},$ (6)
$\displaystyle\frac{\partial\mathcal{L}}{\partial W_{l,i,j}}$
$\displaystyle=\sum_{m=0}^{h_{l}-k_{l,1}}\sum_{n=0}^{w_{l}-k_{l,2}}\frac{\partial\mathcal{L}}{\partial
A_{l,m,n}}\frac{\partial A_{l,m,n}}{\partial
W_{l,i,j}}=\sum_{m=0}^{h_{l}-k_{l,1}}\sum_{n=0}^{w_{l}-k_{l,2}}\frac{\partial\mathcal{L}}{\partial
A_{l,m,n}}H_{l-1,(i+m),(j+n)}.$ (7)

In the following, we write the equations for the exact Hessian diagonals with
respect to weights
$\nicefrac{{\partial^{2}\mathcal{L}}}{{\partial{W^{2}_{l,i,j}}}}$, which
requires the calculation of
$\nicefrac{{\partial^{2}\mathcal{L}}}{{\partial{A^{2}_{l,i,j}}}}$ first:

$\displaystyle\frac{\partial^{2}\mathcal{L}}{\partial{A^{2}_{l,i,j}}}$
$\displaystyle=\frac{\partial}{\partial
A_{l,i,j}}\Bigg{[}\sigma^{\prime}(A_{l,i,j})\sum_{m=0}^{k_{l+1,1}-1}\sum_{n=0}^{k_{l+1,2}-1}\frac{\partial\mathcal{L}}{\partial
A_{l+1,(i-m),(j-n)}}W_{l+1,m,n}\Bigg{]}$
$\displaystyle=\sigma^{\prime}(A_{l,i,j})\sum_{m,p=0}^{k_{l+1,2}-1}\sum_{n,q=0}^{k_{l+1,2}-1}\frac{\partial^{2}\mathcal{L}}{\partial
A_{l+1,(i-m),(j-n)}\partial A_{l+1,(i-p),(j-q)}}\frac{\partial
A_{l+1,(i-p),(j-q)}}{\partial A_{l,i,j}}W_{l+1,m,n}$
$\displaystyle+\sigma^{\prime\prime}(A_{l,i,j})\sum_{m=0}^{k_{l+1,2}-1}\sum_{n=0}^{k_{l+1,2}-1}\frac{\partial\mathcal{L}}{\partial
A_{l+1,(i-m),(j-n)}}W_{l+1,m,n}$
$\displaystyle\frac{\partial^{2}\mathcal{L}}{\partial{W^{2}_{l,i,j}}}$
$\displaystyle=\frac{\partial}{\partial
W_{l,i,j}}\Bigg{[}\sum_{m=0}^{h_{l}-k_{l,1}}\sum_{n=0}^{w_{l}-k_{l,2}}\frac{\partial\mathcal{L}}{\partial
A_{l,m,n}}H_{l-1,(i+m),(j+n)}\Bigg{]}$
$\displaystyle=\sum_{m,p=0}^{h_{l}-k_{l,1}}\sum_{n,q=0}^{w_{l}-k_{l,2}}\frac{\partial^{2}\mathcal{L}}{\partial
A_{l,m,n}\partial A_{l,p,q}}\frac{\partial A_{l,p,q}}{\partial
W_{l,i,j}}H_{l-1,(i+m),(j+n)}$

Since the calculation of
$\nicefrac{{\partial^{2}\mathcal{L}}}{{\partial{A^{2}_{l,i,j}}}}$ and
$\nicefrac{{\partial^{2}\mathcal{L}}}{{\partial{W^{2}_{l,i,j}}}}$ depend on
the off-diagonal terms, the computation complexity becomes quadratic.
Following Becker and Lecun (1989), we approximate the Hessian diagonals by
ignoring the off-diagonal terms, which leads to a backpropagation rule with
linear computational complexity for our estimates
$\widehat{\frac{\partial^{2}\mathcal{L}}{\partial{W^{2}_{l,i,j}}}}$ and
$\widehat{\frac{\partial^{2}\mathcal{L}}{\partial{A^{2}_{l,i,j}}}}$:

$\displaystyle\widehat{\frac{\partial^{2}\mathcal{L}}{\partial{A^{2}_{l,i,j}}}}$
$\displaystyle\doteq\sigma^{\prime}(A_{l,i,j})^{2}\sum_{m=0}^{k_{l+1,2}-1}\sum_{n=0}^{k_{l+1,2}-1}\widehat{\frac{\partial^{2}\mathcal{L}}{\partial
A^{2}_{l+1,(i-m),(j-n)}}}W^{2}_{l+1,m,n}$
$\displaystyle+\sigma^{\prime\prime}(A_{l,i,j})\sum_{m=0}^{k_{l+1,2}-1}\sum_{n=0}^{k_{l+1,2}-1}\frac{\partial\mathcal{L}}{\partial
A_{l+1,(i-m),(j-n)}}W_{l+1,m,n},$ (8)
$\displaystyle\widehat{\frac{\partial^{2}\mathcal{L}}{\partial{W^{2}_{l,i,j}}}}$
$\displaystyle\doteq\sum_{m=0}^{h_{l}-k_{l,1}}\sum_{n=0}^{w_{l}-k_{l,2}}\widehat{\frac{\partial^{2}\mathcal{L}}{\partial
A^{2}_{l,m,n}}}H^{2}_{l-1,(i+m),(j+n)}.$ (9)

##  Appendix C Optimization plots in the number of epochs

We give the training loss, training accuracy, validation loss, validation
accuracy, test loss, and test accuracy for each of the methods we include in
our comparison in Fig. 6 and Fig. 7. Moreover, we give the sensitivity plots
for $\beta_{1}$, $\beta_{2}$, and $\alpha$ for each method in Fig. 8 and Fig.
9.

(a) MNIST

(b) CIFAR-10

Figure 6: Learning curves of each algorithm on two tasks, MNIST-MLP and
CIFAR-10 3C3D, for 100 epochs. We show the best configuration for each
algorithm on the validation set. The best parameter configuration for each
algorithm is selected based on the area under the curve for the validation
loss.

(a) CIFAR-100 All CNN

(b) CIFAR-100 3C3D

Figure 7: Learning Curves of each algorithm on CIFAR-100 with All-CNN and 3C3D
architectures, for 100 epochs. We show the best configuration for each
algorithm on the validation set. The best parameter configuration for each
algorithm is selected based on the area under the curve for the validation
loss.

(a) MNIST

(b) CIFAR-10

Figure 8: Parameter Sensitivity study for each algorithm on two data sets,
MNIST and CIFAR-10. The range of $\beta_{2}$ is $\\{0.99,0.999,0.9999\\}$ and
the range of $\beta_{1}$ is $\\{0.0,0.9\\}$. Each point for each algorithm
represents the average test loss given a set of parameters.

(a) CIFAR-100 All-CNN

(b) CIFAR-100 3C3D

Figure 9: Parameter Sensitivity study for each algorithm on CIFAR-100 with
All-CNN and 3C3D architectures. The range of step size is
$\\{10^{-5},10^{-4},10^{-3},10^{-2},10^{-1},10^{0}\\}$. We choose $\beta_{1}$
to be equal to $0.9$ and $\beta_{2}$ to be equal to $0.999$. Each point for
each algorithm represents the average test loss given a set of parameters.


# Intriguing properties of neural networks

Christian Szegedy
Google Inc.
&Wojciech Zaremba
New York University
&Ilya Sutskever
Google Inc.
&Joan Bruna
New York University
&Dumitru Erhan
Google Inc.
&Ian Goodfellow
University of Montreal
&Rob Fergus
New York University
Facebook Inc.

###### Abstract

Deep neural networks are highly expressive models that have recently achieved
state of the art performance on speech and visual recognition tasks. While
their expressiveness is the reason they succeed, it also causes them to learn
uninterpretable solutions that could have counter-intuitive properties. In
this paper we report two such properties.

First, we find that there is no distinction between individual high level
units and random linear combinations of high level units, according to various
methods of unit analysis. It suggests that it is the space, rather than the
individual units, that contains the semantic information in the high layers of
neural networks.

Second, we find that deep neural networks learn input-output mappings that are
fairly discontinuous to a significant extent. We can cause the network to
misclassify an image by applying a certain hardly perceptible perturbation,
which is found by maximizing the network’s prediction error. In addition, the
specific nature of these perturbations is not a random artifact of learning:
the same perturbation can cause a different network, that was trained on a
different subset of the dataset, to misclassify the same input.

##  1 Introduction

Deep neural networks are powerful learning models that achieve excellent
performance on visual and speech recognition problems [9, 8]. Neural networks
achieve high performance because they can express arbitrary computation that
consists of a modest number of massively parallel nonlinear steps. But as the
resulting computation is automatically discovered by backpropagation via
supervised learning, it can be difficult to interpret and can have counter-
intuitive properties. In this paper, we discuss two counter-intuitive
properties of deep neural networks.

The first property is concerned with the semantic meaning of individual units.
Previous works [6, 13, 7] analyzed the semantic meaning of various units by
finding the set of inputs that maximally activate a given unit. The inspection
of individual units makes the implicit assumption that the units of the last
feature layer form a distinguished basis which is particularly useful for
extracting semantic information. Instead, we show in section 3 that random
projections of $\phi(x)$ are semantically indistinguishable from the
coordinates of $\phi(x)$. This puts into question the conjecture that neural
networks disentangle variation factors across coordinates. Generally, it seems
that it is the entire space of activations, rather than the individual units,
that contains the bulk of the semantic information. A similar, but even
stronger conclusion was reached recently by Mikolov et al. [12] for word
representations, where the various directions in the vector space representing
the words are shown to give rise to a surprisingly rich semantic encoding of
relations and analogies. At the same time, the vector representations are
stable up to a rotation of the space, so the individual units of the vector
representations are unlikely to contain semantic information.

The second property is concerned with the stability of neural networks with
respect to small perturbations to their inputs. Consider a state-of-the-art
deep neural network that generalizes well on an object recognition task. We
expect such network to be robust to small perturbations of its input, because
small perturbation cannot change the object category of an image. However, we
find that applying an _imperceptible_ non-random perturbation to a test image,
it is possible to arbitrarily change the network’s prediction (see figure 5).
These perturbations are found by optimizing the input to maximize the
prediction error. We term the so perturbed examples “adversarial examples”.

It is natural to expect that the precise configuration of the minimal
necessary perturbations is a random artifact of the normal variability that
arises in different runs of backpropagation learning. Yet, we found that
adversarial examples are relatively robust, and are shared by neural networks
with varied number of layers, activations or trained on different subsets of
the training data. That is, if we use one neural net to generate a set of
adversarial examples, we find that these examples are still statistically hard
for another neural network even when it was trained with different
hyperparameters or, most surprisingly, when it was trained on a different set
of examples.

These results suggest that the deep neural networks that are learned by
backpropagation have nonintuitive characteristics and intrinsic blind spots,
whose structure is connected to the data distribution in a non-obvious way.

##  2 Framework

Notation We denote by $x\in\mathbb{R}^{m}$ an input image, and $\phi(x)$
activation values of some layer. We first examine properties of the image of
$\phi(x)$, and then we search for its blind spots.

We perform a number of experiments on a few different networks and three
datasets :

  * •

For the MNIST dataset, we used the following architectures [11]

    * –

A simple fully connected network with one or more hidden layers and a Softmax
classifier. We refer to this network as “FC”.

    * –

A classifier trained on top of an autoencoder. We refer to this network as
“AE”.

  * •

The ImageNet dataset [3].

    * –

Krizhevsky et. al architecture [9]. We refer to it as “AlexNet”.

  * •

$\sim 10$M image samples from Youtube (see [10])

    * –

Unsupervised trained network with $\sim$ 1 billion learnable parameters. We
refer to it as “QuocNet”.

For the MNIST experiments, we use regularization with a weight decay of
$\lambda$. Moreover, in some experiments we split the MNIST training dataset
into two disjoint datasets $P_{1}$, and $P_{2}$, each with 30000 training
cases.

##  3 Units of: $\phi(x)$

Traditional computer vision systems rely on feature extraction: often a single
feature is easily interpretable, e.g. a histogram of colors, or quantized
local derivatives. This allows one to inspect the individual coordinates of
the feature space, and link them back to meaningful variations in the input
domain. Similar reasoning was used in previous work that attempted to analyze
neural networks that were applied to computer vision problems. These works
interpret an activation of a hidden unit as a meaningful feature. They look
for input images which maximize the activation value of this single feature
[6, 13, 7, 4].

The aforementioned technique can be formally stated as visual inspection of
images $x^{\prime}$, which satisfy (or are close to maximum attainable value):

$\displaystyle
x^{\prime}=\operatorname*{arg\,max}_{x\in\mathcal{I}}\langle\phi(x),e_{i}\rangle$

where $\mathcal{I}$ is a held-out set of images from the data distribution
that the network was not trained on and $e_{i}$ is the natural basis vector
associated with the $i$-th hidden unit.

Our experiments show that any random direction $v\in\mathbb{R}^{n}$ gives rise
to similarly interpretable semantic properties. More formally, we find that
images $x^{\prime}$ are semantically related to each other, for many
$x^{\prime}$ such that

$\displaystyle
x^{\prime}=\operatorname*{arg\,max}_{x\in\mathcal{I}}\langle\phi(x),v\rangle$

This suggests that the natural basis is not better than a random basis for
inspecting the properties of $\phi(x)$. This puts into question the notion
that neural networks disentangle variation factors across coordinates.

First, we evaluated the above claim using a convolutional neural network
trained on MNIST. We used the MNIST test set for $\mathcal{I}$. Figure 1 shows
images that maximize the activations in the natural basis, and Figure 2 shows
images that maximize the activation in random directions. In both cases the
resulting images share many high-level similarities.

Next, we repeated our experiment on an AlexNet, where we used the validation
set as $\mathcal{I}$. Figures 3 and 4 compare the natural basis to the random
basis on the trained network. The rows appear to be semantically meaningful
for both the single unit and the combination of units.

(a) Unit sensitive to lower round stroke.

(b) Unit sensitive to upper round stroke, or lower straight stroke.

(c) Unit senstive to left, upper round stroke.

(d) Unit senstive to diagonal straight stroke.

Figure 1: An MNIST experiment. The figure shows images that maximize the
activation of various units (maximum stimulation in the natural basis
direction). Images within each row share semantic properties.

(a) Direction sensitive to upper straight stroke, or lower round stroke.

(b) Direction sensitive to lower left loop.

(c) Direction senstive to round top stroke.

(d) Direction sensitive to right, upper round stroke.

Figure 2: An MNIST experiment. The figure shows images that maximize the
activations in a random direction (maximum stimulation in a random basis).
Images within each row share semantic properties.

(a) Unit sensitive to white flowers.

(b) Unit sensitive to postures.

(c) Unit senstive to round, spiky flowers.

(d) Unit senstive to round green or yellow objects.

Figure 3: Experiment performed on ImageNet. Images stimulating single unit
most (maximum stimulation in natural basis direction). Images within each row
share many semantic properties.

(a) Direction sensitive to white, spread flowers.

(b) Direction sensitive to white dogs.

(c) Direction sensitive to spread shapes.

(d) Direction sensitive to dogs with brown heads.

Figure 4: Experiment performed on ImageNet. Images giving rise to maximum
activations in a random direction (maximum stimulation in a random basis).
Images within each row share many semantic properties.

Although such analysis gives insight on the capacity of $\phi$ to generate
invariance on a particular subset of the input distribution, it does not
explain the behavior on the rest of its domain. We shall see in the next
section that $\phi$ has counterintuitive properties in the neighbourhood of
almost every point form data distribution.

##  4 Blind Spots in Neural Networks

So far, unit-level inspection methods had relatively little utility beyond
confirming certain intuitions regarding the complexity of the representations
learned by a deep neural network [6, 13, 7, 4]. Global, network level
inspection methods _can_ be useful in the context of explaining classification
decisions made by a model [1] and can be used to, for instance, identify the
parts of the input which led to a correct classification of a given visual
input instance (in other words, one can use a trained model for weakly-
supervised localization). Such global analyses are useful in that they can
make us understand better the input-to-output mapping represented by the
trained network.

Generally speaking, the output layer unit of a neural network is a highly
nonlinear function of its input. When it is trained with the cross-entropy
loss (using the Softmax activation function), it represents a conditional
distribution of the label given the input (and the training set presented so
far). It has been argued [2] that the deep stack of non-linear layers in
between the input and the output unit of a neural network are a way for the
model to encode a _non-local generalization prior_ over the input space. In
other words, it is assumed that is possible for the output unit to assign non-
significant (and, presumably, non-epsilon) probabilities to regions of the
input space that contain no training examples in their vicinity. Such regions
can represent, for instance, the same objects from different viewpoints, which
are relatively far (in pixel space), but which share nonetheless both the
label and the statistical structure of the original inputs.

It is implicit in such arguments that $local$ generalization—in the very
proximity of the training examples—works as expected. And that in particular,
for a small enough radius $\varepsilon>0$ in the vicinity of a given training
input $x$, an $x+r$ satisfying $||r||<\varepsilon$ will get assigned a high
probability of the correct class by the model. This kind of smoothness prior
is typically valid for computer vision problems. In general, imperceptibly
tiny perturbations of a given image do not normally change the underlying
class.

Our main result is that for deep neural networks, the smoothness assumption
that underlies many kernel methods does not hold. Specifically, we show that
by using a simple optimization procedure, we are able to find adversarial
examples, which are obtained by imperceptibly small perturbations to a
correctly classified input image, so that it is no longer classified
correctly.

In some sense, what we describe is a way to traverse the manifold represented
by the network in an efficient way (by optimization) and finding _adversarial
examples_ in the input space. The adversarial examples represent low-
probability (high-dimensional) “pockets” in the manifold, which are hard to
efficiently find by simply randomly sampling the input around a given example.
Already, a variety of recent state of the art computer vision models employ
input deformations during training for increasing the robustness and
convergence speed of the models [9, 13]. These deformations are, however,
statistically inefficient, for a given example: they are highly correlated and
are drawn from the same distribution throughout the entire training of the
model. We propose a scheme to make this process adaptive in a way that
exploits the model and its deficiencies in modeling the local space around the
training data.

We make the connection with hard-negative mining explicitly, as it is close in
spirit: hard-negative mining, in computer vision, consists of identifying
training set examples (or portions thereof) which are given low probabilities
by the model, but which should be high probability instead, cf. [5]. The
training set distribution is then changed to emphasize such hard negatives and
a further round of model training is performed. As shall be described, the
optimization problem proposed in this work can also be used in a constructive
way, similar to the hard-negative mining principle.

###  4.1 Formal description

We denote by $f:\mathbb{R}^{m}\longrightarrow\\{1\dots k\\}$ a classifier
mapping image pixel value vectors to a discrete label set. We also assume that
$f$ has an associated continuous loss function denoted by
$\textrm{loss}_{f}:\mathbb{R}^{m}\times\\{1\dots
k\\}\longrightarrow\mathbb{R}^{+}$. For a given $x\in\mathbb{R}^{m}$ image and
target label $l\in\\{1\dots k\\}$, we aim to solve the following box-
constrained optimization problem:

  * •

Minimize $\|r\|_{2}$ subject to:

    1. 1.

$f(x+r)=l$

    2. 2.

$x+r\in[0,1]^{m}$

The minimizer $r$ might not be unique, but we denote one such $x+r$ for an
arbitrarily chosen minimizer by $D(x,l)$. Informally, $x+r$ is the closest
image to $x$ classified as $l$ by $f$. Obviously, $D(x,f(x))=f(x)$, so this
task is non-trivial only if $f(x)\neq l$. In general, the exact computation of
$D(x,l)$ is a hard problem, so we approximate it by using a box-constrained
L-BFGS. Concretely, we find an approximation of $D(x,l)$ by performing line-
search to find the minimum $c>0$ for which the minimizer $r$ of the following
problem satisfies $f(x+r)=l$.

  * •

Minimize $c|r|+\textrm{loss}_{f}(x+r,l)$ subject to $x+r\in[0,1]^{m}$

This penalty function method would yield the exact solution for $D(X,l)$ in
the case of convex losses, however neural networks are non-convex in general,
so we end up with an approximation in this case.

###  4.2 Experimental results

Figure 5: Adversarial examples generated for AlexNet [9].(Left) is a correctly
predicted sample, (center) difference between correct image, and image
predicted incorrectly magnified by 10x (values shifted by 128 and clamped),
(right) adversarial example. All images in the right column are predicted to
be an “ostrich, Struthio camelus”. Average distortion based on 64 examples is
0.006508. Plase refer to http://goo.gl/huaGPb for full resolution images. The
examples are strictly randomly chosen. There is not any postselection
involved.

Figure 6: Adversarial examples for QuocNet [10]. A binary car classifier was
trained on top of the last layer features without fine-tuning. The randomly
chosen examples on the left are recognized correctly as cars, while the images
in the middle are not recognized. The rightmost column is the magnified
absolute value of the difference between the two images.

Our “minimimum distortion” function $D$ has the following intriguing
properties which we will support by informal evidence and quantitative
experiments in this section:

  1. 1.

For all the networks we studied (MNIST, QuocNet [10], AlexNet [9]), for each
sample, we have always managed to generate very close, visually hard to
distinguish, adversarial examples that are misclassified by the original
network (see figure 5 and http://goo.gl/huaGPb for examples).

  2. 2.

Cross model generalization: a relatively large fraction of examples will be
misclassified by networks trained from scratch with different hyper-parameters
(number of layers, regularization or initial weights).

  3. 3.

Cross training-set generalization a relatively large fraction of examples will
be misclassified by networks trained from scratch on a disjoint training set.

The above observations suggest that adversarial examples are somewhat
universal and not just the results of overfitting to a particular model or to
the specific selection of the training set. They also suggest that back-
feeding adversarial examples to training might improve generalization of the
resulting models. Our preliminary experiments have yielded positive evidence
on MNIST to support this hypothesis as well: We have successfully trained a
two layer 100-100-10 non-convolutional neural network with a test error below
$1.2\%$ by keeping a pool of adversarial examples a random subset of which is
continuously replaced by newly generated adversarial examples and which is
mixed into the original training set all the time. We used weight decay, but
no dropout for this network. For comparison, a network of this size gets to
$1.6\%$ errors when regularized by weight decay alone and can be improved to
around $1.3\%$ by using carefully applied dropout. A subtle, but essential
detail is that we only got improvements by generating adversarial examples for
each layer outputs which were used to train all the layers above. The network
was trained in an alternating fashion, maintaining and updating a pool of
adversarial examples for each layer separately in addition to the original
training set. According to our initial observations, adversarial examples for
the higher layers seemed to be significantly more useful than those on the
input or lower layers. In our future work, we plan to compare these effects in
a systematic manner.

For space considerations, we just present results for a representative subset
(see Table 1) of the MNIST experiments we performed. The results presented
here are consistent with those on a larger variety of non-convolutional
models. For MNIST, we do not have results for convolutional models yet, but
our first qualitative experiments with AlexNet gives us reason to believe that
convolutional networks may behave similarly as well. Each of our models were
trained with L-BFGS until convergence. The first three models are linear
classifiers that work on the pixel level with various weight decay parameters
$\lambda$. All our examples use quadratic weight decay on the connection
weights: $\textrm{loss}_{decay}=\lambda\sum w_{i}^{2}/k$ added to the total
loss, where $k$ is the number of units in the layer. Three of our models are
simple linear (softmax) classifier without hidden units (FC10($\lambda$)). One
of them, FC10($1$), is trained with extremely high $\lambda=1$ in order to
test whether it is still possible to generate adversarial examples in this
extreme setting as well.Two other models are a simple sigmoidal neural network
with two hidden layers and a classifier. The last model, AE400-10, consists of
a single layer sparse autoencoder with sigmoid activations and 400 nodes with
a Softmax classifier. This network has been trained until it got very high
quality first layer filters and this layer was not fine-tuned. The last column
measures the minimum average pixel level distortion necessary to reach $0\%$
accuracy on the training set. The distortion is measure by
$\sqrt{\frac{\sum(x_{i}^{\prime}-x_{i})^{2}}{n}}$ between the original $x$ and
distorted $x^{\prime}$ images, where $n=784$ is the number of image pixels.
The pixel intensities are scaled to be in the range $[0,1]$.

In our first experiment, we generated a set of adversarial instances for a
given network and fed these examples for each other network to measure the
proportion of misclassified instances. The last column shows the average
minimum distortion that was necessary to reach 0% accuracy on the whole
training set. The experimental results are presented in Table 2. The columns
of Table 2 show the error (proportion of misclassified instances) on the so
distorted training sets. The last two rows are given for reference showing the
error induced when distorting by the given amounts of Gaussian noise. Note
that even the noise with stddev 0.1 is greater than the stddev of our
adversarial noise for all but one of the models. Figure 7 shows a
visualization of the generated adversarial instances for two of the networks
used in this experiment The general conclusion is that adversarial examples
tend to stay hard even for models trained with different hyperparameters.
Although the autoencoder based version seems most resilient to adversarial
examples, it is not fully immune either.

(a) Even columns: adversarial examples for a linear (FC) classifier
(stddev=0.06)

(b) Even columns: adversarial examples for a 200-200-10 sigmoid network
(stddev=0.063)

(c) Randomly distorted samples by Gaussian noise with stddev=1. Accuracy: 51%.

Figure 7: Adversarial examples for a randomly chosen subset of MNIST compared with randomly distorted examples. Odd columns correspond to original images, and even columns correspond to distorted counterparts. The adversarial examples generated for the specific model have accuracy 0% for the respective model. Note that while the randomly distorted examples are hardly readable, still they are classified correctly in half of the cases, while the adversarial examples are never classified correctly.  Model Name | Description | Training error | Test error | Av. min. distortion
---|---|---|---|---
FC10($10^{-4}$) | Softmax with $\lambda=10^{-4}$ | 6.7% | 7.4% | 0.062
FC10($10^{-2}$) | Softmax with $\lambda=10^{-2}$ | 10% | 9.4% | 0.1
FC10($1$) | Softmax with $\lambda=1$ | 21.2% | 20% | 0.14
FC100-100-10 | Sigmoid network $\lambda=10^{-5},10^{-5},10^{-6}$ | 0% | 1.64% | 0.058
FC200-200-10 | Sigmoid network $\lambda=10^{-5},10^{-5},10^{-6}$ | 0% | 1.54% | 0.065
AE400-10 | Autoencoder with Softmax $\lambda=10^{-6}$ | 0.57% | 1.9% | 0.086
Table 1: Tests of the generalization of adversarial instances on MNIST. | FC10($10^{-4}$) | FC10($10^{-2}$) | FC10($1$) | FC100-100-10 | FC200-200-10 | AE400-10 | Av. distortion
---|---|---|---|---|---|---|---
FC10($10^{-4}$) | 100% | 11.7% | 22.7% | 2% | 3.9% | 2.7% | 0.062
FC10($10^{-2}$) | 87.1% | 100% | 35.2% | 35.9% | 27.3% | 9.8% | 0.1
FC10($1$) | 71.9% | 76.2% | 100% | 48.1% | 47% | 34.4% | 0.14
FC100-100-10 | 28.9% | 13.7% | 21.1% | 100% | 6.6% | 2% | 0.058
FC200-200-10 | 38.2% | 14% | 23.8% | 20.3% | 100% | 2.7% | 0.065
AE400-10 | 23.4% | 16% | 24.8% | 9.4% | 6.6% | 100% | 0.086
Gaussian noise, stddev=0.1 | 5.0% | 10.1% | 18.3% | 0% | 0% | 0.8% | 0.1
Gaussian noise, stddev=0.3 | 15.6% | 11.3% | 22.7% | 5% | 4.3% | 3.1% | 0.3
Table 2: Cross-model generalization of adversarial examples. The columns of
the Tables show the error induced by distorted examples fed to the given
model. The last column shows average distortion wrt. original training set.

Still, this experiment leaves open the question of dependence over the
training set. Does the hardness of the generated examples rely solely on the
particular choice of our training set as a sample or does this effect
generalize even to models trained on completely different training sets?

Model | Error on $P_{1}$ | Error on $P_{2}$ | Error on Test | Min Av. Distortion
---|---|---|---|---
FC100-100-10: 100-100-10 trained on $P_{1}$ | 0% | 2.4% | 2% | 0.062
FC123-456-10: 123-456-10 trained on $P_{1}$ | 0% | 2.5% | 2.1% | 0.059
FC100-100-10’ trained on $P_{2}$ | 2.3% | 0% | 2.1% | 0.058
Table 3: Models trained to study cross-training-set generalization of the generated adversarial examples. Errors presented in Table correpond to original not-distorted data, to provide a baseline. | FC100-100-10 | FC123-456-10 | FC100-100-10’
---|---|---|---
Distorted for FC100-100-10 (av. stddev=0.062) | 100% | 26.2% | 5.9%
Distorted for FC123-456-10 (av. stddev=0.059) | 6.25% | 100% | 5.1%
Distorted for FC100-100-10’ (av. stddev=0.058) | 8.2% | 8.2% | 100%
Gaussian noise with stddev=$0.06$ | 2.2% | 2.6% | 2.4%
Distorted for FC100-100-10 amplified to stddev=$0.1$ | 100% | 98% | 43%
Distorted for FC123-456-10 amplified to stddev=$0.1$ | 96% | 100% | 22%
Distorted for FC100-100-10’ amplified to stddev=$0.1$ | 27% | 50% | 100%
Gaussian noise with stddev=$0.1$ | 2.6% | 2.8% | 2.7%
Table 4: Cross-training-set generalization error rate for the set of
adversarial examples generated for different models. The error induced by a
random distortion to the same examples is displayed in the last row.

To study cross-training-set generalization, we have partitioned the 60000
MNIST training images into two parts $P_{1}$ and $P_{2}$ of size 30000 each
and trained three non-convolutional networks with sigmoid activations on them:
Two, FC100-100-10 and FC123-456-10, on $P_{1}$ and FC100-100-10 on $P_{2}$.
The reason we trained two networks for $P_{1}$ is to study the cumulative
effect of changing the hypermarameters and the training sets at the same time.
Models FC100-100-10 and FC100-100-10 share the same hyperparameters: both of
them are 100-100-10 networks, while FC123-456-10 has different number of
hidden units. In this experiment, we were distorting the elements of the test
set rather than the training set. Table 3 summarizes the basic facts about
these models. After we generate adversarial examples with $100\%$ error rates
with minimum distortion for the test set, we feed these examples to the each
of the models. The error for each model is displayed in the corresponding
column of the upper part of Table 4. In the last experiment, we magnify the
effect of our distortion by using the examples
$x+0.1\frac{x^{\prime}-x}{\|x^{\prime}-x\|_{2}}$ rather than $x^{\prime}$.
This magnifies the distortion on average by 40%, from stddev $0.06$ to $0.1$.
The so distorted examples are fed back to each of the models and the error
rates are displayed in the lower part of Table 4. The intriguing conclusion is
that the adversarial examples remain hard for models trained even on a
disjoint training set, although their effectiveness decreases considerably.

###  4.3 Spectral Analysis of Unstability

The previous section showed examples of deep networks resulting from purely
supervised training which are unstable with respect to a peculiar form of
small perturbations. Independently of their generalisation properties across
networks and training sets, the adversarial examples show that there exist
small additive perturbations of the input (in Euclidean sense) that produce
large perturbations at the output of the last layer. This section describes a
simple procedure to measure and control the additive stability of the network
by measuring the spectrum of each rectified layer.

Mathematically, if $\phi(x)$ denotes the output of a network of $K$ layers
corresponding to input $x$ and trained parameters $W$, we write

$\phi(x)=\phi_{K}(\phi_{K-1}(\dots\phi_{1}(x;W_{1});W_{2})\dots;W_{K})~{},$

where $\phi_{k}$ denotes the operator mapping layer $k-1$ to layer $k$. The
unstability of $\phi(x)$ can be explained by inspecting the upper Lipschitz
constant of each layer $k=1\dots K$, defined as the constant $L_{k}>0$ such
that

$\forall\,x,\,r~{},~{}\|\phi_{k}(x;W_{k})-\phi_{k}(x+r;W_{k})\|\leq
L_{k}\|r\|~{}.$

The resulting network thus satsifies $\|\phi(x)-\phi(x+r)\|\leq L\|r\|$, with
$L=\prod_{k=1}^{K}L_{k}$.

A half-rectified layer (both convolutional or fully connected) is defined by
the mapping $\phi_{k}(x;W_{k},b_{k})=\max(0,W_{k}x+b_{k})$. Let $\|W\|$ denote
the operator norm of $W$ (i.e., its largest singular value). Since the non-
linearity $\rho(x)=\max(0,x)$ is contractive, i.e. satisfies
$\|\rho(x)-\rho(x+r)\|\leq\|r\|~{}$ for all $x,r$; it follows that

$\|\phi_{k}(x;W_{k})-\phi_{k}(x+r;W_{k})\|=\|\max(0,W_{k}x+b_{k})-\max(0,W_{k}(x+r)+b_{k})\|\leq\|W_{k}r\|\leq\|W_{k}\|\|r\|~{},$

and hence $L_{k}\leq\|W_{k}\|$. On the other hand, a max-pooling layer
$\phi_{k}$ is contractive:

$\forall\,x\,,\,r\,,~{}\|\phi_{k}(x)-\phi_{k}(x+r)\|\leq\|r\|~{},$

since its Jacobian is a projection onto a subset of the input coordinates and
hence does not expand the gradients. Finally, if $\phi_{k}$ is a contrast-
normalization layer

$\phi_{k}(x)=\frac{x}{\Big{(}\epsilon+\|x\|^{2}\Big{)}^{\gamma}}~{},$

one can verify that

$\forall\,x\,,\,r\,,~{}\|\phi_{k}(x)-\phi_{k}(x+r)\|\leq\epsilon^{-\gamma}\|r\|$

for $\gamma\in[0.5,1]$, which corresponds to most common operating regimes.

It results that a conservative measure of the unstability of the network can
be obtained by simply computing the operator norm of each fully connected and
convolutional layer. The fully connected case is trivial since the norm is
directly given by the largest singular value of the fully connected matrix.
Let us describe the convolutional case. If $W$ denotes a generic $4$-tensor,
implementing a convolutional layer with $C$ input features, $D$ output
features, support $N\times N$ and spatial stride $\Delta$,

$Wx=\left\\{\sum_{c=1}^{C}x_{c}\star
w_{c,d}(n_{1}\Delta,n_{2}\Delta)\,;d=1\,\dots,D\right\\}~{},$

where $x_{c}$ denotes the $c$-th input feature image, and $w_{c,d}$ is the
spatial kernel corresponding to input feature $c$ and output feature $d$, by
applying Parseval’s formula we obtain that its operator norm is given by

$\|W\|=\sup_{\xi\in[0,N\Delta^{-1})^{2}}\|A(\xi)\|~{},$ (1)

where $A(\xi)$ is a $D\times(C\cdot\Delta^{2})$ matrix whose rows are

$\forall~{}d=1\dots
D~{},~{}A(\xi)_{d}=\Big{(}\Delta^{-2}\widehat{w_{c,d}}(\xi+l\cdot
N\cdot\Delta^{-1})\,;\,c=1\dots C\,,\,l=(0\dots\Delta-1)^{2}\Big{)}~{},$

and $\widehat{w_{c,d}}$ is the 2-D Fourier transform of $w_{c,d}$:

$\widehat{w_{c,d}}(\xi)=\sum_{u\in[0,N)^{2}}w_{c,d}(u)e^{-2\pi
i(u\cdot\xi)/N^{2}}~{}.$

Layer | Size | Stride | Upper bound
---|---|---|---
Conv. $1$ | $3\times 11\times 11\times 96$ | $4$ | $2.75$
Conv. $2$ | $96\times 5\times 5\times 256$ | $1$ | $10$
Conv. $3$ | $256\times 3\times 3\times 384$ | $1$ | $7$
Conv. $4$ | $384\times 3\times 3\times 384$ | $1$ | $7.5$
Conv. $5$ | $384\times 3\times 3\times 256$ | $1$ | $11$
FC. 1 | $9216\times 4096$ | N/A | $3.12$
FC. 2 | $4096\times 4096$ | N/A | $4$
FC. 3 | $4096\times 1000$ | N/A | $4$
Table 5: Frame Bounds of each rectified layer of the network from [9].

Table 5 shows the upper Lipschitz bounds computed from the ImageNet deep
convolutional network of [9], using (1). It shows that instabilities can
appear as soon as in the first convolutional layer.

These results are consistent with the exsitence of blind spots constructed in
the previous section, but they don’t attempt to explain why these examples
generalize across different hyperparameters or training sets. We emphasize
that we compute upper bounds: large bounds do not automatically translate into
existence of adversarial examples; however, small bounds guarantee that no
such examples can appear. This suggests a simple regularization of the
parameters, consisting in penalizing each upper Lipschitz bound, which might
help improve the generalisation error of the networks.

##  5 Discussion

We demonstrated that deep neural networks have counter-intuitive properties
both with respect to the semantic meaning of individual units and with respect
to their discontinuities. The existence of the adversarial negatives appears
to be in contradiction with the network’s ability to achieve high
generalization performance. Indeed, if the network can generalize well, how
can it be confused by these adversarial negatives, which are indistinguishable
from the regular examples? Possible explanation is that the set of adversarial
negatives is of extremely low probability, and thus is never (or rarely)
observed in the test set, yet it is dense (much like the rational numbers),
and so it is found near every virtually every test case. However, we don’t
have a deep understanding of how often adversarial negatives appears, and thus
this issue should be addressed in a future research.

## References

  * [1] David Baehrens, Timon Schroeter, Stefan Harmeling, Motoaki Kawanabe, Katja Hansen, and Klaus-Robert Müller.  How to explain individual classification decisions.  The Journal of Machine Learning Research, 99:1803–1831, 2010.
  * [2] Yoshua Bengio.  Learning deep architectures for ai.  Foundations and trends® in Machine Learning, 2(1):1–127, 2009.
  * [3] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei.  Imagenet: A large-scale hierarchical image database.  In Computer Vision and Pattern Recognition, 2009. CVPR 2009. IEEE Conference on, pages 248–255. IEEE, 2009.
  * [4] Dumitru Erhan, Yoshua Bengio, Aaron Courville, and Pascal Vincent.  Visualizing higher-layer features of a deep network.  Technical Report 1341, University of Montreal, June 2009.  Also presented at the ICML 2009 Workshop on Learning Feature Hierarchies, Montréal, Canada.
  * [5] Pedro Felzenszwalb, David McAllester, and Deva Ramanan.  A discriminatively trained, multiscale, deformable part model.  In Computer Vision and Pattern Recognition, 2008. CVPR 2008. IEEE Conference on, pages 1–8. IEEE, 2008.
  * [6] Ross Girshick, Jeff Donahue, Trevor Darrell, and Jitendra Malik.  Rich feature hierarchies for accurate object detection and semantic segmentation.  arXiv preprint arXiv:1311.2524, 2013.
  * [7] Ian Goodfellow, Quoc Le, Andrew Saxe, Honglak Lee, and Andrew Y Ng.  Measuring invariances in deep networks.  Advances in neural information processing systems, 22:646–654, 2009\.
  * [8] Geoffrey E. Hinton, Li Deng, Dong Yu, George E. Dahl, Abdel rahman Mohamed, Navdeep Jaitly, Andrew Senior, Vincent Vanhoucke, Patrick Nguyen, Tara N. Sainath, and Brian Kingsbury.  Deep neural networks for acoustic modeling in speech recognition: The shared views of four research groups.  IEEE Signal Process. Mag., 29(6):82–97, 2012.
  * [9] Alex Krizhevsky, Ilya Sutskever, and Geoff Hinton.  Imagenet classification with deep convolutional neural networks.  In Advances in Neural Information Processing Systems 25, pages 1106–1114, 2012.
  * [10] Quoc V Le, Marc’Aurelio Ranzato, Rajat Monga, Matthieu Devin, Kai Chen, Greg S Corrado, Jeff Dean, and Andrew Y Ng.  Building high-level features using large scale unsupervised learning.  arXiv preprint arXiv:1112.6209, 2011.
  * [11] Yann LeCun and Corinna Cortes.  The mnist database of handwritten digits, 1998.
  * [12] Tomas Mikolov, Kai Chen, Greg Corrado, and Jeffrey Dean.  Efficient estimation of word representations in vector space.  arXiv preprint arXiv:1301.3781, 2013.
  * [13] Matthew D Zeiler and Rob Fergus.  Visualizing and understanding convolutional neural networks.  arXiv preprint arXiv:1311.2901, 2013.