Stable Architectures for Deep Neural Networks

TL;DR

This paper introduces ODE-inspired forward propagation techniques to improve the stability and training of very deep neural networks, addressing issues like exploding and vanishing gradients.

Contribution

It proposes a novel interpretation of deep learning as a nonlinear dynamical system and develops new architectures based on stability analysis to enhance deep network training.

Findings

New architectures improve stability in very deep networks

Numerical experiments show competitiveness with state-of-the-art models

Strategies effectively address exploding and vanishing gradients

Abstract

Deep neural networks have become invaluable tools for supervised machine learning, e.g., classification of text or images. While often offering superior results over traditional techniques and successfully expressing complicated patterns in data, deep architectures are known to be challenging to design and train such that they generalize well to new data. Important issues with deep architectures are numerical instabilities in derivative-based learning algorithms commonly called exploding or vanishing gradients. In this paper we propose new forward propagation techniques inspired by systems of Ordinary Differential Equations (ODE) that overcome this challenge and lead to well-posed learning problems for arbitrarily deep networks. The backbone of our approach is our interpretation of deep learning as a parameter estimation problem of nonlinear dynamical systems. Given this formulation, we analyze stability and well-posedness of deep learning and use this new understanding to develop new network architectures. We relate the exploding and vanishing gradient phenomenon to the stability of the discrete ODE and present several strategies for stabilizing deep learning for very deep networks. While our new architectures restrict the solution space, several numerical experiments show their competitiveness with state-of-the-art networks.

Tables

Table 1. Table 1: Multi-level training result for MNIST data set (see Section  6.4 ). We show the Training Error (TE) and Validation Error (VE) for the standard and antisymmetric ResNet and the Hamiltonian inspired network with Verlet forward propagation for increasing number of layers in the network. The different forward propagation methods yield comparable results with the antisymmetric ResNet giving slightly lower validation errors at each level.

	ResNet		anti symmetric ResNet		Hamiltonian Verlet
layers	TE	VE	TE	VE	TE	VE
4	0.96%	1.71%	1.13%	1.70%	1.49%	2.29%
8	0.80%	1.59%	0.92%	1.46%	0.82%	1.60%
16	0.73%	1.53%	0.91%	1.38%	0.35%	1.58%

Equations

$${\mathbf{Y}}_{j+1}={\mathbf{Y}}_j+h\sigma ({\mathbf{Y}}_j{\mathbf{K}}_j+b_j)\text{ for }j=0,\dots ,N-1.$$

$${\sigma }_{\mathrm{ht}}(\mathbf{Y})=\tanh (\mathbf{Y})\text{ and }{\sigma }_{\mathrm{ReLU}}(\mathbf{Y})=\max (0,\mathbf{Y}).$$

$$\mathbf{h}(\mathbf{x})=\exp (\mathbf{x})./(1+\exp (\mathbf{x})),$$

$$\mathbf{h}(\mathbf{X})=\exp (\mathbf{X})./(\exp (\mathbf{X}){\mathbf{e}}_m).$$

$$\min$$

$$s.t.$$

$$\frac{1}{s}S(\mathbf{h}({\mathbf{Y}}_N\mathbf{W}+{\mathbf{e}}_s{\mu }^{\top }),\mathbf{C})$$

$$\dot{\mathbf{y}}(t)=\sigma ({\mathbf{K}}^{\top }(t)\mathbf{y}(t)+b(t)),\text{ with }\mathbf{y}(0)={\mathbf{y}}_0,$$

$$\underset{i=1,2,\dots ,n}{\max }Re({\lambda }_i(\mathbf{J}(t)))\le 0,\forall t\in [0,T],$$

$$\mathbf{J}(t)$$

$$\underset{i=1,2,\dots ,n}{\max }Re({\lambda }_i(\mathbf{K}(t)))\le 0,\forall t\in [0,T].$$

$$\underset{i=1,2,\dots ,n}{\max }|1+h{\lambda }_i({\mathbf{J}}_j)|\le 1,\forall j=0,1,\dots ,N-1.$$

$${\bf K}_{+}=\left(\begin{array}[]{rr}2&-2\\ 0&2\end{array}\right),\;{\bf K}_{-}=\left(\begin{array}[]{rr}-2&0\\ 2&-2\end{array}\right),\;{\bf K}_{0}=\left(\begin{array}[]{rr}0&-1\\ 1&0\end{array}\right),$$

$$Re\left({\lambda }_i(\mathbf{K}(t))\right)\approx 0,\forall i=1,2,\dots ,n,t\in [0,T].$$

$${\mathbf{Y}}_{j+1}={\mathbf{Y}}_j+h\sigma \left(\frac{1}{2}{\mathbf{Y}}_j\left({\mathbf{K}}_j-{\mathbf{K}}_j^{\top }-\gamma \mathbf{I}\right)+b_j\right),j=0,\dots ,N-1,$$

$$\dot{\mathbf{y}}(t)=\sigma \left(\frac{1}{2}(\mathbf{K}(t)-\mathbf{K}(t{)}^{\top }-\gamma \mathbf{I})\mathbf{y}(t)+b(t)\right),\forall t\in [0,T].$$

$$\dot{\mathbf{y}}(t)=-{\nabla }_{\mathbf{z}}H(\mathbf{y},\mathbf{z},t)\text{ and }\dot{\mathbf{z}}(t)={\nabla }_{\mathbf{y}}H(\mathbf{y},\mathbf{z},t),\forall t\in [0,T].$$

$$H(\mathbf{y},\mathbf{z})=\frac{1}{2}{\mathbf{z}}^{\top }\mathbf{z}+f(\mathbf{y}),$$

$$\dot{\mathbf{z}}(t)=-{\nabla }_{\mathbf{y}}f(\mathbf{y}(t)),\dot{\mathbf{y}}(t)=\mathbf{z}(t),\forall t\in [0,T].$$

$$\ddot{\mathbf{y}}(t)={\nabla }_{\mathbf{y}}f(\mathbf{y}(t)),\forall t\in [0,T].$$

$$\ddot{\mathbf{y}}(t)=\sigma ({\mathbf{K}}^{\top }(t)\mathbf{y}(t)+b(t)),\mathbf{y}(0)={\mathbf{y}}_0,\dot{\mathbf{y}}(0)={\dot{\mathbf{y}}}_0,$$

$$\mathbf{K}(\mathbf{C})=-{\mathbf{C}}^{\top }\mathbf{C},\text{ for some }\mathbf{C}\in \text{mathds}{R}^{n\times n}.$$

$$\dot{\mathbf{y}}(t)=\sigma (\mathbf{K}(t)\mathbf{z}(t)+b(t))\text{ and }\dot{\mathbf{z}}(t)=-\sigma (\mathbf{K}(t{)}^{\top }\mathbf{y}(t)+b(t)).$$

