A FRAMEWORK OF DEEP NEURAL NETWORKS VIA THE SOLUTION OPERATOR OF PARTIAL DIFFERENTIAL EQUATIONS
Year: 2021

Abstract:
There is a close connection between deep neural networks (DNN) and partial differential equations (PDEs). Many DNN architectures can be modeled by PDEs and have been proposed in the literature. However, their neural network design space is restricted due to the specific form of PDEs, which prevents the design of more effective neural network structures. In this paper, we attempt to derive a general form of PDEs for the design of ResNet-like DNN. To achieve this goal, we first formulate DNN as an adjustment operator applied on the base classifier. Then based on several reasonable assumptions, we show the adjustment operator for ResNet-like DNN is the solution operator of PDEs. To show the effectiveness of the general form of PDEs, we show that several effective networks can be interpreted by our general form of PDEs and design a training method motivated by PDEs theory to train DNN models for better robustness and less chance of overfitting. Theoretically, we prove that the robustness of DNN trained with our method is certifiable and our training method reduces the generalization gap for DNN. Furthermore, we demonstrate that DNN trained with our method can achieve better generalization and is more resistant to adversarial perturbations than the baseline model.

heading: 1 INTRODUCTION
text: Deep neural networks (DNN) have achieved success in tasks such as image classification (Simonyan & Zisserman, 2015), speech recognition (Dahl et al., 2012), video analysis (Bo et al., 2011), and action recognition (Wang et al., 2016). Residual learning (He et al., 2016a;b) has revolutionized DNN architecture design, making it practical to train ultra-deep DNN, and has ability to avoid gradient vanishing. The idea of residual learning has motivated the development of powerful DNN such as WideResNet (Zagoruyko & Komodakis, 2016), ResNeXt (Xie et al., 2017), and DenseNet (Huang et al., 2017).
So how to understand Residual learning from a theoretical perspective? In (Weinan, 2017), the author formulated the forward propagation of ResNet as a dynamic system and in (Li & Shi, 2017; Haber et al., 2018; Li et al., 2017) bridge the connection between partial differential equations(PDE) and ResNet. In fact, ResNet is composed of a series of residual mapping. Each residual block is realized by adding a shortcut connection to connect the input and output of the original mapping (F ). Mathematically, the n-th residual mapping can be formulated as
xn+1 = xn + F(xn,wn) (1)
with x0 ∈ Rd be a data point. xn is the input to the n-th residual mapping and F(xn,wn) is a non-linear function parameterized bywn which can be learned by back propagating the training error. The prediction of the data point x0 can be expressed by ŷ = softmax(wFCxL). By introducing a temporal partition: tn = nT/L, for n = 0, 1, .., L with the time interval ∆t = T/L and let x(tn) = xn and w(tn) = wn. Without considering the dimensional consistency, ResNet can be formulated as :
x(tn+1) = x(tn) + F(x(tn),w(tn))
∆t ∆t, n = 0, 1, · · · , L,
which is the forward Euler discretization of the following equations with time step ∆t, that is,
dx(t) = v(x(t), t)dt, x(0) = x, t ∈ [0, T ], (2)
where v(x(tn), tn) = F(x(tn),w(tn))/∆t and L is the number of resdiual mapping in a ResNet architecture. The prediction can be expressed by ŷ = softmax(wFCx(T )). The equation 2 defines the characteristic curves of the following TE
∂u ∂t (x, t) = v(x, t)∇u(x, t), x ∈ Rd. (3)
So the prediction ŷ of ResNet can be seen as the solution u(x, T ) of following Transport equation(TE):{ ∂u ∂t (x, t) = v(x, t)∇u(x, t), x ∈ R d.
u(x, 0) = softmax(wFCx). (4)
In fact, not only ResNet, similar to the process above, many ResNet-like DNN can correspond to PDEs. For example, it has been shown (Lu et al., 2018) that PolyNet (Zhang et al., 2017) and FractalNet (Larsson et al., 2016) can be respectively considered numerical schemes to solve a PDEs system. Based on the link of PDEs and ResNet-like DNN, many theoretical results to understand ResNet-like DNN were proved. In (Chen et al., 2018), a continuous form for ResNet called Neural ODE was proposed. After that, the property of Neural ODE has been theoretically analyzed (Thorpe & van Gennip, 2018; Avelin & Nyström, 2021; Zhang et al., 2019a). Based on the stability theory of differential equations, many architectures has been designed to build more robust deep networks (Zhang et al., 2019b; Yan et al., 2019; Reshniak & Webster, 2021). Not only stability theory, inspired by the numerical scheme of differential equations, in (Zhu et al., 2018; Tao et al., 2018), various deep network architectures to improve testing accuracy have been proposed in success.
Despite the success of the application of PDEs in residual learning, these neural network design space in the literature mentioned above is restricted due to the specific form of PDEs, which prevents the design of more effective neural network structures. Existed works can not meet these requirements. So, a natural question is that can we find the general form of the PDEs corresponding to the ResNetlike DNN which can unify existed model as well as guide new effective DNN architecture design?
To achieve this goal, in this paper, we follow the setting of paragraph 2, the softmax model and ResNet are the initial value and the terminal value of PDEs equation 4 respectively, so we can view ResNet as an evolution from a softmax model. Without loss of generality, we indicate this evolutive relationship by an adjustment operator denoted by Tt, then we have
Tt : softmax(w FCx)→ ResNet(x).
If we denote the softmax model by f , then the ResNet-like DNN can be expressed as Tt(f) and Tt is the solution operator of PDEs equation 4. To obtain the general form of the PDEs corresponding to these ResNet-like DNN, we list a series of formal properties likely to be satisfied by the operator Tt, and we deduce an explicit formulation from them. We see that under reasonable assumptions, the new model u(x, t) = (Ttf)(x) is the solution of convection diffusion equations. Then, motivated by the theory of PDEs and the framework of our general form of PDEs, we design a training method for ResNet-like DNN models for better robustness and less overfitting risk. Our training method is summarized as follows. First, we train a ResNet-like DNN model using natural training. Because calculating the real solution Tt(f)(x) is difficult, we next train a new ResNet-like DNN denoted by g(x, t) to approximate Tt(f)(x). At last, we choose g(x, T ) as our final classifier.
We make the following contributions:
• Under reasonable assumptions, we build a PDEs framework for ResNet-like DNN with the help of operator Tt and show some effective models that can be modeled by our framework.
• We design a training method for DNN based on the theory of PDEs. This is an essential contribution of this paper due to certified adversarial robustness and theoretical analysis of generalization ability for our training method.
• We verify our training method on a synthetic dataset and the CIFAR-10/CIFAR-100 datasets, then demonstrate that a ResNet trained by our method is more resistant to multistep `2 and `∞ adversarial attacks than natural trained ResNet.

heading: 2 MAIN RESULTS
text: We show under several reasonable assumptions that the adjustment operator Tt is the solution operator of convection diffusion equations. Based on the regularity theory of PDEs, we design a training method for DNN models and theoretically illustrate its benefits for improving the robustness of models and reducing the generalization gap.

heading: 2.1 GENERAL FORM OF PDES
text: We denote the base DNN model by fw(x), which is uniformly continuous. We assume data points are distributed inD ⊂ Rd, and label function l(x) is defined inD but known only in training set SN ∈ D. Our method formulates the DNN as a continuous dynamic system, namely u(x, t) = Tt(f), t ∈ [0, T ]. To get the expression of the adjustment operator Tt, we assume it has some fundamental properties.
[Locality] (Tt(f)− Tt(g))(x) = o(t) as t→ 0+, for all f, g such that Dαf(x) = Dαg(x) for all |α| ≥ 0. [Comparison Principle] Tt(f) ≤ Tt(g) for all t ≥ 0 and f, g such that f ≤ g.
[Markov Property] Tt+s = Tt+s,s ◦Ts, for all s, t ≥ 0 and t+ s ≤ T . Tt+s,sdenotes the flow from time s to time t+ s.
[Spatial Regularity] There exist a positive constant C such that ‖Tt(τhf)− τh(Ttf)‖L∞ ≤ Ch for all h in Rd, t ≥ 0, where (τhf)(x) = f(x+ h) and ‖h‖2 = h. [Temporal Regularity] The following two properties hold: 1) For all t, s, t+ s ∈ [0, T ] and all g ∈ Q, there exist a constant C ≥ 0 depending only on Q such that ‖Tt+s,s(g)− g‖L∞ ≤ Ct, where Q is the subset of C∞b defined by
Q = {f ∈ C∞b ,∀n ≥ 0, ‖Dαf‖L∞ ≤ C for all |α| = n}.
2) For all t, s, t+ s ∈ [0, T ], there exist a constant C ≥ 0 such that ‖Tt+s,t(f)−Ts,0(f)‖L∞ ≤ Cst uniformly in s ∈ (0, T ]. [Ensemble Additivity] Tt(f + g) = Tt(f) + Tt(g), and if C is a constant, then Tt(C) = C.
To illustrate the reasonableness of these assumptions above, we give the following explanation. We first require that Tt can be expressed as a solution operator of PDEs and for small enough t, the value of Tt(f) at any point x is determined by the behavior of the base classifier f near x, which implies [Locality]. We then require the operator Tt is provided with an order-preserving property which means that no enhancement is made. Thus if the confidence level for an event of one base classifier g is higher than another base classifier f , this ordering is preserved, which implies [Comparison Principle]. Next we require Tt+s can be computed from Tt for any s ≤ t, and T0 is of course the identity. This is natural, since the future behaviour of the adjustment process at terminal time t+ s is likely to be deduced from the state at intermediate time t without any dependence upon the original model f , which implies [Markov Property]. We need function Tt(f) to be spatial stable. If we add a perturbation h to data point x, the output Tt(f)(x+ h) will not change a lot, which implies [Spatial Regularity]. We also need an additional property for giving some temporal stability to the operator Tt. Namely, in any time interval, the adjustment will not be rapid, which implies [Temporal Regularity]. Assume an ensemble classifier F consists of base classifiers f and g, and another ensemble classifier consists of base classifiers Tt(f) and Tt(g). We require the prediction Tt(F ) to equal Tt(f) + Tt(g), which implies [Ensemble Additivity]. Theorem 1. Under the above conditions, there exists Lipschitz continuous function v : Rd× [0, T ]→ Rd and Lipschitz continuous positive function σ : Rd × [0, T ]→ Rd×d such that for any bounded and uniformly continuous base model fw(x), u(x, t) = Tt(fw)(x) is the unique viscosity solution of the following convection diffusion equation:{
∂u(x,t) ∂t = v(x, t) · ∇u(x, t) + 1 2 ∑ i,j σi,j ∂2u ∂xi∂xj (x, t), u(x, 0) = fw(x), (5)
where x ∈ Rd, t ∈ [0, T ]. Here σi,j is the i, j-th element of matrix function σ(x, t).
We will provide the proof of Theorem 1 in Appendix A. In this subsection, we introduce a PDEs framework for DNN models. The framework is quite general, many existing models with residual connections can be modeled by our framework. As examples, we briefly list some of them in the next subsection.

heading: 2.2 SEVERAL MODELS
text: In this subsection, we will show Gaussian noise injection(Wang et al., 2020; Liu et al., 2019), randomized smoothing(Cohen et al., 2019; Li et al., 2018; Salman et al., 2019) and ResNet with stochastic dropping out the hidden state of residual block(Srivastava et al., 2014; Sun et al., 2018) can be interpreted by our framework.
Gaussian noise injection: Gaussian noise injection is an effective regularization mechanism for a DNN model. For a base ResNet with L residual mapping, the n-th residual mapping with Gaussian noise injected can be written as
xn+1 = xn + F(xn,wn) + aεn, εn ∼ N (0, I), where the parameter a is a noise coefficient. By introducing a temporal partition: tn = nT/L, for n = 0, 1, .., L with the time interval ∆t = T/L and let x(tn) = xn and w(tn) = wn. And let a = σ √ ∆t and F(xn,wn)/∆t = v(x, t). This noise injection technique in a discrete neural network can be viewed as the approximation of continuous dynamic
dx(t) = v(xn, t)dt+ σdB(t) (6)
where B(t) is multidimensional Brownian motion. The output of L-th residual mapping can be written as It̂o process equation 6 at terminal time T , x(T ). So, an ensemble prediction over all the possible sub-networks with shared parameters can be written as
ŷ = E(softmax(wFCx(T ))|x(0) = x0). (7) According to Feynman-Kac formula(Mao, 2007), equation equation 7 is known to solve the following convection diffusion equation{
∂u(x,t) ∂t = v(x, t) · ∇u+ 1 2σ 2∆u,x ∈ Rd, t ∈ [0, T ] u(x, 0) = softmax(wFCx)
Randomized Smoothing: Consider to transform a base classifier into a new smoothed classifier by adding Gaussian noise to the input when inference time. If we denote the base classifier by f(x) and denote the new smoothed classifier by g(x). Then f(x) and g(x) have the following relation:
g(x) = 1
N N∑ i=1 f(x+ εi) ≈ Eε∼N (0,I)[f(x+ ε)] (8)
According to Feynman-Kac formula, g(x) can be viewed as the solution of the following PDEs{ ∂u(x,t) ∂t = 1 2∆u, t ∈ [0, 1]
u(x, 0) = f(x). (9)
Especially, when f(x) is ResNet, the smoothed classifier g(x) = u(x, T + 1) can be expressed as ∂u(x,t) ∂t = 1 2∆u, t ∈ [T, T + 1] ∂u(x,t) ∂t = v(x, t) · ∇u(x, t),x ∈ R
d, t ∈ [0, T ] u(x, 0) = softmax(wFCx).
Dropout of Hidden Units: Consider the case that disables every hidden units independently from a Bernoulli distribution B(p) with p ∈ (0, 1) in each residual mapping
xn+1 = xn + F(xn,wn) zn p
= xn + F(xn,wn) + F(xn,wn) ( zn p − I),
where zn ∼ B(1, p) namely P(zn = 0) = 1 − p, P(zn = 1) = p and indicates the Hadamard product. If the number of the ensemble is large enough, according to Central Limit Theorem, we have
F(xn,wn) ( zn p − I) ≈ F(xn,wn) N (0, 1− p p ).
The similar way with Gaussian noise injection, the ensemble prediction ŷ can be viewed as the solution u(x, T ) of following equation:{
∂u(x,t) ∂t = v(x, t) · ∇u(x, t) + 1−p 2p ∑ i(v T v)i,i ∂2u ∂x2i (x, t),
u(x, 0) = softmax(wFCx) (10)
where x ∈ Rd, t ∈ [0, T ]. Remark 1. In fact, similar to dropout, shake-shake regularization (Gastaldi, 2017; Huang & Narayanan, 2018) and ResNet with stochastic depth (Huang et al., 2016) can also be formulated by our PDEs framework.
The viewpoint that forward propagation of a ResNet is the solving process of a TE enables us to interpret poor robustness and generalization ability as the irregularity of the solution. Morever we can also interpret the effectiveness of the models above-mentioned in improving generalization as well as robustness as the action of diffusion term. To mitigate this issue, based on our PDEs framework and the knowledge that the terminal value of convection diffusion equations (whose diffusion term is isotropic, namely σ(x, t) is σ2I) is more regular than initial value, we use the ResNet fw(x) as the base classifier and use the solution operator of the diffusion equation{
∂u(x,t) ∂t = v(x, t) · ∇u(x, t) + σ 2∆u,x ∈ Rd, t ∈ [0, T ] u(x, 0) = fw(x).
(11)
as the adjustment operator Tt to smooth the ResNet for improving the robustness and generalization ability. The case that diffusion term σ(x, t) is anisotropic (namely σ(x, t) 6= σ2I) is difficult for theoretical analysis and practical experiment and we take the case as our future work. In next subsection, we will illustrate the benefits of using the adjustment operator Tt from a theoretical perspective.

heading: 2.3 ROBUSTNESS GUARANTEE
text: Consider a classification problem from Rd to blue label class Y . Let G be a classifier defined by G(x) = arg maxi∈Y u
i(x, T ), where ui(x, T ) is the i-th element of u(x, T ). Suppose that the new DNN classifies x the most probable class cA is returned with probability pA = ucA(x, T ), and the “runner-up” class is returned with probability pB . Our main result of this subsection is to estimate around the data point x how large `p radius can make the classifier G(x) be robust, where 1 < p ≤ ∞. Theorem 2. Let the base model fw(x) be a compactly supported function and velocity v(x, t) ∈ C1(Rd × [0, T ]). Let the new DNN u(x, T ) be the solution of equation 11. Suppose cA ∈ Y and
ucA(x, T ) = pA ≥ pB = max i 6=cA ui(x, T ),
then G(x+ δ) = cA for all ‖δ‖p ≤ R, where
R = eσ 2T (pA − pB) 2d1/qeγT maxi ‖f iw(x)‖C1 ,
1/p+ 1/q = 1 and γ is a constant depending on∇v.
We provide the proof of Theorem 2 in Appendix B. According to Theorem 2, we can obtain the certified radius R is large when the diffusion coefficient σ2 and the probability of the top class is high. Not only robustness, comparing to the base classifier, but we will also show the generalization ability of model 11 will be better in next subsection.

heading: 2.4 GENERALIZATION GAP
text: For simplicity, we consider binary classification problems. Suppose SN ={(xi, yi)}Ni=1 is drawn from X × Y ⊂ D × {−1,+1} with X and Y being the input data and label spaces, respectively. Assume p(x) is the underlying distribution ofX , which is defined inD. We also assume label function Label : D → {−1,+1}, which is unknown. Let H ⊂ V X be the hypothesis class of the DNN model, where V is another space that might be different from Y . Let ` : V × Y → [0, B] be the loss function and B is a positive constant. Denote function class `H := {(x, y)→ `(h(x), y) : h ∈ H}. Rademacher complexity is one of the classic measures of generalization error. We first recap on the definition of Rademacher complexity.
Definition 1. (Ledoux & Talagrand, 2013)LetH : X → R be the space of real-valued functions on the space X . For a given sample SN = {(xi, yi)}Ni=1, the empirical Rademacher complexity ofH is defined as
RSN (H) := 1
N Eσ[sup
h∈H N∑ i=1 σih(xi)],
where σ1, σ2, · · · , σN are i.i.d. Rademacher random variables with P(σi = 1) = P(σi = −1) = 12 .
Rademacher complexity is a tool to bound the generalization error. The smaller the generalization gap is, the less overfitting the model is. For ∀xi ∈ Rd and constant c ≥ 0, we consider the following two function classes:
F := { f(x)|‖f‖C1 ≤ c } ,
and Gσ := { g(x) = u(x, 1)|∂u(x, t)
∂t = v(x, t)∇u(x, t) + σ2∆u, u(x, 0) = f(x) where f ∈ F } The function class F represents the continuous analogue of ResNet, and function class Gσ denotes F adjusted by operator Tt. Then we have the following theorem.
Theorem 3. Given a training set SN = {(xi, yi)}Ni=1. With probability at least 1− 1/(2N) we have
RSN (F) ≤ C√ N [(log c+ logN)1/2 + c]
RSN (Gσ) ≤ C√ N [(−σ2 + ln c+ lnN)1/2 + exp(−σ2)c]
Here C is a positive constant depend on D and dimension d.
We provide the proof of Theorem 3 in Appendix C. According to Theorem 3, we can obtain that for function class Gσ , the upper bound of Rademacher complexity is low when the diffusion coefficient σ2 is high. We next present a practical algorithm for real applications to verify our main results.

heading: 3 ALGORITHMS
text: Let SN = {(x1, y1), · · · , (xn, yn)} be training set, where xi is a data point and yi is the corresponding label. Then denote the ResNet with natural training by fw(x) whose parameters can be learned by minimized the loss lCE(yi, fw(xi)), where lCE is cross-entropy loss. . For convenience, we assume the velocity v(x, t) is zero and the diffusion term σ(x, t) is σ2I. We can rewrite equation 5 as {
∂u(x,t) ∂t = σ 2∆u(x, t), x ∈ Rd, t ∈ [0, T ] u(x, 0) = fw(x)
(12)
where σ2 is a hyperparameter. In this model, we set the natural trained ResNet as the initial value. Instead of calculating the solution of equation 12, we train a DNN gw(x, t) to approximate the real solution. The parameters of new DNN gw(x, t) are different from the initial ResNet fw(x). Therefore, to train the DNN model gw(x, t), we design loss terms on training set SN for the initial condition and differential equation.
We use lCE(gw(xi, 0), fw(xi)) to fit the initial condition of equation 12, where lCE is crossentropy loss. To improve the natural accuracy, we enhance the terminal value to fit the label, i.e., lCE(gw(xi, T ), yi). We can write the first loss term as
L1(gw, SN ) = 1
N N∑ i=1 [lCE(gw(xi, 0), fw(xi)) + l CE(gw(xi, T ), yi)]. (13)
We use mean square error loss to fit the differential equation and improve the robustness of our model,
L2(gw, SN ) = 1
N N∑ i=1 ∫ T 0 ( ∂gw ∂t (xi, s)− σ2∆gw(xi, s) )2 ds. (14)
According to the trapezoid formula, the integral term equation 14 is approximated by
L2(gw, SN ) ≈ 1
2N N∑ i=1 ∑ s=0,T ( ∂gw ∂t (xi, s)− σ2∆gw(xi, s) )2 . (15)
Using Taylor formula, the ∆gw(xi, s) can be approximated as
∆gw(xi, s) ≈ 1
M M∑ j=1 ( gw(xi + hεi,j , s) + gw(xi − hεi,j , s)− 2gw(xi, s) h2 ) . (16)
where {εi,j} is i.i.d and M is the average number. In practice, we choose the number of average M equals 1 to reduce the computation cost, and set h to 0.1 because training converges with difficulty when h is too small. We denote the right hand of equation 16 by ∆h,εgw(xi, s). Combining equation 15 and equation 16, the integral term equation 14 can be approximated by
Lh,M2 (gw, SN ) = 1
N N∑ i=1 ∑ s=0,T ( ∂gw ∂t (xi, 0)− σ2∆h,εgw(xi, s) )2 . (17)
Putting the two objective functions together, our training loss combines L1(gw, SN ) and Lh,M2 (gw, SN ), i.e., Loss(gw, SN ) = L1(gw, SN ) + λL h,M 2 (gw, SN ), where λ > 0 is a coefficient to trade off the two loss terms. In practice, we use stochastic gradient descent to optimize gw(x, t), and we use gw(x, T ) as the final model.

heading: 4 EXPERIMENTS
text: In this section, we will numerically verify that 1) can our method improve the robustness and reduce overfitting risk of DNN? 2) what is the advantage of our method over the Gaussian noise injection and randomized smoothing?

heading: 4.1 PRELIMINARIES
text: Datasets. In our experiments, we consider both Half-moon and CIFAR10/CIFAR100 (Krizhevsky et al., 2009) datasets. Half-moon dataset is a randomly generated 2d synthetic dataset in which we randomly generate 500 points and 1000 points with a standard deviation of 0.3 as training set and testing set, respectively. Both CIFAR10 and CIFAR100 contain 60K 32× 32 color images with 50K and 10K of them used for training and test, respectively. CIFAR10 and CIFAR100 contain 10 and 100 different classes, respectively.
Project gradient descent attack (PGD). To verify the efficiency of our method for improving the robustness of DNN, we consider the project gradient descent (PGD)(Madry et al., 2017) in all the experiments below. Before introducing PGD, we first introduce fast gradient sign method attack (FGSM). The FGSM attack searches the adversarial image x′ by maximizing the linearized loss function function L(x′, y) ≈ L(x, y) +∇L(x, y)T (x′ − x) with constraint ‖x′ − x‖p ≤ , where is the maximum perturbation and p = 2,∞. PGD iterates FGSM with step size α and clips the perturbed image to generate the enhanced adversarial attack,
`∞ : x(m) = Clipx, {x(m−1) + αsign(∇L(x(m−1), y))},
`2 : x(m) = Clipx, {x(m−1) + α( ∇L(x(m−1), y) ‖∇L(x(m−1), y)‖2) )},
where m = 1, · · · ,K, x(0) = x, and let the adversarial image be x′ = x(K) with K being the total number of iterations. For the Half-moon dataset, we apply 20 iterations PGD attack (PGD20) with step size 0.01. The perturbation is constrained to be 0.2 under l∞ norm. For the CIFAR-10/CIFAR-100 dataset, when under `∞ norm attack, we apply 20 iterations PGD attack with step size 1/255 and when under `2 norm attack and we apply 50 iterations PGD attack (PGD50) with step size 0.1.
Performance evaluations. For the Half-moon dataset, we consider three indicators to evaluate the performance of the model: natural testing accuracy (Atestnat ), robust testing accuracy (Atestrob ), and the gap between training accuracy and testing accuracy(Atrainnat −Atestnat ). The natural accuracy and robust accuracy are measured on clean and adversarial images, respectively. The robust accuracy is used to evaluate the robustness of the models, and the difference between training and testing accuracy is used to evaluate the overfitting risk of the models. For the CIFAR-10/CIFAR-100 datasets, we consider Atestnat and Atestrob .

heading: 4.2 EXPERIMENTS ON SYNTHETIC DATASET
text: In this subsection, we numerically verify that the efficacy of our training method in improving the robustness and generalization ability of ResNet on the Half-moon dataset.
We first varied hyperparameters λ and σ2 and present the final results in Figure 2. The three indicators (natural accuracy, robust accuracy, and difference between training accuracy and testing accuracy) for ResNet with natural training were 88.6%, 62.3%, and 10.4%, respectively. Comparing natural trained ResNet, the robustness and generalization ability of ResNet trained with our method were improved. On the other hand, when λ and σ2 increased, the natural and robust accuracy increased, while the gap between training and test accuracies decreased. Thus we can conclude that increasing λ and σ2 are helpful for improving the performance of the models.
We plot the decision boundary of ResNet trained with different λ, σ2 and natural trained ResNet in Figure 1. From Figure 1, we can see that the decision boundary of natural training is irregular, while that of our models is smoother. These experimental results are consistent with our theory.

heading: 4.3 EXPERIMENT ON CIFAR-10/CIFAR-100 DATASET
text: We further test the performance of ResNet20/ResNet56 trained with our method on the CIFAR10/CIFAR100 dataset. During training, we apply data augmentation including random crops and flips. We run 200 epochs. The batch size is 128 and the initial learning rate is 0.1, which decays by a factor of 10 at the 100th, and 150th epochs.
There are two hyperparameters λ, σ2 in our algorithm that need to be determined by a large number of experiments. In order to demonstrate the influence of the two hyperparameters, we train ResNet by using our algorithm with different λ, σ2, then list the result in Table 1. On one hands, from Table 1, we can see that the parameter σ2 is fixed, the natural accuracy decreases and robust accuracy increases as the increasing of parameters λ and the converse to be true. On the other hand, comparing to the result of ResNet20 and ResNet56, we can see that the performance of our proposed method will be better in deeper networks.
In order to quantitatively measure the performance of our proposed method, we compare the proposed method with other defensive method without adversarial training. The result of comparison is shown in Table 2. For our method, we choose σ2 = 2 and λ = 0.5 to train ResNet110. For the Gaussian noise injection mentioned in subsection 2.2, the standard error of the Gaussian noise equals to 0.5. The experiment result shows that our method is more resistant to adversarial attack than other methods.

heading: 5 CONCLUDING REMARKS
text: In this paper, we derived a general form of PDE that can correspond to a DNN model. Motivated by PDE theory, we proposed a training method for ResNet. We have theoretically shown that using our training methods, we obtain a model with better robustness and lower overfitting risk. Based on these theoretical results, we developed an algorithm to train DNN models, and we verified them through experiments. There are numerous avenues for future work: 1) Scale our training method to networks that are large and expressive enough to solve problems like ImageNet. 2) Study other more complex equations for DNN models and solve more challenging real tasks.

heading: A PROOF OF THEOREM 1
text: Proof. (Alvarez et al., 1993) The proof of Theorem 1 consists of several steps. First we set
δt,s(f) = Tt(f)− Ts(f)
t− s , δt(f) = δt,0(f).
Using [Temporal Regularity], for all h, t ∈ [0, T ] and all f, g ∈ Q, we can get ‖δt(f + hg)− δt(f)‖L∞ ≤ Ch. (18)
Step 1. We begin the proof of Theorem 1 by proving the following conclusion
(A1) the functions δt(f) are Lipschitz continuous uniformly for t ∈ (0, 1](and uniformly for f
belonging to a set Q). Let z be in Rd and ‖z‖2 = 1, we wish to bound ‖τhz[δt(f)]− δt(f)‖L∞ for h ∈ [0, 1]. According to [Spatial Regularity], for f ∈ Q, we have
‖τhz[δt(f)]− δt(f)‖L∞ ≤ ‖τhz[δt(f)]− δt[f(x+ hz)]‖L∞ + ‖δt[f(x+ hz)]− δt[f(x)]‖L∞ ≤ Ch.
Then the conclusion of (A1) is proved. Step 2. From the preceding steps, we can already deduce a compactness property for δt(f) as t→ 0. Now we need more, since we want the whole sequence to converge to the same limit. So we need a Cauchy estimate. Our next step is to prove the following conclusion.
(A2) For all t, s in [0, 12 ] we have ‖δt+s,t(f)− δs(f)‖L∞ ≤ m(t), where m(t) is some continuous,
nonnegative, nondecreasing function on [0, 12 ] such that m(0) = 0, and m(t) depends only on bounds
of derivatives of f . First, we can rewrite δs(f) as δs(f) ∗ Kε − [δs(f) ∗ Kε − δs(f)], where K ≥ 0, ∫ Rd Kdy = 1, K ∈ C∞0 (Rd) and Kε = ε−dK(·/ε) and ∗ denote the convolution. Using the conclusion of (A1), we can obtain for all ε > 0, there exist a positive constant C1 such that
‖δs(f) ∗Kε − δs(f)‖∞ ≤ C1ε. (19) Then because of [Markov Property], [Temporal Regularity] and [Comparison Principle], we have
‖Tt+s(f)− Tt(f + sδs ∗Kε) s ‖L∞
≤ ‖Tt+s,s ◦ Ts(f)− Tt ◦ Ts(f) s ‖L∞ + ‖ Tt ◦ Ts(f)− Tt(f + sδs ∗Kε) s ‖L∞ ≤ Ct+ ‖Ts(f)− f − sδs ∗Kε s ‖L∞
≤ Ct+ C2ε (20) On the other hand, since δs(f) ∗Kε ∈ C∞b , [Temporal Regularity] implies that
‖Tt(f + sδs(f) ∗Kε)− Tt(f)− sδs(f) ∗Kε‖L∞ ≤ Cεst (21) for some positive constant Cε depending only on ε. Combining equation 19, equation 20 and equation 21, we finally deduce that
‖δs+t,t(f)− δs(f)‖∞ ≤ 2C1ε+ Cεt. The (A2) holds by setting m(t) = infε∈(0,1](2C1ε+ Cεt). Step 3. Now we give a Cauchy estimate in 0 < h ≤ t ≤ 12 for δs(f)
‖δt(f)− δh(f)‖L∞ ≤ 2 C0r t +m(t) where r = t− [ t h ]h. (22)
Observe that δt(f) = Nh
t δNh(f) +
r t δNh+r,Nh(f).
Using (A2) with s = r, t = Nh and let N = [ th ], we have
‖δNh+r(f)− Nh
Nh+ r δNh(f)−
r
Nh+ r δr(f)‖L∞ ≤
r
Nh+ r m(Nh) (23)
Notice that
δNh(f) = (N − 1)h Nh δNh(f) + h Nh δNh,(N−1)h(f).
Using (A2) with s = h, t = (N − 1)h, we have
‖δNh(f)− N − 1 N δ(N−1)h(f)− 1 N δh(f)‖L∞ ≤ 1 N m((N − 1)h) (24)
Combining this inequality with equation 24, we obtain
‖δt(f)− (N − 1) h
t δ(N−1)h(f)−
h t δh(f)− r t δr(f)‖L∞
≤ r t m(Nh) + h t m((N − 1)h) (25)
Reiterating the argument which lends to equation 24, we obtain after (N − 1) more steps that
‖δt(f)− Nh
t δh(f)−
r t δr(f)‖∞ ≤ r t m(Nh) + h t N−1∑ j=1 m(jh) (26)
Since m(t) is nondecreasing and equation 25, we deduce
‖δt(f)− δh(f)‖L∞ ≤ ‖δt(f)− Nh
t δh(f)−
r t δr(f)‖L∞ + r t ‖δh(f)− δr(f)‖L∞
≤ r t m(Nh) + (N − 1)h t m(t) + r t (‖δh(f)‖L∞ + ‖δr(f)‖L∞)
≤ m(t) + 2C0r t
(27)
Because of equation 27, we can pick hn going to 0 and δhn(f) converges uniformly on compact sets to a bounded Lipschitz function on Rd which we denote by A[f ]
lim n→∞ ‖δt(f)− δhn(f)‖L∞ ≤ lim n→∞
m(t) + 2 C0r
t ⇒ ‖δt(f)−A[f ]‖L∞ ≤ m(t).
So δt(f) converges uniformly to A[f ] when t goes to 0. Similarly, there exist an operator At such that δs,t(f) converges uniformly to At[f ] when s goes to t. Step 4. The goal of this step is to show the expression of A[f ]. Let f, g ∈ C∞b and satisfy f(0) = g(0) = 0 (if not equal to 0, we replace f(x), g(x) by f(x) − f(0), g(x) − g(0)), Df(0) = Dg(0) = p ∈ Rd, D2f(0) = D2g(0) = A ∈ Rd×d. Then we set fε = f + ε‖x‖22. Using Taylor formula, there exist a positive constant c such that for ‖x‖2 ≤ cε we have fε ≥ g. Let wε = w(x/ε) where w ∈ C∞b (Rd) and 0 ≤ w ≤ 1 on Rd,{
w = 0 ‖x‖2 ≤ c/2 w = 1 ‖x‖2 ≥ c
(28)
Then we introduce fε0 = wεf ε+(1−wε)g and it is obviously that fε0 ≥ g. Because of [Comparison Principle], Tt(fε0 ) ≥ Tt(g). Due to
fε0 (0) = f ε(0) = f(0) = g(0) = 0 (29)
we can also get
A(fε0 )(0) ≥ A(g)(0), A(fε)(0) ≥ A(g)(0) (30)
Because there exists a neighborhood of 0 that fε0 = f ε, we have Dαfε0 (0) = D αfε(0) for ∀ |α| ≥ 0. In view of [Locality] we have A(fε0 )(0) = A(f
ε)(0). And considering the continuity of A, we can deduce A(fε0 ) converges to A(f) in L
∞ when ε goes to 0. This means A(f)(0) ≥ A(g)(0). By similar method, we can get A(f)(0) ≤ A(g)(0) which means A(f)(0) = A(g)(0). Hence, for any x0 ∈ Rd, if we replace 0 by x0, the similar result can be obtained. So the value of A(f)(x) only depends on x, f,Df,D2f . Observe that for any constant C, A(f +C) = A(f), we can get A(f)(x) only depends on x, Df,D2f . At last, we finally get
A(f) = F (Df,D2f,x). (31)
Here F is a continuous function. In the same way, we can get
At(f) = F (Df,D 2f,x, t). (32)
According to [Ensemble Additivity], since
F (Df,D2f,x, t) = lim t→0 (Ttf − f)/t,
F therefore satisfies
F (rDf + sDg, rD2f + sD2g,x, t) = rF (Df,D2f,x, t) + sF (Dg,D2g,x, t)
for any real numbers r and s and any functions f and g and at any point (x, t). Since the values of Df,Dg,D2f,D2g are arbitrary and can be independently taken to be 0, we obtain for any vectors v1, v2 and symmetric matrices A1, A2 and any fixed point (x0, t0) that
F (rv1 + sv2, rA1 + sA2,x0, t0) = rF (v1,A1,x0, t0) + sF (v2,A2,x0, t0)
F (v1,A1,x0, t0) = F (v1, 0,x0, t0) + F (0,A1,x0, t0).
Thus we can finally there exist Lipschitz continuous functions v : Rd × [0, T ]→ Rd and Lipschitz continuous positive definite functions σ : Rd × [0, T ] → Rd×d such that u(x, t) = Tt(f) is the solution of the equation{
∂u(x,t) ∂t = v(x, t) · ∇u(x, t) + ∑ i,j σi,j(x, t) ∂u xi,xj
(x, t),x ∈ Rd, t ∈ [0, 1] u(x, 0) = f(x),
where σi,j(x, t) is the i, j-th element of matrix function σ(x, t).
The existence and uniqueness of equation 5 are guaranteed by the following theorem Theorem 4. (Kolmogorov’s backward equation)(Theorem 8.1.1 in (Oksendal, 2013)) Assume It̂o process x(t) satisfies SDE
dx(t) = v(x(t), t)dt+ σ(x(t))dB(t), x(0) = x,
where B(t) is d-m Brownian motion and v : Rd → Rd, σ : Rd → Rd×d are all Lipschitz continuous function. Let fw(x) ∈ C20 (Rd). Then u(x, t) = E[fw(x(t)|x(0) = x] is the unique solution of equation 5.

heading: B PROOF OF THEOREM 2
text: To begin the proof of Theorem 2, we first provide the maximum principle for operator L which has the following form
Lu = ∂u
∂t − σ2∆u+ v(x, t)∇u,
where the coefficients v(x, t) are continuous. Theorem 5. (Evans, 1998) (Maximum principle) Assume U to be an open bounded subset of Rn and as before set UT = U × (0, T ] for some fixed time T > 0. Let u ∈ C2(UT ) ∩ C(ŪT ) and ŪT is the closure of UT . Denote ŪT − UT by ΓT . If
Lu ≤ 0 in UT , then
max ŪT u = max ΓT u.
After the introduction of maximum principle, we next illustrate the proof of the theorem 2.
Proof. Let w(x, t) = (µu2(x, t) + ‖∇u(x, t)‖2)e−2λt, where µ and λ are constants which will be defined later. Note that u2(x, t) satisfies
∂(u2)
∂t + v(x, t)∇(u2) = σ2∆(u2)− 2σ2‖∇u‖2,
and ‖∇u‖2 satisfies
∂‖∇u‖2
∂t + v(x, t)∇‖∇u‖2 = −2∇u∇v(x, t)∇u+ σ2∆‖∇u‖2 − 2σ2(∆u)2,
therefore,
∂w ∂t + v(x, t)∇w − σ2∆w
=e−2λt[−2λ(µu2 + ‖∇u‖2)− 2µσ2‖∇u‖2 − 2∇u∇v(x, t)∇u− 2σ2(∆u)2].
Next, let γ(x, t) = min‖ξ‖=1 ξT∇v(x, t)ξ and γ = minx,t γ(x, t), and let Lw := ∂w∂t + v(x, t)∇w − σ2∆w then we have
Lw ≤ −2e−2λt[λµu2 + (λ+ µσ2 − γ)‖∇u‖2].
If we choose λ and µ large enough, such that λ+ µσ − γ ≥ 0, then
Lw ≤ 0.
From the maximum principle (Theorem 5), we know maxx w(x, 1) ≤ maxx w(x, 0). Let µ = 1 and λ = γ − σ2, we can easily get
‖u(x, T )‖C1 ≤ e−σ 2T eγT ‖fw(x)‖C1 , (33)
According to the Taylor formula and Hölder inequality, we can get that for any ‖δ‖p = δ, we have
|fw(x)− fw(x+ δ)| ≤ ‖∇fw(x)‖qδ ≤ d1/q‖fw(x)‖C1δ, (34)
and
|u(x, T )− u(x+ δ, T )| ≤ ‖∇u(x, T )‖qδ ≤ d1/q‖u(x, T )‖C1δ, (35)
According to equation 34 and equation 35, we can get
ucA(x+ δ, T ) ≥ pA − d1/q‖u(x, T )‖C1δ max i 6=cA ui(x+ δ, T ) ≤ pB + d1/q‖u(x, T )‖C1δ.
Thus ucA(x+ δ, T ) ≥ maxi6=cA ui(x+ δ, T ) for all
δ ≤ pA − pB 2d1/q‖u(x, T )‖C1
Then combining equation 33, equation 34 and equation 35 we can easily prove the theorem.

heading: C PROOF OF THEOREM 3
text: The goal of the learning problem is to find h ∈ H such that population risk E(h) := E(x,y)[`(h(x), y)] is minimized. We consider the supervised learning setting where one has access to n i.i.d. training examples SN . A learning algorithm maps the N training examples to a hypothesis h ∈ H. In this section, we are interested in the gap between the empirical risk EN (h) := 1 N ∑N i=1 `(h(xi), yi) and the population risk E(h), known as the generalization error. We start with the following theorem which connects the population and empirical risks via Rademacher complexity.
Lemma 1. (Theorem 3.5 in (Mohri et al., 2018)) Let SN = {xi, yi}Ni=1 be samples chosen i.i.d. according to the distribution p(x). If the loss function ` is bounded byB > 0. Then for any δ ∈ (0, 1), with probability at least 1− δ, the following holds for all h ∈ H,
E (h) ≤ ESN (h) + 2RSN (`H) + 3B √ (log(2/δ)/(2N).

heading: In addition, according to the Ledoux-Talagrand contraction inequality(Maurer, 2016) and assume
text: loss function is L-Lipschitz, we have
RSN (`H) ≤ LRSN (H) .
Next, we begin to prove Theorem 3. Theorem 6. (Theorem 2.3 in (Mendelson, 2003)) Let H be a class of functions from D to [−1, 1] and set p to be a probability measure on D. Let (xi)ni=1 be independent random variables distributed according to µ. Let (Y, d) be a metric space and set F ⊂ Y . For every > 0 and any n ≥ 8/ 2,
P ( sup f∈H | 1 n n∑ i=1 f(xi)− ∫ D f(x)µ(x)dx| > ) ≤ 8Eµ[Num( /8,F , L1(µn))] exp(−n 2/128) (36)
In above theorem, Num( ,H, Lp(µn) denotes the covering numbers of H at scale with respect to the Lp(µn) norm. µn is the empirical measure supported on one sample of (xi)ni=1. For every > 0, denote by Num( ,H, d) the minimal number of open balls (with respect to the metric d) needed to coverH. That is, the minimal cardinality of the set {y1, · · · , ym} ⊂ Y with the property that every f ∈ H has some yi such that d(f, yi) < . The set {y1, · · · , ym} is called an -cover of H . The logarithm of the covering numbers is called the entropy of the set. For every sample {x1, · · · ,xn} let µn be the empirical measure supported on that sample. For 1 ≤ p < ∞ and a function f , ‖f‖Lp(µn) = ( 1 n ∑n i=1 |f(xi)|p )1/p and ‖f‖∞ = max1≤i≤n |f(xi)|.
Notice that
L1(µn) ≤ L∞(µn) ≤ L∞.
We get one immediate corollary of Theorem 6. Corollary 1. LetH be a class of functions fromD to [−1, 1] and set p(x) to be a probability measure on D. Let {xi}Ni=1 be independent random variables distributed according to p(x). For every > 0 and any n ≥ 8/ 2,
P ( sup f∈H | 1 N n∑ i=1 f(xi)− ∫ D f(x)p(x)dx| > ) ≤ 8Num( /8,H, L∞) exp(−N 2/128) (37)
where Num( ,H, L∞) be the covering numbers ofH at scale with respect to the L∞ norm
Then, we get an upper bound of supf∈H | 1N ∑N i=1 f(xi)− ∫ D f(x)p(x)dx|.
Corollary 2. LetH be a class of functions from D to [−1, 1]. Let {xi}Ni=1 be independent random variables distributed according to p, where p is the probability distribution whose C1 norm is bounded. Then with probability at least 1− δ,
sup f∈H |p(f)− pN (f)| ≤ √√√√128 N ( ln Num( √ 2 N ,H, L∞) + ln 8 δ ) ,
where
p(f) = ∫ D f(x)p(x)dx, pN (f) = 1 N N∑ i=1 f(xi). (38)
Proof. Using Corollary 1, with probability at least 1− δ,
sup f∈H |p(f)− pN (f)| ≤ δ,
where δ is determined by
δ =
√ 128
N
( ln Num( δ/8,H, L∞) + ln 8
δ
) .
Obviously,
δ ≥ √ 128
N = 8
√ 2
N which gives that
Num( δ/8,H, L∞) ≤ Num( √ 2
N ,H, L∞)
Then, we have
δ ≤ √√√√128 N ( ln Num( √ 2 N ,H, L∞) + ln 8 δ ) which proves the corollary.
If the entropy bound of H is known, the upper bound of supf∈H |p(f) − pN (f)| follows from Corollary 2. Now, the key point left becomes bounding the entropy of some given function classH. Let us start from the function class F . To apply above corollary, we need to normalize F to make it lie in [−1, 1]. Here we also use F to denote the normalized function class and absorb the bound of F into the generic constant C.
According to Taylor formula, we have for any x, z ∈ D |f(x)− f(z)| ≤
√ d sup f∈F ‖f‖C1‖x− z‖2.
Where d is the dimension of x. This gives an easy bound of Num( ,F , L∞), Num( ,F , L∞) ≤ C ( supf∈F ‖f‖C1 )d (39)
Using Corollary 2, with probability at least 1− 1/(2N), we have
sup f∈F |p(f)− pN (f)| ≤ C√ N
( lnN + ln(sup
f∈F ‖f‖C1) + 1
)1/2 (40)
and
sup f∈Gσ |p(f)− pN (f)| ≤ C√ N
( lnN + ln( sup
f∈Gσ ‖f‖C1) + 1 )1/2 According to equation 33, we can get the following theorem. Theorem 7. With probability at least 1− 1/(2N),
sup f∈F |p(f)− pN (f)| ≤ C√ N (lnN + ln c+ 1) 1/2 ,
sup f∈Gσ |p(f)− pN (f)| ≤ C√ N
( lnN + ln c− σ2 + 1 )1/2 p(f) =
∫ D f(x)p(x)dx, pN (f) = 1 N ∑ xi∈SN f(xi),

heading: F and Gσ is a function class defined as
text: F := { f(x)|‖f‖C1 ≤ c } ,
and Gσ := { g(x) = u(x, 1)|∂u(x, t)
∂t = v(x, t) · ∇u(x, t) + 1/2σ2∆u(x, t), u(x, 0) = f(x) where f ∈ F
} .
For the purpose of proving Theorem 3, we give a inequality about Rademacher complexity as following. Theorem 8. (Theorem 2.23 in (Mendelson, 2003)) Let p(x) be a probability distribution function and setH to be a class of functions on D. Denote
Z = sup f∈H | 1 N N∑ i=1 f(xi)− Ex∼p(x)(f)|,
where {xi}Ni=1 are independent random variables distributed according to p(x). Then we have 1
2 Ex∼p(x)Z ≤ RSN (H) ≤ 2Ex∼p(x)Z + 1√ N | sup f∈H Ex∼p(x)f |.
Then combining Theorem 7, Theorem 8, we have
RSN (Gσ) ≤ 2Ex∼p(x)Z + 1√ n | sup f∈Gσ Ex∼p(x)f |.
≤ C√ N
( lnN + ln c− σ2 + 1 )1/2 +
1√ N | sup f∈Gσ Ex∼p(x)f |
≤ C√ N
( lnN + ln c− σ2 + 1 )1/2 +
1√ N sup f∈Gσ ‖f‖C1
≤ C√ N [(−σ2 + ln c+ lnN)1/2 + exp(−σ2)c].
Similarly, we can estimate the upper bound of RSN (F).

heading: D MORE EXPERIMENTAL DETAILS ON HALF-MOON DATASET
text: In this subsection, we plot the decision boundary of ResNet trained with different hyperparameters in Figure 3. As drawn in Figure 1, we can see as the hyperparameters λ and σ2 increasing, the decision boundary becomes more smooth, which implies the overfitting risk reduce. These experimental results are also consistent with the conclusion of Theorem 3.

heading: E ARCHITECTURES OF THE USED DNNS
text: Figure 4 shows the architectures of ResNets used in this paper. We also plot the inputs for the ResNet used on CIFAR-10 dataset in Figure 5 and every element of the matrix tt equals to t.

Review summaries:
This paper studied a ResNet-like DNN model that can be expressed as a discretization of PDEs. First, this paper showed that, under some assumptions, any adjust operator is a solution of a second-order convection-diffusion PDE (Theorem 1). ResNets and ResNets with Gaussian noise injection are special cases of this theorem. Next, for the specific PDE (Eq. (7)), this paper derived generalization guarantees regarding the Rademacher complexity. This paper also derived robustness guarantees in terms of input perturbations. Finally, this paper analyzed the predictive performance of ResNet trained by the proposed method with various model hyperparameters on clean and adversarial datasets.

This paper considers DNNs as transformations of the input data which can be seen as discrete approximations to the solutions of PDEs. It reformulates a point of view which has been trending for a few years now. More precisely, this paper represents the network as an operator, which is actually solution to a PDE, which acts on a "base classifier". Using facts from PDE theory, they then propose improvements for standard neural architectures.

This paper connects the architecture of a resnet and a gaussian injected resnet to that of a transport equation and a diffusion equation respectively. They show that under this framework they provide robustness guarantees of their PDE framework scale with the $\sigma$ parameter of the diffusion equation. Furthermore, they also show that larger the $\sigma$ better the generalization gap for the resulting neural network.   They introduce a learning algorithm, to attain the network from their PDE framework starting from an initial NN (also a resnet).  They also empirically verify that their method improves robustness, and that it performs better than Gaussian Noise injection.

The authors cast DNN classifiers as the push-forward of a base classifier under a flow map at some fixed final time. Under some reasonable assumptions on the flow, they show that, given any base classifier, the flow map can be obtained as the solution to a convection-diffusion equation. They show that ResNets and Gaussian noise injection can be viewed as special cases of their model and give a robustness guarantee for any  classifier defined as a solution to their PDE. Experiments on the 2-d half moon data set as well CIFAR 10 and 100 show better robustness  of their model to adversarial attacks  when compared to standard ResNet(s).