$${\frac{\partial}{\partial t}}\left(\begin{array}[]{r}{\bf y}\\ {\bf z}\end{array}\right)(t)=\sigma\left(\left(\begin{array}[]{rr}0&{\bf K}(t)\\ -{\bf K}(t)^{\top}&0\\ \end{array}\right)\left(\begin{array}[]{r}{\bf y}\\ {\bf z}\end{array}\right)(t)+b(t)\right),\;\left(\begin{array}[]{r}{\bf y}\\ {\bf z}\end{array}\right)(0)=\left(\begin{array}[]{r}{\bf y}_{0}\\ 0\end{array}\right).$$

$${\bf y}_{j+1}=\left\{\begin{array}[]{ll}2{\bf y}_{j}+h^{2}\sigma({\bf K}_{j}{\bf y}_{j}+b_{j}),&j=0\\ 2{\bf y}_{j}-{\bf y}_{j-1}+h^{2}\sigma({\bf K}_{j}{\bf y}_{j}+b_{j}),&j=1,2,\ldots,N-1\end{array}\right..$$

$${\mathbf{z}}_{j+\frac{1}{2}}={\mathbf{z}}_{j-\frac{1}{2}}-h\sigma ({\mathbf{K}}_j^{\top }{\mathbf{y}}_j+b_j)\text{ and }{\mathbf{y}}_{j+1}={\mathbf{y}}_j+h\sigma ({\mathbf{K}}_j{\mathbf{z}}_{j+\frac{1}{2}}+b_j).$$

$${\bf K}_{v}=\left(\begin{array}[]{rr}2&1\\ -1&2\\ 0&1\end{array}\right).$$

$$\frac{\partial {\mathbf{z}}_{k+\frac{1}{2}}}{\partial {\mathbf{K}}_k}$$

$$\frac{\partial {\mathbf{y}}_{k+1}}{\partial {\mathbf{K}}_k}$$

$$\frac{\partial {\mathbf{z}}_{j+\frac{1}{2}}}{\partial {\mathbf{K}}_k}=$$

$$\frac{\partial {\mathbf{y}}_{j+1}}{\partial {\mathbf{K}}_k}=$$

$$\left(\begin{array}[]{rrrrrrr}{\bf I}&&&&&&\\ &{\bf I}&&&&&\\ -{\bf I}&{\bf B}_{k+1}&{\bf I}&&&&\\ &-{\bf I}&{\bf A}_{k+1}&{\bf I}&&&\\ &&\ddots&\ddots&\ddots&&\\ &&&-{\bf I}&{\bf A}_{N}&{\bf I}&\\ &&&&-{\bf I}&{\bf B}_{N}&{\bf I}\\ \end{array}\right)\left(\begin{array}[]{r}{\frac{\partial{\bf z}_{k+{\frac{1}{2}}}}{\partial{\bf K}_{k}}}\\ {\frac{\partial{\bf y}_{k+1}}{\partial{\bf K}_{k}}}\\ \vdots\\ \\ {\frac{\partial{\bf z}_{N+{\frac{1}{2}}}}{\partial{\bf K}_{k}}}\\ {\frac{\partial{\bf y}_{N+1}}{\partial{\bf K}_{k}}}\\ \end{array}\right)=\left(\begin{array}[]{c}{\bf C}_{1}\\ {\bf C}_{2}\\ 0\\ \vdots\\ \\ \\ 0\end{array}\right),$$

Full text

The backbone of our approach is our interpretation of deep learning as a parameter estimation problem of nonlinear dynamical systems. Given this formulation, we analyze stability and well-posedness of deep learning and use this new understanding to develop new network architectures. We relate the exploding and vanishing gradient phenomenon to the stability of the discrete ODE and present several strategies for stabilizing deep learning for very deep networks. While our new architectures restrict the solution space, several numerical experiments show their competitiveness with state-of-the-art networks.

Keywords. Machine Learning, Deep Neural Networks, Dynamic Inverse Problems, PDE-Constrained Optimization, Parameter Estimation, Image Classification.

1 Introduction

In this work, we propose new architectures for Deep Neural Networks (DNN) and exemplarily show their effectiveness for solving supervised Machine Learning (ML) problems; for a general overview about DNN and ML see, e.g., [40, 21, 1, 22] and reference therein. We consider the following classification problem: Assume we are given training data consisting of $s$ feature vectors, ${\bf y}_{1},\ldots,{\bf y}_{s}\in\mathds{R}^{n}$, and label vectors, ${\bf c}_{1},\ldots,{\bf c}_{s}\in\mathds{R}^{m}$, whose $k$th components represent the likelihood of an example belonging to class $k$. The goal is to learn a function that approximates the data-label relation on the training data and generalizes well to similar unlabeled data. Our goals in this work are to highlight the relation of the learning problem to dynamic inverse problems, analyze its stability and ill-posedness for commonly used architectures, and derive new architectures that alleviate some of these difficulties for arbitrarily deep architectures.

We are particularly interested in deep learning, i.e., machine learning using neural networks with many hidden layers. DNNs have been successful in supervised learning, particularly when the relationship between the data and the labels is highly nonlinear; see, e.g., [34, 30, 5, 35] and reference therein. Their depths (i.e., their number of layers) allow DNNs to express complex data-label relationships since each layer nonlinearly transforms the features and therefore effectively filters the information content.

Given the training data $({\bf y}_{1},{\bf c}_{1}),\ldots({\bf y}_{s},{\bf c}_{s})$, an inverse problem needs to be solved in order to train a given network architecture. This problem, also called the learning problem, aims at finding a parameterization of the DNN that explains the data-label relation and generalizes well to new unlabeled data. Clearly, using deeper network architectures increases the capacity of the network but also the dimensionality, and thus the computational complexity, of the parameter estimation problem. Additionally, more labeled data is required to calibrate very deep networks reliably. Therefore, despite the fact that neural networks have been used since the early 70’s, deep learning has only recently revolutionized many applications fueled by advances in computational hardware and availability of massive data sets.

Well-known sources of difficulty in deep learning are the dimensionality and non-convexity of the associated optimization problem. Traditionally, stochastic gradient descent methods have been used [43]. It has been observed that the performance of the DNN can be highly dependent on the choice of optimization algorithm and sample size [14, 11]. Furthermore, it has been noted that some optimization algorithms yield DNNs that generalize poorly to new unlabeled data [33].

Additional difficulties in deep learning stem from instabilities of the underlying forward model, most importantly, the propagation of features through the DNN. As has been shown in [11], the output of some networks can be unstable with respect to small perturbations in the original features. A related problem is the observation of vanishing or exploding gradients [6]. These results are unsettling since predictions made by networks with unstable forward propagation are very sensitive to small perturbations of the input features (as is common, e.g., in adversarial attacks), which may render the network useless in practice.

The main goal of this work is to gain new insight into the stability of the forward propagation and the well-posedness of the learning problem summarized in the following two questions:

Given a network architecture and parameters obtained by some optimization process, is the forward propagation problem well-posed? 2. 2.

Is the learning problem well-posed? In other words, given sufficient training are there parameters such that the DNN generalizes well or can generalization be improved by adding appropriate regularization?

The first question is important because, while it may be possible to fit the training data even for unstable forward propagation models, the trained network is unlikely to generalize. In other words, small deviations in the data, e.g., due to noise, may be drastically amplified by the forward propagation, resulting in incorrect labels. We show that the forward problem can be thought of as a discretization of an Ordinary Differential Equation (ODE). Therefore, the stability of the network corresponds to the stability of its underlying ODE. Based on this observation we develop stability criteria for a simplified version of the commonly used Residual Network (ResNet) architecture [27] and develop new network architectures that are ensured to be stable and lead to well-posed learning problems.

The paper is organized as follows. In Section 2 we give a brief mathematical derivation of the deep learning problem illustrated using the ResNet architecture. In Section 3 we analyze the stability of the simplified ResNet forward propagation and the well-posedness of the resulting learning problem. Our analysis and examples suggest stability criteria, and as a result, in Section 4, we propose three new stable architectures. In Section 5 we propose regularization functions favoring smoothness of parameters and multi-level ideas to initialize deep networks. In Section 6 we demonstrate the effectiveness of our new architectures on a number of small-scale model problems and explore their effectiveness for solving image classification problems. Finally, in Section 7 we summarize the paper.

2 Mathematical Formulation of the Deep Learning Problem

In this section we briefly describe three main ingredients of deep learning relevant to our work; for a comprehensive introduction see, e.g., [40, 21, 1, 22]. First, we outline forward propagation techniques, which transforms the input features in a nonlinear way to filter their information. Second, we describe classification, which predicts the class label probabilities using the features at the output layer (i.e., the output of the forward propagation). Finally, we formulate the learning problem, which aims at estimating parameters of the forward propagation and classification that approximate the data-label relation.

For notational convenience we stack the training features and labels row-wise into matrices ${\bf Y}_{0}=[{\bf y}_{1},{\bf y}_{2},\ldots,{\bf y}_{s}]^{\top}\in\mathds{R}^{s\times n}$ and ${\bf C}=[{\bf c}_{1},{\bf c}_{2},\ldots,{\bf c}_{s}]^{\top}\in\mathds{R}^{s\times m}$.

To exemplify our discussion of forward propagation we consider a simplified version of the Residual Neural Network (ResNet) [27] model that has been very successful in classifying images using deep network architectures; see [22] for other options. In ResNets, the forward propagation of the input values, ${\bf Y}_{0}\in\mathds{R}^{s\times n}$, through a network consisting of $N$ layers is given by

[TABLE]

The propagation in Eq. 2.1 is parametrized by the nonlinear activation function $\sigma:\mathds{R}^{s\times n}\to\mathds{R}^{s\times n}$ and affine transformations represented by their weights, ${\bf K}_{0},{\bf K}_{1},\ldots,{\bf K}_{N-1}\in\mathds{R}^{n\times n}$, and biases, $b_{0},b_{1},\ldots,b_{N-1}\in\mathds{R}$. We augmented the original formulation in [27] by the parameter $h>0$ in order to increase the stability of the forward propagation and allow for a continuous interpretation of the process; see also Section 3. The values ${\bf Y}_{1},\ldots,{\bf Y}_{N-1}$ are also called hidden layers and ${\bf Y}_{N}$ is called the output layer. The activation function is applied element-wise and is typically (piecewise) smooth and monotonically non-decreasing. As two commonly used examples, we consider the hyperbolic tangent and the Rectified Linear Unit (ReLU) activations

[TABLE]

The class label probabilities are predicted using the values at the output layers, ${\bf Y}_{N}$, a hypothesis function ${\bf h}({\bf Y}_{N}{\bf W}+{\bf e}_{s}\mu^{\top})$, and its associated weights, ${\bf W}\in\mathds{R}^{n\times m}$, and bias, $\mu\in\mathds{R}^{m}$. Here ${\bf e}_{k}\in\mathds{R}^{k}$ denotes the $k$-dimensional vector of all ones. For Bernoulli variables (i.e., ${\bf C}\in\{0,1\}^{s\times m}$) it is natural to consider the logistic regression function

[TABLE]

where the exponential and the division are applied element-wise. For Multinomial distributions we use the softmax function

[TABLE]

The learning problem aims at estimating the parameters of the forward propagation (i.e., ${\bf K}_{0},{\bf K}_{1},\ldots,{\bf K}_{N-1}$ and $b_{0},b_{1},\ldots,b_{N-1}$) and the classifier (${\bf W}$ and $\mu$) so that the DNN accurately approximates the data-label relation for the training data and generalizes to new unlabeled data. As we show below, the learning problem can be cast as a dynamic inverse problem, which provides new opportunities for applying theoretical and computational techniques from parameter estimation to deep learning problems. We phrase learning as an optimization problem

[TABLE]

where the loss function $S$ is convex in its first argument and measures the quality of the predicted class label probabilities, the convex regularizer $R$ penalizes undesirable (e.g., highly oscillatory) parameters, and the parameter $\alpha>0$ balances between minimizing the data fit and regularity of the parameters. A simple example for a loss function is the sum-of-squared difference function $S({\bf C}_{\rm pred},{\bf C})=\frac{1}{2}\|{\bf C}_{\rm pred}-{\bf C}\|^{2}_{F}$. Since our numerical experiments deal with classification we use cross entropy loss functions. Choosing an ”optimal” regularizer, $R$, and regularization parameter, $\alpha$, is both crucial and nontrivial. Commonly Tikhonov regularization [19, 26, 48], also referred to as weight decay, has been used [23], although, other possibilities enforcing sparsity or other structure have been proposed [41]. We introduce novel regularization functions in Section 5. For simplicity, we assume that a suitable value of $\alpha>0$ is chosen by the user or that it is done dynamically as suggested in [12].

There are numerous approaches to solving the learning problem. In this work we use a simple block coordinate descent method to demonstrate the properties of the forward propagation. Our method alternates between updating the parameters of the classifier $({\bf W},\mu)$ fixing the current value of the propagated features ${\bf Y}_{N}$ and then updating the parameters of the forward propagation while keeping the updated weights of the classifier fixed. The first problem is typically convex and the latter problem is generally non-convex due to the forward propagation process. Both steps are based on subsampling the training data. To this end note that most common loss functions can be written as a sum over all examples, i.e.,

[TABLE]

where $\mathcal{T}\subset\{1,2,\ldots,s\}$ is a randomly chosen set updated in each iteration of the block coordinate descent method. The size of the batches is a parameter in our algorithm whose choice depends on the size and complexity of the problem and resource considerations. In the following, we assume the sample size is constant; for adaptive selection of sample size see, e.g., [14]. In each iteration of the block coordinate descent scheme, our algorithm approximately solves the resulting classification problem using an Newton-PCG method, i.e., an inexact Newton method that uses a Preconditioned Conjugate Gradient (PCG) method to determine a search direction (see, e.g., [42, 29, 47]). Subsequently, the weights of the forward propagation are updated using a Gauss-Newton-PCG method. Note that gradient and approximated Hessian computations require matrix-vector products with the derivative matrices of the values of the output layer, ${\bf Y}_{N}$, with respect to ${\bf K}_{0,\ldots,N-1}$ and $b_{0,\ldots,N-1}$. The matrix-vector products with the derivative matrix can be computed without its explicit construction through forward and backward propagation, respectively; see also [38]. However, this requires storing (or recomputing) the values at the output layers; see, for instance, Section 4.4 for derivative computation. Therefore, as also suggested in [13], we subsample $\mathcal{T}$ further to reduce the cost of the Hessian matrix-vector products in our PCG scheme.

Our implementation includes computing the validation error in each iteration of the block coordinate descent method. The final output of our algorithm are the parameters that achieve the lowest validation error.

3 Stability and well-posedness of the forward propagation

In this section we analyze the stability of the ResNet forward problem 2.1 and illustrate why some choices of transformation weights may generate instabilities or prohibit effective learning altogether.

It is well-known that any parameter estimation problem requires a well-posed forward problem, i.e., a problem whose output is continuous with respect to its input. For example, practical image classification algorithms need to be robust against noisy or slightly shifted input images. Ill-posedness of the forward problem implies that even if the estimated parameters lead to a small training error they will probably fail or will do poorly on a perturbation of that data. In other words a network whose forward propagation is ill-posed will generalize poorly. Thus, well-posedness of the forward propagation is a necessary condition to obtain DNNs that generalize well.

The following discussion also gives new perspectives on two well-known phenomena in the deep learning community: Vanishing and exploding gradients; see, e.g., [6]. These phenomena refer to the gradient of the objective function in Eq. 2 and pose a severe challenge for very deep architectures. Note that the gradient represents the sensitivity of the output with respect to a perturbation in the input. Thus, an exploding gradient implies that the output is unstable with respect to the input. Similarly a vanishing gradient implies that the output is insensitive with respect to the input. Clearly both cases prohibit effective training, but more importantly, may not provide DNNs that generalize well.

To understand the phenomena we consider a simplified version of the forward propagation in ResNets given in Eq. 2.1. As pointed out in [25] the forward propagation can be seen as an explicit Euler discretization of the nonlinear Ordinary Differential Equation (ODE)

[TABLE]

over a time interval $t=[0,T]$. The final time $T>0$ and the magnitude of ${\bf K}(t)$ control the depth of the network. The ODE is stable if ${\bf K}(t)$ is changing sufficiently slow and

[TABLE]

where $Re(\cdot)$ denotes the real part, $\lambda_{i}({\bf J}(t))$ is the $i$th eigenvalue of the Jacobian of the right hand side in Eq. 3.5, denoted by ${\bf J}(t)\in\mathds{R}^{n\times n}$. A more accurate statement that uses kinematic eigenvalues of the Jacobian can be found in [3]. Here, the Jacobian is

[TABLE]

Since the activation function $\sigma$ is typically monotonically non-decreasing, i.e., $\sigma^{\prime}(\cdot)\geq 0$, Eq. 3.6 is satisfied if ${\bf K}$ changes sufficiently slowly and

[TABLE]

Controlling the smoothness of ${\bf K}$ can be done by regularization as described in Section 5. To ensure the stability of the overall discrete forward propagation, we also require the discrete version of the ODE to have a sufficiently small $h$ as summarized in the following well-known lemma.

Lemma 1** **(Stability of Forward Euler Method)

The forward propagation in Eq. 2.1 is stable if

[TABLE]

Proof 1

See, e.g., [2, 4] for the proof of stability criteria for the forward Euler method.

The above discussion suggests that the stability of the continuous forward propagation in Eq. 3.8 and its discrete analog Eq. 3.9 need to be added to the optimization problem 2 as constraints. Otherwise, one may obtain some transformation weights, ${\bf K}_{0},{\bf K}_{1},\ldots,{\bf K}_{N-1}$, that may fit the training data but generate an unstable process. As discussed above these solutions cannot be expected to generalize well for other data.

We illustrate the stability issues in ResNet using a simple example.

Example 1** **(Stability of ResNet)

For $s=3$ and $n=2$ we consider the forward propagation through a ResNet as given by Eq. 2.1. We consider three networks consisting of $N=10$ identical layers, i.e., on each layer we use the activation $\sigma_{\rm ht}=\tanh$, $h=0.1$, $b=0$, and a constant weight matrix. To illustrate the impact of the eigenvalues of the weight matrix on the propagation, we consider three ResNets parameterized by

[TABLE]

where $\lambda_{1}({\bf K}_{+})=\lambda_{2}({\bf K}_{+})=2,\ \ \lambda_{1}({\bf K}_{-})=\lambda_{2}({\bf K}_{-})=-2$ and $\lambda_{1}({\bf K}_{0})=\imath,\ \ \lambda_{2}({\bf K}_{0})=-\imath$. We consider the feature vectors ${\bf y}_{1}=[0.1,0.1]^{\top},\ {\bf y}_{2}=-{\bf y}_{1},\ {\bf y}_{3}=[0,0.5]^{\top}$. After propagating the features through the layers, we illustrate the different propagations in Figure 1. We represent the values at the hidden layers as colored lines in the 2D plane where each color is associated with one feature vector. To highlight the differences in the dynamics in all three cases, we also depict the force field using black arrows in the background. This plot is often referred to as the phase plane diagram.

As can be seen in left subplot, the features diverge away from the origin and each other using ${\bf K}_{+}$. Note that ${\bf y}_{1}$ and ${\bf y}_{2}$, which are close together initially, depart into opposite directions. This clearly suggests an unstable forward propagation that cannot be expected to generalize well. In contrast to that, the center subplot shows that ${\bf K}_{-}$ yields an accumulation point at the origin. While the forward propagation satisfies Eq. 3.8 and is thus stable, the learning problem, which requires inversion of this process, is ill-posed. Finally, the antisymmetric matrix ${\bf K}_{0}$ yields a rotation in the feature space, which preserves the distances between the features and leads to a stable forward propagation and a well-posed learning problem.

Clearly the effects depicted here are more pronounced in deeper networks with more layers and/or larger values for $h$.

While the first case in Example 1 indicates that Eq. 3.8 is a necessary condition for successful learning, the second case suggests that it is not sufficient. In general, when $Re(\lambda({\bf K}(t)))<0$ for most times and the network is deep (e.g., long time integration) differences in the initial feature vectors decay. In other words, for all initial conditions ${\bf y}_{0}$ we have that ${\bf y}(t)\rightarrow 0$ as $t\to\infty$; compare with center plot in Figure 1. Hence, even though the forward problem is stable the inverse problem is highly ill-posed as it comparable to an inverse heat equation. In these situations, the gradients of the objective function in Eq. 2 will vanish since small changes in the kernels will have no noticeable impact on the values of the outputs.

If we are inspired from the propagation of signals through neurons, then a ResNet as in Eq. 2.4 with $\max_{i}Re\left(\lambda_{i}({\bf K})\right)>0$ for all $i$ consists of neurons that amplify the signal with no upper bound (which is not biological) and a ResNet with $\max_{i}Re\left(\lambda_{i}({\bf K})\right)\ll 0$ can be seen as a lossy network. A moderately lossy network may be advantageous when the input is noisy, since it tends to decay high order oscillations. However, having too much signal loss is clearly harmful as it also annihilates relevant differences in the input features.

In summary, our discussion of ResNets illustrates that stable forward propagation and well-posed learning problems can be obtained for deep networks when

[TABLE]

In this case, the forward propagation causes only moderate amplification or loss and thus even deep networks preserve features in the input data and allow for effective learning.

4 Stable Forward Propagation for DNNs

Motivated by the discussion in the previous section, we introduce three new forward propagation methods that are stable for arbitrarily deep neural networks and lead to well-posed learning problems. Our approaches, presented in Section 4.1 and 4.2, use different means to enforce Jacobians whose eigenvalues have very small real part; see Eq. 3.11. The methods in Section 4.2 are inspired by Hamiltonian systems and we propose leapfrog and Verlet integration techniques for forward propagation in Section 4.3. Finally, we compute the derivative of the Verlet method using back propagation [46] and discuss the relation between back propagation and the older and more general adjoint method [32] in Section 4.4.

4.1 Antisymmetric Weight Matrices

Perhaps the simplest way to obtain a stable forward propagation is to construct force fields whose Jacobians are antisymmetric. For example, consider the forward propagation

[TABLE]

where $\gamma\geq 0$ is a small constant and ${\bf I}$ is an identity matrix. The resulting forward propagation can be seen as a forward Euler discretization of the ODE

[TABLE]

Since ${\bf K}(t)-{\bf K}(t)^{\top}$ is anti-symmetric its eigenvalues are imaginary, which also holds for the Jacobian in Eq. 3.7. Thus, the continuous ResNet parameterized by the antisymmetric weight matrix is stable and preserves information given an appropriate integration technique and sufficiently small time steps.

To stabilize the discrete forward propagation, which is based on the forward Euler method, we have added diffusion to the system in Eq. 4.12. The amount of diffusion depends on the parameter $\gamma\geq 0$. While small values of $\gamma$ might improve the robustness against noise in the feature vectors, too large values can render the learning problem ill-posed. Alternatively, more advanced time integration methods can be used with $\gamma=0$, i.e., without adding diffusion.

4.2 Hamiltonian Inspired Neural Networks

Restricting the parameter space to antisymmetric kernels is only one way to obtain a stable forward propagation. Alternatively we can recast forward propagation as a Hamiltonian system, which has the structure

[TABLE]

The function $H:\mathds{R}^{n}\times\mathds{R}^{n}\times[0,T]\to\mathds{R}$ is the Hamiltonian. In our application Hamiltonian systems for forward propagation are attractive due to their property to conserve rather than increase or dissipate the energy of the system; see [2, Ch.6] for a general introduction.

The Hamiltonian function $H$ measures the energy of the system, which in the autonomous case gets conserved. A simple approach can be derived from the Hamiltonian

[TABLE]

where $f:\mathds{R}^{n}\to\mathds{R}$ is at least twice continuously differentiable. This leads to the ODE system

[TABLE]

Eliminating ${\bf z}$ we obtain a second order system

[TABLE]

Inspired by the continuous version of the ResNet forward propagation given in Eq. 3.5 we propose to use the following second order ODE

[TABLE]

where in the following we assume that $\dot{{\bf y}}_{0}=0$. The dynamical system in Eq. 4.13 is stable for all weight matrices with non-positive real eigenvalues, i.e., that satisfy Eq. 3.8. This constraint can either be added to the optimization problem 2 or, similar to the previous section, be enforced by design. The latter approach can be obtained, e.g., using the parametrization

[TABLE]

Note that in contrast to the antisymmetric model in Eq. 4.12, the negative definite model is based on a nonlinear parametrization which might lead to a more challenging optimization problem. Alternatively Eq. 3.8 can be added as a constraint to the optimization problem 2 requiring expensive eigenvalue computations.

The forward propagations Eqs. 2.1, 4.12, or 4.13 require additional constraints on ${\bf K}$ and are limited to square weight matrices. We therefore propose a network that is intrinsically stable and supports non-square weight matrices, i.e., a network whose forward propagation is stable independent of the choice of the weights. To this end, we introduce symmetry into the Hamiltonian system by defining

[TABLE]

Note that Eq. 4.15 can be seen as a special case of ResNet with an augmented variable ${\bf z}\in\mathds{R}^{m}$ and the associated ODE

[TABLE]

This system is stable regardless of the spectrum or size of the matrices ${\bf K}(t)$ since the overall matrix in the linear transformation is antisymmetric. To ensure stability of the discrete version the time step size needs to be sufficiently small; see Lemma 3.9.

4.3 Symplectic Forward Propagation

In this section we use symplectic integration techniques for solving the discrete versions of the Hamiltonian-inspired networks in Eq. 4.13 and Eq. 4.15. Symplectic methods have been shown to capture the long time features of Hamiltonian systems and thus provide a slightly different approach as compared to the forward Euler method used in ResNet 2.4; we refer to [2, Ch.6] for a detailed discussion. For the second-order ODE in Eq. 4.13 we recommend the conservative leapfrog discretization

[TABLE]

To discretize the augmented network in Eq. 4.15 we use a Verlet integration where for $j=0,1,\ldots,N-1$ we have

[TABLE]

Since both discretizations are symplectic, the respective Hamiltonians are preserved if the transformation weights are time invariant and the step size is sufficiently small.

We demonstrate the long time dynamics of the two new Hamiltonian-inspired forward propagation methods using a simple example with two-dimensional features and identical layers.

Example 2

Let $s=2$ and $n=2$ and consider the features ${\bf y}_{1}=[0.1,0.1]^{\top}$ and ${\bf y}_{2}=[-0.1,-0.1]^{\top}$. We consider networks with identical layers featuring $b=0$ and hyperbolic tangent as activation function. For the leapfrog integration we use the matrix ${\bf K}={\bf K}_{-}$ defined in Eq. 3.10 (recall that $\lambda_{1}({\bf K})=\lambda_{2}({\bf K})=-2$) and a step size of $h=1$. To illustrate the Verlet integration (which can handle non-square kernels) we use a time step size of $h=0.1$ and

[TABLE]

To expose the short term behavior we use $N=500$ layers and to demonstrate the non-trivial long term characteristics we use $N=5,000$ layers. We depict the phase plane diagrams for both networks and both depths in Figure 2. Even though we use identical layers with constant kernels the features neither explode nor vanish asymptotically.

We proposed three new ways of forward propagation to ensure stability for arbitrary numbers of layers. Among the three approaches, the Verlet method is the most flexible as it can handle non-square weighting matrices and does not require additional constraints on ${\bf K}$.

4.4 Derivatives of the Verlet Method

We now discuss the computation of derivatives for the Verlet method. The derivatives of the ResNet are standard and can be found, e.g., in [28]. Given the ResNet derivatives the antisymmetric model can be differentiated using chain rule and the time stepping. The leapfrog method can be treated in a similar way to usual ResNets and is therefore omitted.

Our presentation follows the standard approach used in machine learning known as * back propagation *[28]. However, we note that back propagation is a special case of the adjoint method [7] used in time-dependent optimal control; see, e.g.,[16, 9]. The adjoint method is more general than back propagation as it can also be used compute the gradient of less “standard” time stepping methods. Thus, we discuss the back propagation method in the context of the adjoint method as applied to the Verlet method. We start by the usual sensitivity equation and then discuss how to efficiently compute matrix-vector products.

Exemplarily we show how to differentiate the values at the output layer in Eq. 4.17 with respect to the weight matrix ${\bf K}_{k}$ at an arbitrary layer $k$. For simplicity we assume that ${\bf K}$ is square. Derivatives for non-square weight matrices and with respect to the bias, $b_{k}$, can be computed along the same lines. Clearly, the values of ${\bf z}_{j}$ and ${\bf y}_{j}$ at the hidden layers $j\leq k$ are independent of ${\bf K}_{k}$. Applying the chain rule to Eq. 4.17 we see that

[TABLE]

where $\otimes$ denotes the Kronecker product and ${\bf I}$ is the $n\times n$ identity matrix. We can now differentiate layer by layer, where for $k<j<N-1$ we obtain

[TABLE]

Combining equations 4.18, 4.19, and 4.20 the derivatives can be written compactly as a block linear system

[TABLE]

where

[TABLE]

The linear system in Eq. 4.22 is block triangular with identity matrices on its diagonal and thus can be solved explicitly using forward substitution. Such systems commonly arise in optimal control of time-dependent PDEs and in dynamic inverse problems. For more details involving notes on implementation see, e.g., [45].

In the optimal control literature the matrices ${\frac{\partial{\bf z}_{j+{\frac{1}{2}}}}{\partial{\bf K}_{k}}}$ and ${\frac{\partial{\bf y}_{j}}{\partial{\bf K}_{k}}}$ are often referred to as the sensitivity matrices. While the matrices can be dense and large, computing matrix-vector products can be done efficiently by forward propagation. To this end, the right hand side of Eq. 4.22 is multiplied by the vector and the linear system is then solved with a vector right hand side. The multiplication with the transpose, also known as the adjoint method, is done by solving the equation backwards. It is important to note that the back-propagation algorithm [46] is nothing but a particular implementation of the adjoint method discussed much earlier in [32].

5 Regularization

In this section we present derivative-based regularization functions and a multi-level learning approach to ensure smoothness of parameters in deep learning. By biasing the learning process towards smooth time dynamics we aim at improving the stability of the forward propagation (e.g., to ensure that ${\bf K}(t)$ changes sufficiently slow as assumed in Sec. 3) and ultimately the generalization. Intuitively, we cannot expect the network to generalize well in the absence of smoothness. While the presented regularization techniques are commonly employed in other areas of inverse problems [26, 48], e.g., imaging their application to deep learning is, to the best of our knowledge, rather novel.

5.1 Regularizing the Forward Propagation

The perhaps most common regularization strategy used in deep learning is referred to as * weight decay*, which is equivalent to standard Tikhonov regularization of the form

[TABLE]

where $\|\cdot\|_{F}$ denotes the Frobenius norm; see, e.g., [22]. While this regularization reduces the magnitude of the weights, it is not sensitive to rapidly changing weights between adjacent layers. To illustrate why this may be problematic, consider removing a single layer from a deep network. Since the network is deep, we should not expect large changes in the values of the output layer and thus similar classification errors. However, if this layer performs, e.g., a 90 degree rotation of the features while the adjacent layers keep the features unchanged, the effect will be dramatic.

Our interpretation of forward propagation as a system of non-autonomous ODEs and the stability analysis in Sec. 3 motivate that ${\bf K}$ should be smooth, or at least piecewise smooth in time. To this end we propose the following new regularization for the transformation weights

[TABLE]

This regularization favors weights that vary smoothly between adjacent layers. Furthermore, as we see in our numerical experiments, the regularization adds robustness to the process. We can easily add or subtract steps without significantly changing the final result, thus adding more generalizing power to our network.

5.2 Regularizing the Classification Weights

We propose using smoothness regularization on the classification weights for image classification problems, e.g., in Convolution Neural Networks (CNN); see; e.g., [22, Ch.9]. For motivation, let the examples in ${\bf Y}_{0}$ represent vectorized $n_{1}\times n_{2}\times n_{3}$ images. The network propagates the images in time and generates perturbed images of the same size. The hypothesis function predicts the class label probabilities based on affinely transformed output images, i.e.,

[TABLE]

The number of rows in ${\bf W}$ is $n=n_{1}\cdot n_{2}\cdot n_{3}$ and the operation ${\bf Y}_{N}(j,:)^{\top}{\bf W}(:,k)$ is a dot product between the $j$th output image and the classification weights for the $k$th class. Noting that the rows of ${\bf Y}_{N}$ and the columns of ${\bf W}$ can be interpreted as discrete images sets the stage for developing regularization approaches commonly used in image processing; see, e.g., [17].

To further motivate the importance of regularization, note that if the number of examples is much larger than the number of features (pixels in the image) and if there is no significant redundancy, finding the optimal ${\bf W}$ given ${\bf Y}_{N}$ and ${\bf C}$ is an over-determined problem. Otherwise the problem is ill-posed and there may be infinitely many weights that yield the same classification on the observed data. The risk in both cases is overfitting since the optimal classification weights can be highly irregular and generalize poorly. Thus, one way to enforce uniqueness and improve generalization is to regularize the weights.

A standard approach in machine learning also known as pooling aims at achieving this goal; see, e.g., [22]. The main idea of pooling is to coarsen the output images and thereby to decrease the dimensionality of the classification problem. The simplest approach is skipping, i.e., subsampling the image ${\bf Y}_{N}$, e.g., at every other pixel. Other common options are average pooling and max-pooling, where image patches are represented by their average and maximum value, respectively. From an inverse problems perspective, pooling can be thought of as a subspace regularization [26]. For example, average pooling is similar to requiring ${\bf W}$ to be constant over each image patch.

An alternative approach to regularization is to interpret the $j$th feature, ${\bf y}_{j}\in\mathds{R}^{n}$, and the $k$th classification weight, ${\bf w}_{k}\in\mathds{R}^{n}$, of the CNN as discretizations of image functions $y_{j}:\Omega\to\mathds{R}$ and $w_{k}:\Omega\to\mathds{R}$, respectively, where $\Omega\subset\mathds{R}^{2}$ is the image domain. Assuming $y_{j}$ is sufficiently regular (which is to be enforced by the regularization) we can see that the probability of the $j$th example belonging to the $k$th class can be viewed as

[TABLE]

where ${\rm vol}(\Omega)$ denotes the volume of the image domain. To obtain weights that are insensitive to small displacements of the images it is reasonable to favor spatially smooth parameters by using regularization of the form

[TABLE]

where ${\bf L}$ is a discretized differential operator. This regularization also embeds the optimal $w_{k}$ into a suitable function space. For example, using an image gradient ensures that $w_{k}$ is in the Sobolev space $H^{1}(\Omega,\mathds{R})$ [20]. Most importantly for our application at hand, derivative-based regularization yields smooth classification weights that, as we see next, can be interpreted by visual inspection.

5.3 Multi-level Learning

As another means of regularizing the problem, we exploit a multi-level learning strategy that gradually increases the number of layers in the network. Our idea is based on the continuous interpretation of the forward propagation in which the number of layers in the network corresponds to the number of discretization points. Our idea is closely related to cascadic multigrid methods [8] and ideas in image processing where multi-level strategies are commonly used to decrease the risk of being trapped in local minima; see, e.g., [39]. More details about multi-level methods in learning can be found in [25].

The basic idea is to first solve the learning problem using a network with only a few layers and then prolongate the estimated weights of the forward propagation to a network with twice as many layers. The prolongated weights are then used to initialize the optimization problem on the finer level. We repeat this process until a user-specified maximum number of layers is reached.

Besides realizing some obvious computational savings arising from the reduced size of the networks, the main motivation behind our approach is to obtain good starting guesses for the next level. This is key since, while deeper architectures offer more flexibility to model complicated data-label relation, they are in our experience more difficult to initialize. Additionally, the Gauss-Newton method used to estimate the parameters of the forward propagation benefits in general from good starting guesses.

6 Numerical Examples

In this section we present numerical examples for classification problems of varying level of difficulty. We begin with three examples aiming at learning classification functions in two variables, which allow us to easily assess and illustrate the performance of the DNNs. In Section 6.4 we show results for the MNIST data set [36, 37], which is a common benchmark problem in image classification.

6.1 Concentric Ellipses

As a first test we consider a small-scale test problem in two dimensions. The test data consists of 1,200 points that are evenly divided into two groups that form concentric ellipsoids; see left subplot in Figure 4. The original data is randomly divided into 1,000 training examples and 200 examples used for validation.

We train the original ResNet, the antisymmetric ResNet, and the Hamiltonian network with Verlet propagation using the block coordinate descent method and a multi-level strategy with $4,8,16,32,64,128,256,1024$ layers. Each block coordinate descent iteration consists of a classification step ($\leq$2 iterations of Newton-PCG with $\leq 2$ PCG iterations) and a Gauss-Newton-PCG step to update the propagation weights ($\leq$20 PCG iterations preconditioned by regularization operator). In all examples we use the logistic regression hypothesis function in Eq. 2.2 and $\tanh$ activation function. The final time is $T=20$ and the width of the network is $n=2$. To enforce smoothness of the propagation weights in time we employ the regularizer in Eq. 5.23 weighted by a factor of $\alpha=10^{-3}$. No regularization is used in the classification.

We show the performance of the multi-level scheme in the left subplot of Figure 3. For this simple data set, all forward propagation methods achieve an optimal validation accuracy of $100\%$ at some level. As to be expected, the validation accuracy increases with increasing depth of the network, which also results in more degrees of freedom. The results for the Verlet method at the final level are shown in Figure 4. The two steps of deep learning (propagation and classification) can be seen in the center plot. The propagation transforms the feature such that they can be linearly separated. The result of the learning process is a network that predicts the class for all points in the 2D plane, which is illustrated in the right subplot of Figure 4.

6.2 Swiss Roll

We consider another small-scale test problem that is inspired by the swiss roll example. The data is obtained by sampling the vector functions

[TABLE]

for $r\in[0,1]$ and $\theta\in[0,4\pi]$ at 513 points each. Every other point along the curve is removed from the data set and used for validation; see left subplot in Figure 5.

Given the data we use the same parameterization of the multi-level and block coordinate descent method as in the previous example. For all networks, except the Hamiltonian network with Verlet forward propagation, we increase the dimensionality of the input feature space to $n=4$. As before we use the $\tanh$ activation, logistic regression function in Eq. 2.2, and choose a final time of $T=20$. To enforce smoothness of the propagation weights in time we employ the regularizer in Eq. 5.23 weighted by a factor of $\alpha=5\cdot 10^{-3}$. No regularization is used in the classification.

We plot the validation accuracy for each network and the different steps of the multi-level strategy in the center subplot of Figure 3. The standard ResNet and the Hamiltonian NN with Verlet forward propagation both achieve an optimal validation accuracy for $N=1,024$. However, the convergence considerably faster for the Hamiltonian network that reaches the optimal accuracy with $N=32$ layers. We visualize the results obtained using the Verlet method in Figure 5.

6.3 Peaks

We propose a new challenging test problem for classification into multiple classes using the peaks function in MATLAB®, which reads,

[TABLE]

where $x\in[-3,3]^{2}$. The peaks function is smooth but has some nonlinearities and most importantly non-convex level sets. We discretize the function on a regular $256\times 256$ grid and divide the points into 5 different classes based on their function value. The points in each class are then randomly subsampled such that the training data approximately represents the volumes of the level sets. In our case the $s=5,000$ sample points are divided evenly into five classes of size 1,000. We illustrate the test data in the left subplot of Figure 6.

We randomly choose 20% of the data for validation and train the networks using the remaining examples. We use $\tanh$ as activation, the softmax hypothesis function in Eq. 2.3, and choose a final time of $T=5$. To enforce smoothness of the propagation weights in time we employ the regularizer in Eq. 5.23 weighted by a factor of $\alpha=5\cdot 10^{-6}$. No regularization is used in the classification. We use the same multi-level strategy and identical parameters in the block coordinate descent methods as in the previous examples.

We first train the standard ResNet and the antisymmetric variant where we use a width of $n=8$ by duplicating the features. The performance of both models is approximately the same; see right subplot in Figure 3. The optimal accuracy, achieved for $N=1,024$, is around 98.8%.

For the Hamiltonian network with Verlet forward propagation we use a narrower network containing only the original features (i.e., $n=2$ and no duplication). As in the previous experiments the accuracy generally improves with increasing depth (with one exception between levels 4 and 5). The optimal accuracy is obtained at the final level ($N=1,024$) with a validation error of $99.1\%$. We illustrate the results for the Verlet network in Figure 6. The center subplot shows how the forward propagation successfully rearranges the features such that they can be labeled using a linear classifier. The right subplot shows that the prediction function fits the training data, but also approximates the true level sets.

6.4 MNIST

We use the MNIST dataset of hand-written digits [36, 37] to illustrate the applicability of our methods in image classification problems. The data set consists of 60,000 labeled digital images of size $28\times 28$ showing hand written digits from 0 to 9. We randomly divide the data into 50,000 training and 10,000 validation examples. We train the network using the standard and antisymmetric ResNet and the first-order Hamiltonian network using Verlet integration.

We use a three-level multi-level strategy where the number of layers is $4,8,$ and $16$. In each step we use the block coordinate descent method to approximately solve the learning problem and prolongate the forward propagation parameters to the next level using piecewise linear interpolation. The width of the network is 6 (yielding $n=4,704$ features used in the classification) and we use $3\times 3$ convolution operators that are fully connected within a given layer to model the linear transformation matrices ${\bf K}_{0},\ldots,{\bf K}_{N-1}$. The final time is set to $T=6$. To compute the Gauss-Newton step we first compute the full gradient over all 50,000 examples and then randomly subsample 5,000 terms for Hessian computations in the PCG step. The maximum number of iterations is set to 20 at each layer.

Within each step of the block coordinate descent we solve the classification problem using at most 5 iterations of Newton-PCG with up to 10 inner iterations. We use a discretized Laplacian as a regularization operator in Eq. 5.25 and its shifted symmetric product as a preconditioner to favor smooth modes; see [15]. The regularization parameter used in the classification is $0.01$. The output values for all examples are used at the final layer and no pooling is performed.

Smooth time dynamics are enforced by penalizing the time derivatives as outlined in Section 5.1 with a manually calibrated regularization parameter of $\alpha=0.005$. The regularization operator is used to precondition the PCG iterations in the Gauss-Newton method.

To show the performance of the multi-level strategy, we summarize the training and validation errors at the respective levels in Table 1. The table shows similar performance for all forward propagation methods. Note that both the validation and training error are reduced for larger number of layers, but no overfitting is observed. In this experiment, the multi-level approach considerably simplified the initialization of the deep networks. For example, the initial parameters of the standard ResNet at the final level ($N=16$ layers) already gave a validation accuracy of 98.41%, which was improved to 98.47%. We illustrate the results of the antisymmetric ResNet, which yields slightly superior performance in terms of validation error, in Figure 7. The smoothness of the classification weights enforced by our regularizer can be seen in the right column.

7 Summary and conclusions

In this paper, we expose the relation between deep learning and dynamic inverse problems and lay the foundation for fruitful research at the interface of inverse problems and data science. Specifically, we propose new architectures for deep neural networks (DNN) that improve the stability of the forward propagation. Using derivative-based regularization, we also improve the well-posedness of the learning task. We also propose a multi-level strategy that simplifies the choice of initial parameters and thereby simplifies the training of very deep networks. Our forward propagation methods are inspired by Hamiltonian systems and motivated by a stability analysis that is based on a continuous formulation of the learning problem.

Our stability discussion in Section 3 provides new insights why ResNets [27], a commonly used forward propagation scheme in deep learning, can be unstable. Our result is based on a continuous interpretation of a simplified version of ResNets as a system of nonlinear ODEs and a standard stability argument. Using intuitive examples, we show the impact of the spectral properties of the transformation matrices on the stability of the forward propagation and state conditions under which gradients explode, vanish, or are stable. We also note that stability depends on the smoothness of the network parameters over time. Based on our findings and numerical experience, we argue that it is desirable to restrict ResNets to matrices that lead to Jacobians whose real parts of the spectrum are close to zero and penalize the time-derivative of network parameters through regularization. We also emphasize that the well-posedness of the learning problem in ResNets relies on sufficiently small time steps. Extending our analysis to more general versions of ResNet [27] is straightforward and will be investigated in future work.

Motivated by our theoretical considerations, we propose three new forward propagation methods in Section 4 that can be used to obtain well-posed learning problems for arbitrarily deep networks. The first one is a simple modification of ResNet that applies a linear transformation to the transformation matrices to annihilate the real parts of their eigenvalues. The other two methods are inspired by Hamiltonian systems and provide alternative ways to preserve in information in the forward propagation. The Hamiltonian networks require special integration techniques; see [2] for details on their treatment. The second-order forward propagation can be discretized using the leapfrog method, and our implementation of the first-order propagation uses the Verlet method. While the stability of the leapfrog method requires matrices with non-positive real eigenvalues, the Verlet method does not require any restrictions and yields the best performance in our numerical experiments.

Our approach to stabilizing the learning problem for very deep architectures differs significantly from existing approaches, most importantly batch normalization [31]. Batch normalization scales values at hidden layers to prevent vanishing or exploding gradients and has been shown to improve the efficiency of training deep networks. In contrast to that our approach does not modify the values of the features, but alleviates the need of normalization by constructing stable forward propagation methods.

In order to improve the generalization quality, improve stability, and simplify training of deep networks we also propose new regularization approaches that depend on our continuous formulation of the problem. We use derivative-based regularizers to favor smooth time dynamics and, for image classification problems, spatially smooth classification weights. To further regularize and simplify the initialization of the learning algorithm, we employ a multi-level learning strategy that gradually increases the depth of the network. In our experiments, this approach has been a simple yet effective way to obtain good initializations for the learning problems. Our regularization methods are commonly used in imaging science, e.g., image registration [39], however, to the best of our knowledge not commonly employed in deep learning. Our numerical examples show that approximately solving the regularized learning problem yields works that generalize well even when the number of network parameters exceeds the number of training features.

We illustrate our methods using three academic test problems with available ground-truth and the MNIST problem [37], which is a commonly used benchmark problem for image classification. Our experiments show that the proposed new architectures yield results that are competitive with the established ResNet architecture. This is particularly noteworthy for the proposed anti-symmetric ResNet, where we restrict the dimensionality of the search space to ensure stability.

By establishing a link between deep learning and dynamic inverse problems, we are positive that this work will stimulate further research by both communities. An important item of future work is investigating the impact of the proposed architectures on the performance of learning algorithms. Currently, stochastic gradient descent [43] is commonly used to train deep neural networks due to its computational efficiency and empirical evidence supporting its generalizations [10]. A specific question is if the improvements in the forward propagation and regularization proposed in this paper will lead to better generalization properties of subsampled second-order methods such as [13, 44]. Another thrust of future work is the development of automatic parameter selection strategies for deep learning based on the approaches presented, e.g., in [19, 26, 48, 24, 18]. A particular challenge in this application is the nontrivial relationship between the regularization parameters chosen for the classification and forward propagation parameters.

Acknowledgements

This work is supported in part by National Science Foundation award DMS 1522599 and by the NVIDIA Corporation through their donation of a Titan X GPU.