# Simulation of the Kinetic Energy Spectrum of Primary Electron
Generated by a 6 keV X-Ray Source in a Gas Volume

Lucas Maia Rios<EMAIL_ADDRESS>Mauro Rogerio Cosentino
<EMAIL_ADDRESS>Federal University of ABC

(May 18, 2018)

###### Abstract

With the intention of creating a complete simulation of the Thick GEM
detector, this work presents a partial simulation of the interaction of 6 keV
X-Ray radiation with a gas volume. For the gas, it was used pure Argon, Xenon
and CO2 as well as a mixture of each noble gas with CO2 in a 70%-30%
proportion. For this range of energy, the dominant physical effect is the
photoelectric, therefore knowing the binding energies of the atoms we predict
the kinetic energy of the primary electron. A simulated energy spectrum of the
primary electron kinetic energy was obtained, which contains peaks
corresponding to the inner shell’s ionization of the atoms. The spectrum shows
very little Compton effect, transferring only small energies (in the order of
eVs). We concluded that energy spectrum for the mixture is a linear
combination of the spectrum of its individual components.

.

Keywords

electron energy spectrum, X-Ray beam, simulation, Geant4.

###### pacs:

Valid PACS appear here

##  I Introduction

###  I.1 Motivation

Gaseous detectors are present in the fields of nuclear and particle physics
since the beginning, but it was only after the invention of the Multi-Wire
Proportional Chamber (MWPC) that they became present in almost every
experiment. There are applications for them in medical science, engineering
and even in every day use technology. The application we are interested in to
further study is the detection of charged particles created in a particle
accelerator collision. The ALICE experiment (A Large Ion Collider Experiment)
is a multi-detector experiment at the LHC at CERN that is designed to detect
particles generated in heavy-ion collisions. The goal of the experiment is to
study the quark-gluon plasma, which is studied indirectly through the decay
modes of the plasma.

To detect these particles, it is currently in use the MWPC technology, within
a Time Projection Chamber (TPC) tpc_alice , aimed to reconstruct the track of
charged particles. A MWPC is a chamber with anode wires between cathode planes
arranged in such a way to create an electric field that captures the electrons
from ionizing radiation. The electron travels through the drift area, to the
multiplication area where the electric field intensifies and triggers an
avalanche of electrons that induce a signal in the wires.

One problem of using MWPC though, is that it creates a significant ion back
flow (IBF) as a result of the multiplication. As the ions are heavier than the
electrons, they are much slower. Therefore, they accumulate in the drift area.
This concentration of ions distorts the electric field, thus distorting the
tracking reconstruction. To prevent this a gating grid is set in place, in
order to block the passage of the ions, but it also prevents electrons from
the drift area to reach the anode wires. Therefore, the gating grids must be
opened for a period at each event, short enough to prevent the IBF but long
enough to allow the electrons to go through.

In 2019 the LHC will undergo a shutdown for upgrades on the luminosity of the
beam (to increase the number of collisions per second) as well as other
improvements in the experiments. As for the ALICE experiment upgrade
alice_upgrade , it is expected to tackle this IBF issue by switching the MWPC
to GEM (Gas Electron Multiplier) Sauli  in the TPC, since the higher collision
rate makes it impossible for the gating grid to work correctly.

Figure 1: Typical GEM electrode, 50 $\mu m$ thick, 140 $\mu m$ pitch and 70
$\mu m$ diameter. Sauli

GEMs are a pair of foils with a dieletric layer in between, with equally
spaced holes as shown in figure 1. There is an electric potential difference
between these two foils so it creates an intense electric field in the holes.
The charged particles will tend to follow the eletric field lines as shown in
figure 2. The primary electron from an ionization travels through the
equipotential lines in the drift area until it reaches the multiplication area
inside the hole. On the other side of the GEM comes out a greater number of
electrons. Whereas ions follow the opposite path through the equipotential
lines to either return to the drift area or to be deposited in the GEM foil,
reducing the IBF. Also, one can arrange several GEMs in series so the ions
created are captured by the previous GEM. Only a small fraction of the ions go
to the drift area.

Figure 2: Electric field in the region of the holes. Sauli

GEMs are considerably expensive and sensitive (a short-circuit could damage
them). For that reason, we have been developing a Thick GEM prototype.
Conceptually there is no difference between the two. However, Thick GEMs are
more robust. Their holes are significantly larger and the dieletric layer is
thicker (this is why they are called Thick GEM). The disadvantage is the
reduced position resolution.

Thick GEMs, for their low price and better efficiency than MWPC, could be used
for several applications. Here at Federal University of ABC, there are
laboratories that work with X-Ray crystallography that still use MWPC that
could benefit from the use of the Thick GEM.

This work is an ongoing project to develop a complete simulation of an
incident monochromatic X-Ray beam, ionization the gas and then detection of
the electrons by the Thick GEM. In order to have the complete simulation,
three softwares will be needed. Geant4 geant4  is used to simulate the
interaction of the radiation with matter. Then Garfield++ garfieldpp  will be
used to simulate the electron transport, multiplication and signal generation.
At last, ROOT root  is used to analyze the data.

This paper contains the simulation of the radiation interaction with the gas
and a study of the electron spectrum created by it.

###  I.2 Physics

Ionizing radiation can interact with matter through many ways. Some of which
result in the creation of charged particles. The particle of interest in this
study that will be interacting is the photon. A source of X-ray radiation
emits photons that, when compared to charged particles, can better penetrate a
material because of its neutral nature. Nevertheless, it may interact with the
electrons of the material along its path. The physical processes that a photon
can undergo are the photoelectric effect, Compton scattering, pair production,
Thomson and Rayleigh scattering. The latter two are non-ionizing elastic
scattering, thus they will not be discussed. Also, for the pair production to
occur it is necessary that the photon has an energy greater than a certain
threshold (at least the rest mass energy of the produced particles). The rest
energy of an electron and a positron together is 1022 keV. At the energy we
are working there will be no pair production, so we will put it aside as well.
That leaves us with photoelectric effect and Compton scattering, both transfer
energy to the medium that manifests as an electron gaining kinetic energy. In
the first case the photon is completely absorbed by an atom, which ejects an
electron as a result. In the latter the photon only transfers a fraction of
its energy, allowing it to further interact via either one of the processes or
to escape the volume of interest, without further interactions.

The electron bounded to an atom is trapped within a potential barrier. This
barrier is interpreted as the binding energy. In an interaction with ionizing
radiation it is more likely that the innermost electron in the atom is freed.
leo  With this in mind, we can write the kinetic energy of an electron ejected
via photoelectric effect as:

$K_{e}=E_{\gamma}-\phi,$ (1)

where $E_{\gamma}$ is the incoming photon energy and $\phi$ the binding
energy.

Given the atoms in the gaseous medium, and knowing its binding energies, we
can foretell the kinetic energy of the ejected electron. Also note that it
will be a discrete spectrum, since the binding energies are discrete and our
X-ray ray is monochromatic. The binding energy for the innermost electron in
Ar-40, for example, is 3205.9 eV xray2 . Therefore, using equation 1 we expect
our electron to have a kinetic energy of 2794.1 eV. All of the binding
energies for Argon and Xenon can be seen in table 1, as well as the expected
kinetic energy for the electron.

| Ar | Xe | Kinetic E. ($eV$)
---|---|---|---
$K$ 1s | 3205.9 | 34561 | 2794.1 / -
$L_{1}$ 2s | 326.3 | 5453 | 5673.7 / 547
$L_{2}$ 2p1/2 | 250.6 | 5107 | 5749.4 / 893
$L_{3}$ 2p3/2 | 248.4 | 4786 | 5751.6 / 1214
$M_{1}$ 3s | 29.3 | 1148.7 | 5970.7 / 4851.3
$M_{2}$ 3p1/2 | 15.9 | 1002.2 | 5984.1 / 4997.8
$M_{3}$ 3p3/2 | 15.7 | 940.6 | 5984.3 / 5059.4
$M_{4}$ 3d3/2 |  | 689.0 | 5311.0
$M_{5}$ 3d5/2 |  | 676.4 | 5323.6
$N_{1}$ 4s |  | 213.2 | 5786.8
$N_{2}$ 4p1/2 |  | 146.7 | 5853.3
$N_{3}$ 4p3/2 |  | 145.5 | 5854.5
$N_{4}$ 4d3/2 |  | 69.5 | 5930.5
$N_{5}$ 4d5/2 |  | 67.5 | 5932.5
$O_{1}$ 5s |  | 23.3 | 5976.7
$O_{2}$ 5p1/2 |  | 13.4 | 5986.6
$O_{3}$ 5p3/2 |  | 12.1 | 5987.9
Table 1: Binding energy for Ar and Xe and their respective electron kinetic
energy when ejected by 6 $keV$ photon. Data collected from X-Ray Data Booklet
which used xray1  and xray2 .

The Compton scattering, on the other hand, is a collision between the photon
and a free electron. The amount of energy transferred is random, we can not
know it exactly beforehand. However, we know that the loss of energy by the
photon is completely gained by the electron as kinetic energy.

Differently from the photoelectric effect, the electrons will gain kinetic
energy within a continuous range of energy. Also, free electrons in gases are
not as abundant as in metal. According to pdg  the Compton cross section for
this range of energy is smaller than the photoelectric by 2 orders of
magnitude. So it is expect that the latter should be the dominant process for
6 keV photons.

Our interest in this work is to simulate the X-ray interaction with a gas
volume. We chose an energy of 6 keV because of the common use of Fe-55 as a
source of X-Ray radiation, which has energy emissions about 6 keV. The reason
for simulating X-Ray is due to the prospect of implementing the Thick GEM in
the X-Ray crystallography laboratory. We expect to get as an output the
electron energy spectra for different gas mixtures. From that we can compare
with the theoretical predictions to test the reliability of our simulation.

##  II Simulation

We used the Geant4 toolkit geant4 . It is an object-oriented C++ library of
classes for simulation of radiation interaction with matter. With it we can
define the geometry of our detector and its materials. For the primary event,
the first particle created that is, the photon is defined. Its momentum
direction is centered in the origin of the XY surface and it has a normal
distribution along the Z axis of standard deviation of $\sigma=1/12\pi$ rad.

Considering the discussion presented in the previous section, the processes
added to the simulation were the photoelectric effect and the Compton
scattering.

The initial purpose of Geant4 was simulating high energy physics. Because of
that there has been for some time a deficit in quality for low energy
interaction. In recent years, however, new models have been implemented such
as Livermore and Penelope. Both have a low limit applicability of 250 eV. For
the high limits they are 100 GeV and 1 GeV, respectively.

The primary event, i.e. the first particle (photons from our X-Ray source) of
the simulation, is generated in the origin, immerse in air. In 1 mm of
distance of the origin in the Z axis there is a volume filled with a gas
mixture (Argon-CO2 or Xenon-CO2, in 70%-30% or its pure constituents). This
volume is a cube of 5 cm of side. All of the physics of interest is in this
volume. The X-ray ray will interact with the particles of such gas and deposit
energy in it. It will also ionize it, creating electrons and ions.

In the simulation a number of events is set (each event is the creation of the
primary particle until it and all its secondaries are absorbed and there no
longer are particles to track). That number is currently ten thousand. A
particle is tracked step-by-step, in which we know all of its information
(position, momentum, energy, etc) before and after the current step. Every
time a particle penetrates the volume of our interest, information about the
hits (interactions) is collected. If there are electrons created by the X-ray,
their kinetic energy is saved as well.

Figure 3: Perspective view for Argon 70%. Photon trajectory in green.

Geant4 also has visualization tools as we can see in figure 3 and 4. They both
show the primary particle traveling into the volume and interacting with its
gas. In green is the photon. The primary electrons were to be seen in red,
however, because of their low energy they are quickly absorbed by the medium
and barely travel.

Figure 4: Side view for Xenon 70%. Photon trajectory in green.

In the side view (figure 4), we can visualize the collimation of the beam.
Another difference in the two figures is the composition of the gas. We can
see that because Xenon has a higher Z number the X-Ray travels shorter
distances than in Argon.

##  III Results

The output of the Geant4 simulation is a ROOT data file which contains the
information of position, momentum of the interactions as well as the energy
spectrum of the primary electron. We ran the simulation for the pure noble
gases (Ar and Xe), CO2 and their mixtures with 70% of noble gas and 30% of
CO2. In a histogram there are 200 channels per keV. We begin by looking at
spectrum created by Argon alone.

Figure 5: Primary electron energy spectrum for Argon 100%.

In figure 5 we can see that there are two visible peaks. One is in the middle
of the range and the other one is closer to 6 keV. If we take a closer look at
the peaks, figures 6 and 7, we can fit a gaussian in each peak using ROOT. As
a result their mean values and errors are 2.8228 $\pm$ 0.0002 keV and 5.7061
$\pm$ 0.0007 keV. For each peak from here on the fitting method will be used
to get their information.

Figure 6: Peak one for Argon 100%. Figure 7: Peak two for Argon 100%.

We expected that the main energy peaks would correspond to inner shells of the
atoms. From table 1 the energies for the K and L shell for Argon are rather
close to the peaks. That is, K shell has a 1.02% relative difference. For the
L shell however, it seems that there could be an overlap of the L shell peaks
in the spectrum because we see only one peak instead of three. Nonetheless,
the relative differences are respectively 0.57, 0.75 and 0.79%.

Figure 8: Primary electron energy spectrum for Xenon 100%.

As for the Xenon gas alone, we can see two sets of three peaks (figure 8). The
first set is between 0 and 2 keV and their mean values are 0.5830 $\pm$
0.0004, 0.8960 $\pm$ 0.0003 and 1.2252 $\pm$ 0.0002 keV. For the second set
there is some overlapping but we can still see the three different peaks with
mean values of 4.881 $\pm$ 0.002, 5.052 $\pm$ 0.001 and 5.3160 $\pm$ 0.0008
keV.

Figure 9: Lower energy peaks for Xenon 100%. Figure 10: Higher energy peaks
for Xenon 100%.

From table 1 we can infer that the first three peaks correspond to the three L
shell energies. Since the photon does not have enough energy to ionize the K
shell, L shell becomes the innermost ionizable one. The relative difference to
the peaks are respectively 6.58, 0.33 and 0.92%. The second set of peaks
correspond to the next innermost shell, the M shell. Their relative
differences are 0.61 ($M_{1}$), 0.15 ($M_{3}$) and 0.09 % ($M_{4}$) to the
theoretical prediction.

Figure 11: Primary electron energy spectrum for Argon 70%.

Now we can start looking at the mixtures. First for Argon 70% with CO2 30%, in
figure 11. We expect that the same peaks will appear, changing only their
amplitudes. The histogram of mix thus should be a linear combination of their
pure histograms. From figure 11 we can see that the peak for the K shell
remains. However, there is a distortion where the L shell would be due to the
presence of CO2, as we can see in figure 12. We have plotted in figure 12 both
the mixture and its individual components rescaled to show that when summed
they become the mixture histogram.

Figure 12: In black there is the simulated spectrum for the mixture 70% Argon
and 30% CO2. In red the pure Argon spectrum rescaled by 0.9 whereas in green
is the pure CO2 rescaled by 0.2.

Identically we can do the same for the Xenon 70% mixture in figure 13.

Figure 13: In black there is the simulated spectrum for the mixture 70% Xenon
and 30% CO2. In red the pure Xenon spectrum whereas in green is the pure CO2
spectrum rescaled by 0.1.

For Argon, CO2 and their mixture there are nothing else to explore in the
spectrum. For the Xenon, however, there are some counts that do not belong to
the Xenon binding energies because there are no shells between the L and M
that correspond to that energy (3-4 keV) as we can see in figure 8. One could
think those energy counts are due to Compton scattering, but figure 14 show us
that it only happens for energy in the order of eVs.

Figure 14: Compton scattering spectrum for Xenon 70%.

So according to the simulation those middle range counts should be
Photoeletric effect. It could be some mistake in the simulation. More
investigation needs to be done.

It should also be acknowledged that the measurement error depends of the
standard deviation we choose for the energy. Since the Fe-55 has different
energy emissions, we approximated the energy to a Gaussian with mean energy of
6 keV.

##  IV Conclusion

The simulation is able to create the energy peaks corresponding to the binding
energies of the gas filling the volume. The highest relative difference being
6.58%, however most of the rest are less than 1%. As the error associated to
the measurement in the simulation is mostly in the order of $10^{-2}$ or
$10^{-3}$%, the relative difference can not be ignored. Also, we can conclude
that the spectrum of the mixtures of gases is a linear combination of the
spectrum of its individual components.

Further study of the simulation is necessary to grasp the reason why there are
some scattered energy counts where there should not be. Experiments will be
made in the future with the Thick GEM, then it will be possible to confront
experimental data with the simulation.

## References

  * (1) The ALICE TPC, a large 3-dimensional tracking device with fast readout for ultra-high multiplicity events.  Nuclear Instruments and Methods In Physics Research, 622(1):316–367, 2010.
  * (2) S. Agostinelli et al.  GEANT4: A Simulation toolkit.  Nucl. Instrum. Meth., A506:250–303, 2003.
  * (3) J. A. BEARDEN and A. F. BURR.  Reevaluation of x-ray atomic energy levels.  Rev. Mod. Phys., 39:125–142, Jan 1967.
  * (4) Rene Brun and Fons Rademakers.  ROOT — an object oriented data analysis framework.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 389(1):81–86, 1997.
  * (5) M. Cardona and L. Ley.  Photoemission in solids: General principles.  Topics in applied physics. Springer-Verlag, 1978.
  * (6) W. R. Leo.  Techniques for Nuclear and Particle Physics Experiments: A How to Approach.  1987\.
  * (7) Christian Lippmann.  Upgrade of the ALICE Time Projection Chamber.  2014\.
  * (8) C. Patrignani et al.  Review of Particle Physics.  Chin. Phys., C40(10):100001, 2016.
  * (9) Fabio Sauli.  The gas electron multiplier (GEM): Operating principles and applications.  Nucl. Instrum. Meth., A805:2–24, 2016.
  * (10) H. Schindler and R. Veenhof.  Garfield++ - Simulation of Tracking Detectors.


# Stable Neural ODE with Lyapunov-Stable Equilibrium Points for Defending
Against Adversarial Attacks

Qiyu Kang
Continental-NTU Corporate Lab
Nanyang Technological University
50 Nanyang Avenue, 639798, Singapore
<EMAIL_ADDRESS>
&Yang Song11footnotemark: 1
School of Electrical and Electronic Engineering
Nanyang Technological University
50 Nanyang Avenue, 639798, Singapore
<EMAIL_ADDRESS>
&Qinxu Ding
School of Business
Singapore University of Social Sciences
463 Clementi Road, 599494, Singapore
<EMAIL_ADDRESS>
&Wee Peng Tay
School of Electrical and Electronic Engineering
Nanyang Technological University
50 Nanyang Avenue, 639798, Singapore
<EMAIL_ADDRESS>
First two authors contributed equally to this work.

###### Abstract

Deep neural networks (DNNs) are well-known to be vulnerable to adversarial
attacks, where malicious human-imperceptible perturbations are included in the
input to the deep network to fool it into making a wrong classification.
Recent studies have demonstrated that neural Ordinary Differential Equations
(ODEs) are intrinsically more robust against adversarial attacks compared to
vanilla DNNs. In this work, we propose a stable neural ODE with Lyapunov-
stable equilibrium points for defending against adversarial attacks (SODEF).
By ensuring that the equilibrium points of the ODE solution used as part of
SODEF is Lyapunov-stable, the ODE solution for an input with a small
perturbation converges to the same solution as the unperturbed input. We
provide theoretical results that give insights into the stability of SODEF as
well as the choice of regularizers to ensure its stability. Our analysis
suggests that our proposed regularizers force the extracted feature points to
be within a neighborhood of the Lyapunov-stable equilibrium points of the ODE.
SODEF is compatible with many defense methods and can be applied to any neural
network’s final regressor layer to enhance its stability against adversarial
attacks.

##  1 Introduction

Although deep learning has found successful applications in many tasks such as
image classification [1, 2], speech recognition [3], and natural language
processing [4], the vulnerability of deep learning to adversarial attacks
(e.g., see [5]) has limited its real-world applications due to performance and
safety concerns in critical applications. Inputs corrupted with human-
imperceptible perturbations can easily fool many vanilla deep neural networks
(DNNs) into mis-classifying them and thus significantly impact their
performance.

Recent studies [6, 7, 8] have applied neural Ordinary Differential Equations
(ODEs) [9] to defend against adversarial attacks. Some works like [6] have
revealed interesting intrinsic properties of ODEs that make them more stable
than conventional convolutional neural networks (CNNs). The paper [6] proposes
a time-invariant steady neural ODE (TisODE) using the property that the
integral curves from a ODE solution starting from different initial points
(inputs) do not intersect and always preserve uniqueness in the solution
function space. However, this does not guarantee that small perturbations of
the initial point lead to small perturbations of the integral curve output at
a later time $T$. The authors thus proposed a regularizer to limit the
evolution of the curves by forcing the integrand to be close to zero. However,
neither the non-intersecting property nor the steady-state constraint used in
TisODE can guarantee robustness against input perturbations since these
constraints do not ensure that the inputs are within a neighborhood of
Lyapunov-stable equilibrium points. An example is an ODE that serves as an
identity mapping is not robust to input perturbations but satisfies all the
constraints proposed in [6].

In this paper, our objective is to design a neural ODE such that the features
extracted are within a neighborhood of the Lyapunov-stable equilibrium points
of the ODE. We first develop a diversity promoting technique applied in the
final fully connected (FC) layer to improve the ODE’s stability and analyze
the reasons why. We then propose a stable neural ODE with Lyapunov-stable
equilibrium points to eliminate the effects of perturbations in the input.
From linear control theory [10], a linear time-invariant system
$\,\mathrm{d}\mathbf{z}(t)/\,\mathrm{d}t=\mathbf{A}\mathbf{z}(t)$, where
$\mathbf{A}$ is a constant matrix, is exponentially stable if all eigenvalues
of $\mathbf{A}$ have negative real parts. Specifically, we propose to force
the Jacobian matrix of the ODE used in the neural ODE to have eigenvalues with
negative real parts. Instead of directly imposing constraints on the
eigenvalues of the matrix, which lead to high computational complexity when
the Jacobian matrix is large, we instead add constraints to the matrix
elements to implicitly force the real parts of its eigenvalues to be negative.

Our main contributions are summarized as follows:

  1. 1.

Based on the concept of Lyapunov-stable equilibrium points, we propose a
simple yet effective technique to improve the robustness of neural ODE
networks by fixing the final FC layer to be a matrix whose rows have unit norm
and such that the maximum cosine similarity between any two rows is minimized.
Such a FC layer can be constructed off-line.

  2. 2.

We propose a stable neural ODE for deFending against adversarial attacks
(SODEF) to suppress the input perturbations. We derive an optimization
formulation for SODEF to force the extracted feature points to be within a
neighborhood of the Lyapunov-stable equilibrium points of the SODEF ODE. We
provide sufficient conditions for learning a robust feature representation
under SODEF.

  3. 3.

We test SODEF on several widely used datasets MNIST [11], CIFAR-10 and
CIFAR-100 [12] under well-known adversarial attacks. We demonstrate that SODEF
is robust against adversarial white-box attacks with improvement in
classification accuracy of adversarial examples under PGD attack [13] of up to
$44.02\%$, $52.54\%$ and $18.91\%$ percentage points compared to another
current state-of-the-art neural ODE network TisODE [6] on MNIST, CIFAR-10 and
CIFAR-100, respectively. Similar improvements in classification accuracy of
adversarial examples of up to $43.69\%$, $52.38\%$ and $18.99\%$ percentage
points compared to ODE net [9] are also obtained.

The rest of this paper is organized as follows. We provide essential
preliminaries on neural ODE and its stability analysis in Section 2. In
Section 3, we present SODEF model architecture and its training method. We
show how to maximize the distance between stable equilibrium points of neural
ODEs. We propose an optimization and present theoretical results on its
stability properties. We summarize experimental results in Section 4 and
conclude the paper in Section 5. The proofs for all lemmas and theorems
proposed in this paper are given in the supplementary material. We also refer
interested readers to the supplementary material for a more detailed account
of related works [14, 15, 16, 6] and some popular adversarial attacks [17, 13]
that are used to verify the robustness of our proposed SODEF. In the paper, we
use lowercase boldface characters like $\mathbf{z}$ to denote vectors in
$\mathbb{R}^{n}$, capital boldface characters like $\mathbf{A}$ to denote
matrices in $\mathbb{R}^{n\times n}$, and normal characters like $z$ to denote
scalars except that the notation $(x,y)$ are normal characters reserved to
denote the input and label pairs. A vector $\mathbf{z}\in\mathbb{R}^{n}$ is
represented as $(\mathbf{z}^{(1)},\mathbf{z}^{(2)},\ldots,\mathbf{z}^{(n)})$.
The $(i,j)$-th element of a matrix $\mathbf{A}$ is $\mathbf{A}_{ij}$ or
$[\mathbf{A}]_{ij}$. The Jacobian matrix of a function
$f:\mathbb{R}^{n}\mapsto\mathbb{R}^{n}$ evaluated at $\mathbf{z}$ is denoted
as $\nabla f(\mathbf{z})$. The set of functions
$\mathbb{R}^{n}\mapsto\mathbb{R}^{n}$ with continuous first derivatives is
denoted as $C^{1}(\mathbb{R}^{n},\mathbb{R}^{n})$.

##  2 Preliminaries: Neural ODE and Stability

In a neural ODE layer, the relation between the layer input $\mathbf{z}(0)$
and output $\mathbf{z}(T)$ is described as the following differential
equation:

$\displaystyle{\dfrac{\mathrm{d}{\mathbf{z}(t)}}{\mathrm{d}{t}}}=f_{\bm{\theta}}(\mathbf{z}(t),t)$
(1)

where $f_{\bm{\theta}}:\mathbb{R}^{n}\times[0,\infty)\mapsto\mathbb{R}^{n}$
denotes the non-linear trainable layers that are parameterized by weights
$\bm{\theta}$ and $\mathbf{z}:[0,\infty)\mapsto\mathbb{R}^{n}$ represents the
$n$-dimensional state of the neural ODE. Neural ODEs are the continuous analog
of residual networks where the hidden layers of residual networks can be
regarded as discrete-time difference equations
$\mathbf{z}(t+1)=\mathbf{z}(t)+f_{\bm{\theta}}(\mathbf{z}(t),t)$. In this
work, for simplicity, we only consider the time-invariant (autonomous) case
$f_{\bm{\theta}}(\mathbf{z}(t),t)=f_{\bm{\theta}}(\mathbf{z}(t))$, where the
dynamical system does not explicitly depend on $t$. For such non-linear
dynamical systems, the following theorem shows that under mild conditions, its
behaviour can be studied via linearization near special points called
hyperbolic equilibrium points.

######  Theorem 1 (Hartman–Grobman Theorem [18]).

Consider a system evolving in time with state $\mathbf{z}(t)\in\mathbb{R}^{n}$
that satisfies the differential equation
${\dfrac{\mathrm{d}{\mathbf{z}(t)}}{\mathrm{d}{t}}}=f(\mathbf{z}(t))$ for some
$f\in C^{1}(\mathbb{R}^{n},\mathbb{R}^{n})$,
$f(\mathbf{z})=(f^{(1)}(\mathbf{z}),\ldots,f^{(n)}(\mathbf{z}))$. Suppose the
map has a hyperbolic equilibrium state $\mathbf{z}^{*}\in\mathbb{R}^{n}$,
i.e., $f(\mathbf{z}^{*})=0$ and the Jacobian matrix $\nabla f=[\partial
f^{(i)}/\partial\mathbf{z}^{(j)}]_{i,j=1}^{n}$ of $f$ evaluated at
$\mathbf{z}=\mathbf{z}^{*}$ has no eigenvalue with real part equal to zero.
Then there exists a neighbourhood $N_{\mathbf{z}^{*}}$ of the equilibrium
point $\mathbf{z}^{*}$ and a homeomorphism
$g:N_{\mathbf{z}^{*}}\mapsto\mathbb{R}^{n}$, such that $g(\mathbf{z}^{*})=0$
and in the neighbourhood $N_{\mathbf{z}^{*}}$, the flow of
${\dfrac{\mathrm{d}{\mathbf{z}(t)}}{\mathrm{d}{t}}}=f(\mathbf{z}(t))$ is
topologically conjugate by the continuous map
$\bar{\mathbf{z}}(t)=g(\mathbf{z}(t))$ to the flow of its linearization
${\dfrac{\mathrm{d}{\bar{\mathbf{z}}(t)}}{\mathrm{d}{t}}}=\nabla
f(\mathbf{z}^{*})\cdot\bar{\mathbf{z}}(t)$.

The theorem states that when the Jacobian matrix at the zeros of $f$ has no
eigenvalue with zero real part, the behaviour of the original dynamical system
can be studied using the simpler linearization of the system around those
zeros. We next review some definitions and theorems from linear control theory
[10].

######  Definition 1 (Lyapunov Stability [10]).

The linear time-invariant system
${\dfrac{\mathrm{d}{\bar{\mathbf{z}}(t)}}{\mathrm{d}{t}}}=\mathbf{A}\bar{\mathbf{z}}(t)$
with constant matrix $\mathbf{A}$ is marginally stable or stable in the sense
of Lyapunov if every finite initial state $\bar{\mathbf{z}}(0)$ excites a
bounded response. It is asymptotically stable if every finite initial state
excites a bounded response, which, in addition, approaches $0$ as
$t\to\infty$.

######  Theorem 2 (Lyapunov Stability Theorem [10]).

a) The equation
${\dfrac{\mathrm{d}{\bar{\mathbf{z}}(t)}}{\mathrm{d}{t}}}=\mathbf{A}\bar{\mathbf{z}}(t))$
is marginally stable if and only if all eigenvalues of A have zero or negative
real parts and those with zero real parts are simple roots of the minimal
polynomial of $\mathbf{A}$. b) The equation
${\dfrac{\mathrm{d}{\bar{\mathbf{z}}(t)}}{\mathrm{d}{t}}}=\mathbf{A}\bar{\mathbf{z}}(t)$
is asymptotically stable if and only if all eigenvalues of A have negative
real parts.

In Theorem 1, we say that a hyperbolic equilibrium point is _Lyapunov-stable_
if all eigenvalues of the Jacobian matrix evaluated at it have negative real
parts. From Theorems 1 and 2, we see that a small perturbation around the
Lyapunov-stable equilibrium point $\mathbf{z}(0)$ leads to
$\tilde{\mathbf{z}}(t)\to\mathbf{z}(0)$ as $t\to\infty$, i.e.,
$\exists\>\delta>0$ such that for all $\tilde{\mathbf{z}}(0)$ with
$\lVert\mathbf{z}(0)-\tilde{\mathbf{z}}(0)\rVert_{2}<\delta$, we have
$\lVert\tilde{\mathbf{z}}(t)-\mathbf{z}(0)\rVert_{2}\to 0$ as $t\to\infty$,
where $\tilde{\mathbf{z}}(t)$ is the ODE solution for the perturbed input
$\tilde{\mathbf{z}}(0)$. In the context of neural network adversarial attacks,
if the malicious perturbations around the ODE input $\mathbf{z}(0)$ is small,
then the output $\mathbf{z}(T)$ for large enough $T$ will not be affected
significantly by the perturbation. Consequently, the succeeding network layers
after the neural ODE layer can still perform well without being affected by
the input perturbation. The perturbation weakening phenomenon around Lyapunov-
stable equilibrium points works like a noise filter and acts as a defense
against adversarial attacks.

We require the following definition and result in our stability analysis.

######  Definition 2 (Strictly diagonally dominant [19]).

Let $\mathbf{A}\in\mathbb{C}^{n\times n}$ . We say that $\mathbf{A}$ is
strictly diagonally dominant if $|\mathbf{A}_{ii}|>\sum_{j\neq
i}|\mathbf{A}_{ij}|$ for all $i=1,...,n$.

######  Theorem 3 (Levy–Desplanques theorem [19]).

If $\mathbf{A}\in\mathbb{C}^{n\times n}$ is strictly diagonally dominant and
if every main diagonal entry of A is real and negative, then $A$ is non-
singular and every eigenvalue of A has negative real part.

######  Lemma 1.

Given $k$ distinct points $\mathbf{z}_{i}\in\mathbb{R}^{n}$ and matrices
$\mathbf{A}_{i}\in\mathbb{R}^{n\times n}$, $i=1,...,k$, there exists a
function $f\in C^{1}(\mathbb{R}^{n},\mathbb{R}^{n})$ such that
$f(\mathbf{z}_{i})=0$ and $\nabla
f_{\bm{\theta}}(\mathbf{z}_{i})=\mathbf{A}_{i}$.

##  3 SODEF Architecture

We consider a classification problem with $L$ classes. The proposed SODEF
model architecture is shown in Fig. 1. The input $x\in X$ (e.g., an image) is
first passed through a feature extractor
$h_{\bm{\phi}}:X\mapsto\mathbb{R}^{n}$ to obtain an embedding feature
representation $\mathbf{z}(0)$. A neural ODE layer $f_{\bm{\theta}}$ follows
as a nonlinear feature mapping to stabilize the feature representation output
$\mathbf{z}(0)$ from $h_{\bm{\phi}}$. The final FC layer $\mathbf{V}$ serves
as a linear mapping to generate a prediction vector based on the output
$\mathbf{z}(T)$ of the neural ODE layer. The parameters
$\bm{\phi},\bm{\theta}$ and $\mathbf{V}$ are parameterized weights for the
feature extractor, neural ODE layer and FC layer, respectively.

We provide motivation and design guidance for the FC layer $\mathbf{V}$ in
Section 3.1, which attempts to separate Lyapunov-stable equilibrium points
implicitly by maximizing the similarity distance between feature
representations corresponding to the $L$ different classes. Experimental
results demonstrate the advantages of our diversity promoting FC layer in
Section 3.1 with comparisons to traditional neural ODEs without diversity
promoting.

However, the embedded features after using diversity promoting are not
guaranteed to locate near the Lyapunov-stable equilibrium points. In Section
3.2, we formulate an optimization problem to force embedding features to
locate near the Lyapunov-stable equilibrium points. We introduce optimization
constraints to force the Jacobian matrix of the ODE in our neural ODE layer to
have eigenvalues with negative real parts at the Lyapunov-stable equilibrium
points. Instead of directly imposing constraints on the eigenvalue of the
matrix, which may be computationally complex especially when the matrix is
large, we add constraints to the matrix elements instead.

$\displaystyle x$$\displaystyle y$feature extractor $\displaystyle
h_{\bm{\phi}}$ODE
${\dfrac{\mathrm{d}{\mathbf{z}(t)}}{\mathrm{d}{t}}}=f_{\bm{\theta}}(\mathbf{z}(t))$FC
$\mathbf{V}$$\mathbf{z}(0)$$\mathbf{z}(T)$ Fig. 1: SODEF model architecture.

###  3.1 Maximizing the Distance between Lyapunov-Stable Equilibrium Points

From Section 2, we observe that points in a small neighbourhood of a Lyapunov-
stable equilibrium point is robust against adversarial perturbations. We call
this neighborhood a _stable neighborhood_. However Lyapunov-stable equilibrium
points for different classes may very well locate near each other and
therefore each stable neighborhood may be very small, leading to poor
adversarial defense. In this section, we propose to add a FC layer after the
neural ODE layer given by 1 to avoid this scenario. The purpose of the FC
layer is to map the output of the neural ODE layer to a feature vector
$\mathbf{v}_{l}$ if the input $x$ belongs to the class $l=1,\ldots,L$. We
design the FC layer so that the cosine similarities between different
$\mathbf{v}_{l}$’s are minimized.

######  Lemma 2.

Given a set of $k$ unit vectors $\mathbf{v}_{1},\ldots,\mathbf{v}_{k}$ in
$\mathbb{R}^{n}$, where $n\geq k$, let
$a(\mathbf{v}_{1},\ldots,\mathbf{v}_{k})=\max_{i\neq j}\
\mathbf{v}_{i}^{\intercal}\mathbf{v}_{j}$. Then $\min
a(\mathbf{v}_{1},\ldots,\mathbf{v}_{k})=1/(1-k)$, where the minimum is taken
over all possible sets of $k$ unit vectors
$\mathbf{v}_{1},\ldots,\mathbf{v}_{k}$.

######  Corollary 1.

Consider a $k\times k$ matrix $\mathbf{B}=[b_{ij}]_{i,j=1}^{k}$ with
$b_{ii}=1$ and $b_{ij}=1/(1-k),\ \forall\,i\neq j$. Let the eigen
decomposition of $\mathbf{B}$ be
$\mathbf{B}=\mathbf{U}\bm{\Sigma}\mathbf{U}^{\intercal}$. For any $n\geq k$
and $i=1,\ldots,k$, let $\mathbf{v}_{i}$ be the $i$-th column of
$\mathbf{Q}\bm{\Sigma}^{1/2}\mathbf{U}^{\intercal}$, where $\mathbf{Q}$ is any
$n\times k$ matrix such that
$\mathbf{Q}^{\intercal}\mathbf{Q}=\mathbf{I}_{k}$. Then,
$a(\mathbf{v}_{1},\ldots,\mathbf{v}_{k})=\max_{i\neq
j}\mathbf{v}_{i}^{\intercal}\mathbf{v}_{j}=1/(1-k)$.

Corollary 1 suggests a diversity promoting scheme to maximally separate the
equilibrium points of the neural ODE layer. The FC layer is represented by an
$n\times L$ matrix $\mathbf{V}=[\mathbf{v}_{1},\ldots,\mathbf{v}_{L}]$, where
$n$ is the dimension of $\mathbf{z}(T)$, the output from the neural ODE layer.
If $\mathbf{z}(T)$ is generated from an input from class $l$, it is mapped to
$\mathbf{v}_{l}$. By minimizing the maximum cosine similarity
$a(\mathbf{v}_{1},\ldots,\mathbf{v}_{k})=\max_{i\neq
j}\mathbf{v}_{i}^{\intercal}\mathbf{v}_{j}$ between the representations from
two different classes, we ensure that the output of SODEF is robust to
perturbations in the input. Corollary 1 provides a way to choose the FC layer
weights $\mathbf{V}$.

To validate our observations, we conduct experiments to compare the robustness
of ODE net [9] and TisODE [6] with and without our proposed FC layer
$\mathbf{V}$, on two standard datasets: MNIST [2] and CIFAR10 [12] 111Our
experiments are run on a GeForce RTX 2080 Ti GPU.. On the MNIST dataset, all
models consist of four convolutional layers and one fully-connected layer. On
the CIFAR10 dataset, the networks are similar to those for MNIST except the
down-sampling network is a stack of 2 ResNet blocks. In practice, the neural
ODE can be solved with different numerical solvers such as the Euler method
and the Runge-Kutta methods [9]. Here, we use Runge-Kutta of order 5 in our
experiments. Our implementation builds on the open-source neural ODE
codes.222https://github.com/rtqichen/torchdiffeq During training, no Gaussian
noise or adversarial examples are augmented into the training set. We test the
performance of our model in defending against white-box attacks FGSM [17] and
PGD [13] . The parameters for different attack methods used in this paper are
given in the supplementary material. From Tables 1 and 2, we observe that for
both datasets, our fixed FC layer improves each network’s defense ability by a
significant margin. We visualize the features before the final FC layer using
t-SNE [20] in Figs. 2 and 3. We observe that with the FC layer, the features
for different classes are well separated even under attacks.

Table 1: Classification accuracy (%) on adversarial MNIST examples, where the superscript + indicates the last FC layer is fixed to be $\mathbf{V}$. Attack | Para. | ODE | ODE+ | TisODE | TisODE+
---|---|---|---|---|---
None | - | 99.6 | 99.7 | 99.5 | 99.7
FGSM | $\epsilon=0.3$ | 31.4 | 52.8 | 45.9 | 63.5
PGD | $\epsilon=0.3$ | 0.29 | 0.30 | 0.4 | 20.20
Table 2: Classification accuracy (%) on adversarial CIFAR10 examples, where the superscript + indicates the last FC layer is fixed to be $\mathbf{V}$. Attack | Para. | ODE | ODE+ | TisODE | TisODE+
---|---|---|---|---|---
None | - | 87.0 | 85.0 | 87.4 | 81.8
FGSM | $\epsilon=0.1$ | 12.9 | 47.6 | 13.1 | 41.9
PGD | $\epsilon=0.1$ | 7.8 | 14.7 | 7.4 | 16.2

Fig. 2: t-SNE visualization results on the features before the final FC layer.
The input is the test set of MNIST. Left: trained with TisODE, middle: TisODE
using a randomly chosen orthogonal matrix as the final FC, right: TisODE using
proposed $\mathbf{V}$ as the final FC.

Fig. 3: t-SNE visualization results on the features before the final FC layer.
The input is the adversarial examples of the test set of MNIST generated using
FGSM method at $\epsilon=0.3$. Left: trained with TisODE, middle: TisODE using
a randomly chosen orthogonal matrix as the final FC, right: TisODE using
proposed $\mathbf{V}$ as the final FC.

###  3.2 Objective Formulation and Stability

In this subsection, we formulate an optimization framework for SODEF to force
output features to locate within the stable neighborhood of Lyapunov-stable
equilibrium points. We make the following assumption.

######  Assumption 1.

The input $x$ takes values in a compact metric space $X$ and has probability
distribution $\mu$. The feature extractor $h_{\bm{\phi}}$ is injective and
continuous.

The above assumption is satisfied if the input $x$ (e.g., an image) resides in
a bounded and closed set of a Euclidean space. We denote the pushforward
measure (still a probability distribution) of $\mu$ under the continuous
feature extractor mapping $h_{\bm{\phi}}$ as $\nu_{\bm{\phi}}=\mu\circ
h_{\bm{\phi}}^{-1}$, where $\circ$ denotes function composition. The
conditional probability distribution for the embedding of each class
$l\in\\{1,...,L\\}$ has compact support $E_{l}\subset\mathbb{R}^{n}$ since
$E_{l}$ is closed and $h_{\bm{\phi}}(X)$ is bounded in $\mathbb{R}^{n}$. In
Section 3.1, the FC layer $\mathbf{V}$ tries to maximize the distance between
$E_{l}$, $l=1,\ldots,L$. In this section for analysis purposes, we also assume
the following.

######  Assumption 2.

We have $E_{l}\bigcap E_{l^{\prime}}=\emptyset$ if $l\neq l^{\prime}$, i.e.,
the supports of each class are pairwise disjoint.

Our objective function is formulated as follows, which is explained in detail
in the sequel:

$\displaystyle\min_{\bm{\theta},\bm{\phi}}$
$\displaystyle~{}\mathbb{E}_{\mu}\ell(\mathbf{V}^{\intercal}(\mathbf{z}(T)),y_{i})$
(2) $\displaystyle\operatorname{s.t.\ }$
$\displaystyle~{}\mathbb{E}_{\nu_{\bm{\phi}}}\lVert
f_{\bm{\theta}}(\mathbf{z}(0))\rVert_{2}<\epsilon,\ f_{\bm{\theta}}\in
C^{1}(\mathbb{R}^{n},\mathbb{R}^{n}),$ (3)
$\displaystyle~{}\mathbb{E}_{\nu_{\bm{\phi}}}\Big{[}\left\lvert[\nabla
f_{\bm{\theta}}(\mathbf{z}(0))]_{ii}\right\rvert-\sum_{j\neq
i}\left\lvert[\nabla
f_{\bm{\theta}}(\mathbf{z}(0))]_{ij}\right\rvert\Big{]}>0,\
\forall\,i=1,\dots,n,$ (4)
$\displaystyle~{}\mathbf{z}(0)=h_{\bm{\phi}}(x),\text{ and
}\mathbf{z}(T)\text{ is the output of \lx@cref{creftype~refnum}{eq:NODE} with
input }\mathbf{z}(0).$ (5)

Here, $\ell$ is a loss function and $\epsilon>0$ is a positive constant. The
constraints 3, LABEL:eq:re2 and 4 force $\mathbf{z}(0)$ to be near the
Lyapunov-stable equilibrium points with strictly diagonally dominant
derivatives. We limit the $f_{\bm{\theta}}$ to be in
$C^{1}(\mathbb{R}^{n},\mathbb{R}^{n})$ to satisfy the condition in Theorem 1.
From [21], we also know that standard multi-layer feed forward networks with
as few as a single hidden layer and arbitrary bounded and non-constant
activation function are universal approximators for
$C^{1}(\mathbb{R}^{n},\mathbb{R}^{n})$ functions with respect to some
performance criteria provided only that sufficiently many hidden units are
available.

As a comparison, TisODE [6] only includes a constraint similar to 3, which in
general provides no guarantee to force $\mathbf{z}(0)$ near the Lyapunov-
stable equilibrium points. In the extreme case with parameters $\bm{\theta}=0$
for $f_{\bm{\theta}}$ such that $f_{\bm{\theta}}=0$, the ODE degenerates to an
identity mapping. No $\mathbf{z}(0)\in\mathbb{R}^{n}$ can now be a Lyapunov-
stable equilibrium point, and no stability can therefore be guaranteed to
defend against adversarial attacks even though the ODE curves still possess
the non-intersecting property and steady-state constraint, which were cited as
reasons for the stability of TisODE.

Instead of directly optimizing the above objective function, in our
implementation, we optimize the following empirical Lagrangian with a training
set $\\{(x_{k},y_{k}):k=1,...,N\\}$:

$\displaystyle\min_{\bm{\theta},\bm{\phi}}$ $\displaystyle\
\frac{1}{N}\sum_{k=0}^{N-1}\left(\ell\big{(}\mathbf{V}^{\intercal}\mathbf{z}_{k}(T),y_{k}\big{)}+\alpha_{1}\lVert
f_{\bm{\theta}}\big{(}\mathbf{z}_{k}(0)\big{)}\rVert_{2}+\alpha_{2}g_{1}\Big{(}\sum_{i=1}^{n}[\nabla
f_{\bm{\theta}}(\mathbf{z}_{k}(0))]_{ii}\Big{)}\right.$
$\displaystyle\quad+\left.\alpha_{3}g_{2}\Big{(}\sum_{i=1}^{n}(-\left\lvert[\nabla
f_{\bm{\theta}}(\mathbf{z}_{k}(0))]_{ii}\right\rvert+\sum_{j\neq
i}\left\lvert[\nabla
f_{\bm{\theta}}(\mathbf{z}_{k}(0))]_{ij}\right\rvert)\Big{)}\right)$ (6)
$\displaystyle\operatorname{s.t.\ }$ $\displaystyle\
\mathbf{z}_{k}(0)=h_{\bm{\phi}}(x_{k}),\text{ and }\mathbf{z}_{k}(T)\text{ is
the output of \lx@cref{creftype~refnum}{eq:NODE} with input
}\mathbf{z}_{k}(0),\ \forall\,k=1,\ldots,N$ (7)

where $\alpha_{1}$, $\alpha_{2}$ and $\alpha_{3}$ are hyperparameter weights,
$g_{1}$ and $g_{2}$ are chosen monotonically increasing functions bounded
below to eliminate the unbounded impact of the two regularizers that can
otherwise dominate the loss. In this paper, we set
$g_{1}(\cdot)=g_{2}(\cdot)=\text{exp}(\cdot)$. We call these two latter terms
the SODEF regularizers.

Suppose for each class $l=1,\ldots,L$, the embedding feature set
$E_{l}=\\{\mathbf{z}^{(l)}_{1},\ldots,\mathbf{z}^{(l)}_{k}\\}$ is finite. For
each $i=1,\ldots,k$, let $\mathbf{A}_{i}\in\mathbb{R}^{n\times n}$ be strictly
diagonally dominant matrix with every main diagonal entry be negative such
that the eigenvalues for $\mathbf{A}_{i}$ all have negative real part. From
Theorem 3, each $\mathbf{A}_{i}$ is non-singular and every eigenvalue of
$\mathbf{A}_{i}$ has negative real part. Therefore, from Theorem 2 and Lemma
1, there exists a function $f_{\bm{\theta}}$ such that all
$\mathbf{z}^{(l)}_{i}$ are Lyapunov-stable equilibrium points with
corresponding first derivative $\nabla
f_{\bm{\theta}}(\mathbf{z}^{(l)}_{i})=\mathbf{A}_{i}$. This shows that if
there exist only finite representation points for each class, we can find a
function $f_{\bm{\theta}}$ such that all inputs to the neural ODE layer are
Lyapunov-stable equilibrium points for $f_{\bm{\theta}}$ and

  1. (a)

$\mathbb{E}_{\nu_{\bm{\phi}}}\lVert
f_{\bm{\theta}}(\mathbf{z}(0))\rVert_{2}=0$,

  2. (b)

$\mathbb{E}_{\nu_{\bm{\phi}}}\left[{\nabla
f_{\bm{\theta}}(\mathbf{z}(0))}\right]_{ii}<0,\ \forall\,i=1,\dots,n$,

  3. (c)

$\mathbb{E}_{\nu_{\bm{\phi}}}\left[{\left\lvert[\nabla
f_{\bm{\theta}}(\mathbf{z}(0))]_{ii}\right\rvert-\sum_{j\neq
i}\left\lvert[\nabla
f_{\bm{\theta}}(\mathbf{z}(0))]_{ij}\right\rvert}\right]>0,\
\forall\,i=1,\dots,n$.

If the input space $X$ has infinite cardinality, then an injective and
continuous feature extractor $h_{\bm{\phi}}$ results in a $\nu_{\bm{\phi}}$
with non-finite support, i.e., at least one $E_{l}$, $l=1,\ldots,L$, is
infinite. It is not obvious whether we can obtain a $f_{\bm{\theta}}$ where
every point in $E=\bigcup_{l}E_{l}$ is a stable equilibrium point. The
following result gives a negative answer if $\nu_{\bm{\phi}}$ is a continuous
measure (i.e., absolutely continuous with respect to (w.r.t.) Lebesgue
measure) on some subset.

######  Lemma 3.

If the restriction of $\nu_{\bm{\phi}}$ to some open set $E^{\prime}\subset E$
is a continuous measure, there is no continuous function $f_{\bm{\theta}}$
such that for $\nu_{\bm{\phi}}$-almost surely all $\mathbf{z}\in E$,
$f_{\bm{\theta}}(\mathbf{z})=0$ and all the eigenvalues of $\nabla
f_{\bm{\theta}}(\mathbf{z})$ have negative real parts. In other words, there
is no continuous function $f_{\bm{\theta}}$ such that almost surely all
$\mathbf{z}$ in $E$ are Lyapunov-stable equilibrium points.

Lemma 3 indicates that it is too much to hope for all points in $E$ to be
Lyapunov-stable equilibrium points. In the following, we relax this
requirement and show that under mild conditions, for all $\epsilon>0$, we can
find a continuous function $f_{\bm{\theta}}$ with finitely many stable
equilibrium points such that Items b and c above hold and Item a is replaced
by $\mathbb{E}_{\nu_{\bm{\phi}}}\lVert
f_{\bm{\theta}}(\mathbf{z}(0))\rVert_{2}<\epsilon$. This motivates the
optimization constraints in 3, LABEL:eq:re2 and 4.

######  Theorem 4.

Suppose Assumptions 1 and 2. If $\nu_{\bm{\phi}}$ is not a continuous uniform
measure on $E_{l}$ for each $l=1,\ldots,L$, then the following holds: 1) The
function space satisfying the constraints in 3, LABEL:eq:re2 and 4 is non-
empty for all $\epsilon>0$. 2) If additionally the restriction of
$\nu_{\bm{\phi}}$ to any open set $O\subset E_{l}$ is not a continuous uniform
measure, there exist functions in this space such that each support $E_{l}$
contains at least one Lyapunov-stable equilibrium point.

##  4 Experiments

In this section, we evaluate the robustness of SODEF under adversarial attacks
with different attack parameters. We conduct experiments to compare the
robustness of ODE net [9] and TisODE net [6] on three standard datasets: MNIST
[2], CIFAR10 and CIFAR100 [3]. Since SODEF is compatible with many defense
methods, it can be applied to any neural network’s final regressor layer to
enhance its stability against adversarial attacks. Our experiment codes are
provided in https://github.com/KANGQIYU/SODEF.

###  4.1 Setup

We use open-source pre-trained models that achieve the top accuracy on each
dataset as the feature extractor $h_{\bm{\phi}}$. Specifically for simple
MNIST task, we use the ResNet18 model provided in Pytorch. We use the model
provided by [22], which obtains nearly $88\%$ clean accuracy on CIFAR100 using
EfficientNet [23] and the model provided by [24], which has nearly $95\%$
clean accuracy on CIFAR10. In the neural ODE layer, $f_{\bm{\theta}}$ consists
of 2 FC layers. During the trainings of SODEF (except in the experiment
included in Section 4.2), we train the neural network with the fixed FC
introduced in Section 3.1. In the first $30$ epochs, we fixed
$f_{\bm{\theta}}$ to let the feature extractor $h_{\bm{\phi}}$ learn a feature
representation with only the cross-entropy loss $\ell$, and in the remaining
$120$ epochs, we release $h_{\bm{\phi}}$ to further train $f_{\bm{\theta}}$
using Section 3.2 with $\alpha_{1}=1$ and $\alpha_{2}=\alpha_{3}=0.05$. For
CIFAR10 and CIFAR100, the pixel values are normalized by $(x-\mu)/\sigma$
where $\mu=[0.4914,0.4822,0.4465]$ and $\sigma=[0.2023,0.1994,0.2010]$ 333To
test AutoAttack, we have strictly followed the instruction in
https://github.com/RobustBench/robustbench to attack the original images
before any normalization or resizing.. To show that our SODEF is compatible
with many defense methods and can be applied to any neural network’s final
regression layer, we conduct an experiment where we use a recently proposed
robust network TRADES [25] as the feature extractor in our SODEF. The
pretrained model is provided here 444https://github.com/P2333/Bag-of-Tricks-
for-AT, and we choose the model with architecture "WRN-34-10" to conduct our
experiments. Besides the two vanilla white-box attacks FGSM and PGD as
metioned in Section 3.1, we also include a strong ensemble attack AutoAttack
[26], which sequentially performs attack using all of the following four
individual attacks: three white-box attacks APGD${}_{\text{CE}}$,
APGD${}^{\text{T}}_{\text{DLR}}$ and FAB${}^{\text{T}}$[27], and one black-box
Square attack [28]. We refer the reader to the the supplementary material for
more details of the attacks used in this paper, where, in additional, more
experiments are included.

###  4.2 Compatibility of SODEF

Adversarial training (AT) is one of the most effective strategies for
defending adversarial attacks. TRADES [25] is one of the adversarial training
defense methods with combinations of tricks of warmup, early stopping, weight
decay, batch size and other hyper parameter settings. In this experiment we
fix the pretained TRADES model (except the final FC layer (size 640x10)) as
our feature extractor $h_{\phi}$. We then append our (trainable) SODEF with
integration time $T=5$ to the output of the feature extractor. To evaluate
model robustness, we use AutoAttack and attack the models using both the
$\mathcal{L}_{2}$ norm ($\epsilon=0.5$) and $\mathcal{L}_{\infty}$ norm
($\epsilon=8/255$). The results are shown in Table 3. We clearly observe that
our SODEF can enhance TRADES’s robustness under all the four individual
attacks and the strongest ensemble AutoAttack. For the strong
$\mathcal{L}_{2}$ AutoAttack, our SODEF have improved the model robustness
from $59.42\%$ to $67.75\%$. Our experiment show that SODEF can be applied to
many defense models’ regression layer to enhance their stability against
attacks.

Table 3: Classification accuracy (%) using TRADES (w/ and w/o SODEF) under AutoAttack on adversarial CIFAR10 examples with $\mathcal{L}_{2}$ norm ($\epsilon=0.5$) and $\mathcal{L}_{\infty}$ norm ($\epsilon=8/255$). Attack / Model |  TRADES $\mathcal{L}_{\infty}$ |  TRADES+SODEF $\mathcal{L}_{\infty}$ |  TRADES $\mathcal{L}_{2}$ |  TRADES+SODEF $\mathcal{L}_{2}$
---|---|---|---|---
Clean | 85.48 | 85.18 | 85.48 | 85.18
APGD${}_{\text{CE}}$ | 56.08 | 70.90 | 61.74 | 74.35
APGD${}^{\text{T}}_{\text{DLR}}$ | 53.70 | 64.15 | 59.22 | 68.55
FAB${}^{\text{T}}$ | 54.18 | 82.92 | 60.31 | 83.15
Square | 59.12 | 62.21 | 72.65 | 76.02
AutoAttack | 53.69 | 57.76 | 59.42 | 67.75

###  4.3 Influence of Integration Time $T$

From the discussion after Theorems 1 and 2, we know if the malicious
perturbations around the ODE input Lyapunov-stable equilibrium point
$\mathbf{z}(0)$ is small, then the output $\mathbf{z}(T)$ for large enough $T$
will not be affected significantly by the perturbation:
$\lVert\tilde{\mathbf{z}}(t)-\mathbf{z}(0)\rVert_{2}\to 0\text{ as
}t\to\infty$. Consequently, the succeeding network layers after the neural ODE
layer can still perform well without being affected by the input perturbation.
In this section, we test the influence of the SODEF integration time $T$ using
CIFAR100. We use the model EfficientNet provided by [23] as $h_{\phi}$ (Note,
unlike Section 4.2, $h_{\phi}$ is trainable in this experiments). We use
AutoAttack with $\mathcal{L}_{2}$ norm ($\epsilon=0.5$). We observe that for
all the four individual attacks and the strongest ensemble AutoAttack, SODEF
performs generally better for large integration time $T$. We also test larger
integration time $T>10$, but do not see any obvious improvements.

Table 4: Classification accuracy (%) under AutoAttack on adversarial CIFAR100 examples with $\mathcal{L}_{2}$ norm, $\epsilon=0.5$ and different integration time $T$ for SODEF. Attack / $T$ | 1 | 3 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---
Clean | 88.00 | 88.12 | 88.15 | 88.00 | 87.92 | 88.00 | 88.05 | 88.10
APGD${}_{\text{CE}}$ | 17.20 | 21.33 | 21.05 | 23.67 | 69.67 | 85.33 | 87.10 | 86.88
APGD${}^{\text{T}}_{\text{DLR}}$ | 21.02 | 21.00 | 22.00 | 26.00 | 63.30 | 86.90 | 86.20 | 86.54
FAB${}^{\text{T}}$ | 86.33 | 85.10 | 86.36 | 87.70 | 87.67 | 86.55 | 86.22 | 85.93
Square | 84.67 | 86.22 | 87.05 | 87.20 | 86.90 | 86.33 | 87.05 | 86.75
AutoAttack | 2.00 | 3.53 | 4.87 | 4.33 | 30.66 | 78.80 | 78.97 | 79.10

###  4.4 Performance Comparison Under AutoAttack

For a comparison, we provide the results of applying AutoAttack to other
baseline models mentioned in the paper. We set the same integration time for
ODE, TisODE and SODEF. We observe that for the strongest AutoAttack, our SODEF
outperforms the other baseline models by a significant margin. In this case,
SODEF achieves $79.10\%$ accuracy while other models only get less than $3\%$
accuracy.

Table 5: Classification accuracy (%) under AutoAttack on adversarial CIFAR100 examples with $\mathcal{L}_{2}$ norm, $\epsilon=0.5$ and $T=10$. Attack / Model | No ODE | ODE | TisODE | SODEF
---|---|---|---|---
Clean | 88.00 | 87.90 | 88.00 | 88.10
APGD${}_{\text{CE}}$ | 23.30 | 6.75 | 14.32 | 86.88
APGD${}^{\text{T}}_{\text{DLR}}$ | 7.33 | 22.00 | 24.20 | 86.54
FAB${}^{\text{T}}$ | 79.30 | 78.67 | 77.16 | 85.93
Square | 84.52 | 85.67 | 86.32 | 86.75
AutoAttack | 0.00 | 1.33 | 4.06 | 79.10

###  4.5 Performance Under PGD and FGSM Attacks

White-box adversaries have knowledge of the classifier models, including
training data, model architectures and parameters. We test the performance of
our model in defending against the white-box attacks, PGD and FGSM. We set
$T=5$ as the integration time for the neural ODE layer. The parameters for
different attack methods used are given in the supplementary material. The
subsequent experiments use these settings by default, unless otherwise stated.

Table 6: Classification accuracy (%) on adversarial MNIST examples. Attack | Para. | no ode | ODE | TisODE | SODEF
---|---|---|---|---|---
None | - | 99.45 | 99.42 | 99.43 | 99.44
FGSM | $\epsilon=0.3$ | 10.03 | 29.6 | 36.70 | 63.36
PGD | $\epsilon=0.3$ | 0.31 | 1.56 | 1.82 | 45.25

The classification results on MNIST are shown in Table 6. We observe that
while maintaining the state-of-the-art accuracy on normal images, SODEF
improves the adversarial robustness as compared to the other two methods. For
the most effective attack in this experiment, i.e., PGD attack, SODEF shows a
$45.25\%-1.56\%=43.69\%$ improvement over ODE and a $45.25\%-1.23\%=44.02\%$
improvement over TisODE.

Table 7: Classification accuracy (%) on adversarial CIFAR10 examples. Attack | Para. | no ode | ODE | TisODE | SODEF
---|---|---|---|---|---
None | - | 95.2 | 94.9 | 95.1 | 95.0
FGSM | $\epsilon=0.1$ | 47.31 | 45.23 | 43.28 | 68.05
PGD | $\epsilon=0.1$ | 3.09 | 3.21 | 3.80 | 55.59

For CIFAR-10, we see from Table 7 that SODEF maintains high accuracy on normal
examples and makes the best predictions under adversarial attacks. In
particular, SODEF achieves an absolute percentage point improvement over ODE
net up to $52.38\%$ and over TisODE up to $52.54\%$ for PGD attack.

For CIFAR-100, the results in the supplementary material shows that the most
effective attack causes the classification accuracy to drop relatively by
$74.6\%=\frac{88.0-22.35}{88.0}$ for SODEF and by
$97.3\%=\frac{88.3-2.39}{88.3}$ for vanilla EfficientNet, which is pre-trained
on ImageNet to obtain a top clean accuracy. Neither ODE net nor TisODE net can
improve the classification accuracy under PGD attack by a big margin, e.g.
TisODE net only improves the classification accuracy from $2.39\%$ to
$3.44\%$, while SODEF still shows clear defense capability in this scenario.

###  4.6 Ablation Studies

The impact of the ODE with and without the SODEF regularizers in Section 3.2
has been presented in the above comparisons between SODEF and ODE. In this
section, we show the necessity of diversity promoting using the FC introduced
in Section 3.1 and conduct transferability study.

####  4.6.1 Impact of Diversity Promotion

Table 8: Classification accuracy (%) on adversarial MNIST examples, where the superscript - indicates the last FC layer is not fixed to be $\mathbf{V}$ and is set to be a trainable layer.  Attack | Para. | SODEF | SODEF-
---|---|---|---
None | - | 95.0 | 95.1
FGSM | $\epsilon=0.1$ | 63.36 | 51.6
PGD | $\epsilon=0.1$ | 45.25 | 34.9

Table 8 shows the difference of the defense performance when fixing the final
FC be $\mathbf{V}$ or setting it to a trainable linear layer. It can be seen
that having diversity control improves the robustness. One possible reason for
this phenomenon given in Section 3 is that diversity promotion with a fixed
designed FC attempts to make the embedding feature support $E_{l}$ of each
class $l$ disjoint to each other and therefore the Lyapunov-stable equilibrium
points for each $E_{l}$ are well separated.

####  4.6.2 Transferability Study

Transferability study is carried out on CIFAR-10, where the adversarial
examples are generated using FGSM and PGD attacks using ResNet18 without any
ODEs. The classification accuracy drops from $68.05\%$ to $59\%$ for FGSM with
$\epsilon=0.3$, and from $55.59\%$ to $34\%$ for PGD with $\epsilon=0.1$. One
possible reason for this phenomenon is that ODEs have obfuscated gradient
masking effect as discussed in [15], and a transfer attack may deteriorate the
defense effect. However, as we observe from Table 7, even with a transfer
attack on SODEF, it still performs better than other ODEs without transfer
attacks.

##  5 Conclusion

In this paper, we have developed a new neural ODE network, SODEF, to suppress
input perturbations. SODEF is compatible with any existing neural networks and
can thus be appended to the state-of-the-art networks to increase their
robustness to adversarial attacks. We demonstrated empirically and
theoretically that the robustness of SODEF mainly derives from its stability
and proposed a training method that imposes constraints to ensure all
eigenvalues of the Jacobian matrix of the neural ODE layer have negative real
parts. When each classification class converges to its own equilibrium points,
we showed that the last FC layer can be designed in such a way that the
distance between the stable equilibrium points is maximized, which further
improves the network’s robustness. The effectiveness of SODEF has been
verified under several popular while-box attacks.

## Acknowledgments and Disclosure of Funding

This research is supported in part by A*STAR under its RIE2020 Advanced
Manufacturing and Engineering (AME) Industry Alignment Fund – Pre Positioning
(IAF-PP) (Grant No. A19D6a0053) and the RIE2020 Industry Alignment Fund –
Industry Collaboration Projects (IAF-ICP) Funding Initiative, as well as cash
and in-kind contribution from the industry partner(s). The computational work
for this article was partially performed on resources of the National
Supercomputing Centre, Singapore (https://www.nscc.sg).

## Broader Impact

Our work, which contributes to more robust DNNs, is supposed to mitigate the
threat of adversarial attacks. However, on the hand, the reliable deployment
of DNNs in automation of tasks will potentially bring mass-scale unemployment
and social unrest. As DNNs become more robust and more tasks, especially those
whose failures will bring high risks to human lives or large property losses
under adversarial attacks, fall into the automatic task category, massive jobs
could disappear.

##  Appendix A Related Work and Adversarial Attacks

In this section, we give an overview of related work in stable neural ODE
networks. We also give an overview of common adversarial attacks and recent
works that defend against adversarial examples.

### Stable Neural Network

Gradient vanishing and gradient exploding are two well-known phenomena in deep
learning [29]. The gradient of the objective function, which strongly relies
on the training method as well as the neural network architecture, indicates
how sensitive the output is w.r.t. input perturbation. Exploding gradient
implies instability of the output w.r.t. the input and thus resulting in a
non-robust learning architecture. On the other hand, vanishing gradient
implies insensitivity of the output w.r.t. the input, i.e., robustness against
input perturbation. However, this prohibits effective training. To overcome
these issues, [14, 15] proposed a new DNN architecture inspired by ODE
systems. This DNN architecture is stable in the sense that the input-output of
the linearized system is always norm-preserving. This is different from the
Lyapunov stability [10] our paper pursues.

The reference [16] developed an implicit Euler based skip-connection to
enhance ResNets with better stability and adversarial robustness. More
recently, [6] explored the robustness of neural ODE and related the robustness
of neural ODE against input perturbation with the non-intersecting property of
ODE integral curves. In addition, [6] proposed TisODE, which regularizes the
flow on perturbed data via the time-invariant property and the imposition of a
steady-state constraint. However, in general, neither the non-intersecting
property nor the steady-state constraint can guarantee robustness against
input perturbation. In this paper, we propose to impose constraints to ensure
Lyapunov stability during training so that adversarial robustness can be
enhanced.

### Adversarial Examples

Recent works [5, 30, 31] demonstrated that in image classification tasks, the
input images can be modified by an adversary with human-imperceptible
perturbations so that a DNN is fooled into mis-classifying them. In [32, 33],
the authors proposed targeted attacks for audio waveform tasks. They showed
that adversarial audio can be embedded into speech so that DNNs cannot
recognize the input as human speech. Furthermore, adversarial speech can fool
speech-to-text DNN systems into transcribing the input into any pre-chosen
phrase. The references [34, 35] proposed to generate targeted unnoticeable
adversarial examples for 3D point clouds, which threatens many safety-critical
applications such as autonomous driving.

Adversarial attacks [17, 36, 13, 37, 38] can be grouped into two categories:
white-box attacks where adversaries have knowledge of the classifier models,
including training data, model architectures and parameters, and black-box
attacks where adversaries do not know the model’s internal architecture or
training parameters. An adversary crafts adversarial examples based on a
substitute model and then feed these examples to the original model to perform
the attack.

To mitigate the effect of adversarial attacks, many defense approaches have
been proposed such as adversarial training [5, 17], defensive distillation
[39], ECOC based methods [40, 41], and post-training defenses [42, 43, 44].
None of these methods make use of the stability of ODEs.

For a normal image-label pair $(x,y)$ and a trained DNN $f_{\bm{\theta}}$ with
$\bm{\theta}$ being the vector of trainable model parameters, an adversarial
attack attempts to find an adversarial example $x^{\prime}$ that remains
within the $\mathcal{L}_{p}$-ball with radius $\epsilon$ centered at the
normal example $x$, i.e., $\|x-x^{\prime}\|_{p}\leq\epsilon$, such that
$f_{\bm{\theta}}(x^{\prime})\neq y$. In what follows, we briefly present some
popular adversarial attacks that are used to verify the robustness of our
proposed SODEF.

Fast Gradient Sign Method (FGSM) [17] perturbs a normal input $x$ in its
$\mathcal{L}_{\infty}$ neighborhood to obtain

$\displaystyle x^{\prime}=x+\epsilon\cdot{\rm
sign}\left(\nabla_{x}\mathcal{L}(f_{\bm{\theta}}(x),y)\right),$

where $\mathcal{L}(f_{\bm{\theta}}(x),y)$ is the cross-entropy loss of
classifying $x$ as label $y$, $\epsilon$ is the perturbation magnitude, and
the update direction at each pixel is determined by the sign of the gradient
evaluated at this pixel. FGSM is a simple yet fast and powerful attack.

Projected Gradient Descent (PGD) [13] iteratively refines FGSM by taking
multiple smaller steps $\alpha$ in the direction of the gradient. The
refinement at iteration $i$ takes the following form:

$\displaystyle x^{\prime}_{i}=x^{\prime}_{i-1}+{\rm
clip}_{\epsilon,x}\left(\alpha\cdot{\rm
sign}\left(\nabla_{x}\mathcal{L}(f_{\bm{\theta}}(x),y)\right)\right),$

where $x^{\prime}_{0}=x$ and ${\rm clip}_{\epsilon,x}(x^{\prime})$ performs
clipping of the input $x^{\prime}$. For example, in an image $x^{\prime}$, if
a pixel takes values between 0 and 255, ${\rm
clip}_{\epsilon,x}(x^{\prime})=\min\left\\{255,x+\epsilon,\max\\{0,x-\epsilon,x^{\prime}\\}\right\\}$.

AutoAttack: In [26], the authors proposed two extensions of the PGD-attack,
named APGD${}_{\text{CE}}$ and APGD${}^{\text{T}}_{\text{DLR}}$, to overcome
failures of vanilla PGD due to suboptimal step size and problems of the
objective function. They combined these two attacks with two complementary
attacks called (white-box) FAB${}^{\text{T}}$[27] and (black-box) Square[28]
to form a parameter-free ensemble of attacks, called AutoAttack, to test
adversarial robustness.

##  Appendix B Attack Parameters

For two vanilla white-box FGSM and PGD attacks performed in Sections 3.1 and
4, we use $\mathcal{L}_{\infty}$ norm with maximum perturbation $\epsilon=0.3$
for MNIST and $\epsilon=0.1$ for CIFAR10 and we iterates PGD for 20 times with
step size 0.1. For the ensemble AutoAttack, we use the standard
$\mathcal{L}_{2}$ norm with maximum perturbation $\epsilon=0.5$ and
$\mathcal{L}_{\infty}$ norm with maximum perturbation $\epsilon=8/255$ as in
the benchmark555https://robustbench.github.io/

##  Appendix C Influence of Attack Strength

We note that in adversarial attacks, perturbations are assumed to be
imperceptible [5], otherwise attack detection techniques can be used [43].
Therefore, most related literature on robust techniques against adversarial
attacks assume small perturbation. However, it is hard to obtain a theoretical
analysis for an exact perturbation bound, which highly depends on specific
datasets and model parameters. We instead demonstrate it experimentally. To
show how much perturbation over input $x$ the neural network can defend
against, we conduct an experiment here using SODEF where we increase the FGSM
attack strength $\epsilon$ at input $x$. Note that in the paper, for CIFAR10
and CIFAR100, the pixel values are normalized by $(x-\mu)/\sigma$ where
$\mu=[0.4914,0.4822,0.4465]$ and $\sigma=[0.2023,0.1994,0.2010]$. The attack
results are shown in Fig. 4. The classification accuracy decreases as the
attack strength $\epsilon$ increases. As expected, when the perturbation is
smaller, SODEF can provide good adversarial robustness.

Fig. 4: Classification accuracy under FGSM attack with different parameters
$\epsilon$.

Since perturbation in the input $x$ is altered by the feature extractor
$h_{\bm{\phi}}$, in another experiment, we also directly perturb the feature
$\mathbf{z}(0)$ instead of the input $x$. The results are shown in Fig. 5(c),
where we can observe how the accuracy changes with the attack strength
$\epsilon$. Without any neural ODE layer, the classification accuracy decays
dramatically and becomes zero when $\epsilon$ is increased to $0.5$. However,
while both ODE net and SODEF have demonstrate defense ability over the FGSM
attack even when $\epsilon\geq 0.5$, SODEF performs much better than ODE,
i.e., the classification accuracy at $\epsilon=0.5$ for SODEF is $59.9\%$
while ODE net achieves only $18.5\%$ accuracy.

(a) Eigenvalues of $\nabla f_{\bm{\theta}}(\mathbf{z}(0))$ of ODE net.

(b) Eigenvalues of $\nabla f_{\bm{\theta}}(\mathbf{z}(0))$ of SODEF.

(c) Classification accuracy under FGSM attack with different parameters
$\epsilon$.

Fig. 5: Eigenvalue visualization for the Jacobian matrix $\nabla
f_{\bm{\theta}}(\mathbf{z}(0))$ from ODE net and SODEF, and classification
accuracy under FGSM attack with different parameters $\epsilon$.

##  Appendix D Eigenvalues Visualization

To further show the effect of the SODEF regularizers in Section 3.2, we
compute the eigenvalues of $\nabla f_{\bm{\theta}}(\mathbf{z}(0))$ from ODE
net (no stability constraint is imposed) and SODEF, as shown in Figs. 5(a) and
5(b). We observe that the eigenvalues of the Jacobian matrix from SODEF are
located in the left plane, i.e., the eigenvalues all have negative real parts.
As a comparison, a large number of ODE net’s eigenvalues are located in the
right plane, as shown in Fig. 5(a).

##  Appendix E Classification Accuracy on CIFAR10 using AutoAttack

In the paper, we study the influence of the SODEF integration time $T$ using
CIFAR10. Here in this section, we further test the influence of the SODEF
integration time $T$ using CIFAR10. We use the model provided by [24], which
has nearly $95\%$ clean accuracy on CIFAR10, as $h_{\phi}$. We use AutoAttack
with $\mathcal{L}_{2}$ norm ($\epsilon=0.5$). We still observe that for all
the four individual attacks and the strongest ensemble AutoAttack, SODEF has a
good performance for large enough integration time $T$.

Table 9: Classification accuracy (%) under AutoAttack on adversarial CIFAR10 examples with $\mathcal{L}_{2}$ norm, $\epsilon=0.5$ and different integration time $T$ for SODEF. Attack / $T$ | 1 | 3 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---
Clean | 95.10 | 94.95 | 95.05 | 95.00 | 94.92 | 95.10 | 95.11 | 95.02
APGD${}_{\text{CE}}$ | 3.52 | 8.20 | 8.59 | 12.59 | 73.05 | 92.41 | 92.58 | 92.67
APGD${}^{\text{T}}_{\text{DLR}}$ | 8.20 | 8.94 | 9.49 | 11.58 | 71.06 | 92.19 | 91.88 | 91.80
FAB${}^{\text{T}}$ | 26.95 | 93.93 | 95.05 | 95.00 | 94.92 | 95.10 | 95.11 | 95.02
Squre | 73.95 | 76.95 | 80.30 | 80.62 | 81.30 | 80.86 | 83.59 | 85.55
AutoAttack | 0.10 | 2.73 | 3.34 | 4.52 | 24.22 | 80.16 | 81.25 | 82.81

##  Appendix F Classification Accuracy on CIFAR100 using FGSM and PGD

We further test SODEF on the CIFAR100 dataset. Table 10 shows that SODEF is
consistenly more robust than all the competitors, especially when stronger
attack is applied.

Table 10: Classification accuracy (%) on adversarial CIFAR100 examples. Attack | Para. | no ode | ODE | TisODE | SODEF
---|---|---|---|---|---
None | - | 88.3 | 88.1 | 88.14 | 88.0
FGSM | $\epsilon=0.1$ | 25.32 | 29.50 | 31.90 | 37.67
PGD | $\epsilon=0.1$ | 2.39 | 3.36 | 6.82 | 22.35

We then do an experiment for the CIFAR100 dataset where we fix the entire
feature extractor (EffcientNet, pretrained on imagenet and CIFAR100, used in
our paper already) except for the last layer in the extractor. The results are
shown as follows in Table 11. In the first line we use the pretrained
EffcientNet with our SODEF and with fine-tuning. In the second line we use a
fixed pretrained EffcientNet and with trainable SODEF. In the third line we
just attack the pretrained EffcientNet (no ODE and no training since this
opensource model has already been trained on CIFAR100 and get $88\%$ clean
accuracy).

We observe the fine-tuned model performs the best among three methods and
fixed (w/ SODEF) performs better than the vanilla model. One possible reason
for fixed (w/ SODEF) being inferior to the fine-tuned model is that the
pretrained feature extractor does not have a well diversified $\mathbf{z}(0)$.

Table 11: Classification accuracy (%) using SODEF under PGD and FGSM attack on adversarial CIFAR100 examples with parameters $\epsilon=0.1$. Model / Attack | clean | FGSM | PGD
---|---|---|---
Fine-tuned (w/ SODEF) | 88.0 | 37.67 | 22.35
Fixed (w/ SODEF) | 87.1 | 33.65 | 17.32
Fixed (no ODE block) | 88.0 | 25.32 | 2.39

##  Appendix G Further Ablation Studies

We performed more experiments on the CIFAR10 dataset to better understand the
proposed SODEF. We let the feature extractor $h_{\phi}$ be ResNet18 (which
ends up with a feature vector $\mathbf{z}(0)\in\mathbb{R}^{10}$) and let
$f_{\theta}$ be a linear function such that
$f_{\theta}(\mathbf{z}(t))=-\mathbf{C}\mathbf{C}^{\intercal}\mathbf{z}(t)$
where $\theta=\mathbf{C}\in\mathbb{R}^{10\times 10}$. Since
$-\mathbf{C}\mathbf{C}^{\intercal}$ is a negative definite matrix for full
rank $\mathbf{C}$, the ODE block is guaranteed to be asymptotically stable
according to Theorem 2. Note that this setup is a special case of our proposed
SODEF architecture, which allows $f_{\theta}$ to be non-linear. The final FC
is fixed to be $\mathbf{V}$ given in Corollary 1. For comparison, we construct
three other neutral network architectures as described in Table 12.

Table 12: Architectures used in our ablation study. Arch. | $h_{\phi}$ | $f_{\theta}$ | FC
---|---|---|---
SODEF | ResNet18 | $-\mathbf{C}\mathbf{C}^{\intercal}\mathbf{z}(t)$ | $\mathbf{V}$
Net1 | ResNet18 | $-\mathbf{C}\mathbf{C}^{\intercal}\mathbf{z}(t)$ | random orth. matrix
Net2 | ResNet18 | $-\mathbf{C}\mathbf{C}^{\intercal}\mathbf{z}(t)$ | linear (trainable)
Net3 | ResNet18 | $\mathbf{C}^{\intercal}\mathbf{z}(t)+{\bf b}$ | $\mathbf{V}$

Table 13 provides the following insights: 1) Comparing SODEF with Net3
highlights the importance of stability of the ODE block; 2) Comparing SODEF
with Net 1 and Net2 shows the superiority of our designed FC layer; 3) By
comparing Net2 with Net1 and SODEF, we see that enforcing $\mathbf{z}(T)$ to
align with a prescribed feature vector, e.g., a row in $\mathbf{V}$, improves
model robustness.

Table 13: Classification accuracy (%) on adversarial CIFAR10 examples. We apply PGD for 20 iterations to generate adversarial examples. Attack | Para. | SODEF | Net1 | Net2 | Net3
---|---|---|---|---|---
None | - | 91.5 | 91.1 | 91.3 | 91.9
PGD | $\epsilon=0.1$ | 29.2 | 24.7 | 5.6 | 18.8

##  Appendix H Further Transferability Studies

We use the same SODEF setup as in Appendix G. In order to determine if SODEF’s
robustness is a result of gradient masking, we generate adversarial examples
from two substitute networks whose architectures are given in Table 14 and
feed the generated examples into SODEF.

Table 14: Two substitute networks used to generate adversarial examples by running PGD attack with $\epsilon=0.1$ for 20 iterations. Both networks have no ODE block, i.e., $f_{\theta}=0$ Arch. | $h_{\phi}$ | $f_{\theta}$ | FC
---|---|---|---
NetA | ResNet18 | - | linear (trainable)
NetB | ResNet18 | - | $\mathbf{V}$

Table 15 provides the following insights: 1) The SODEF column shows that
performing the white-box attack directly on SODEF is more effective than
performing transferred attacks using substitute models such as NetA and NetB
and thereby implying that the ODE block in SODEF does not mask the gradient
when performing PGD attack. 2) The NetA and NetB columns show that even when
there is no ODE block, NetB using $\mathbf{V}$ as the last FC layer is much
more robust than NetA whose FC layer is a trainable dense layer.

Table 15: Classification accuracy (%) on adversarial CIFAR10 examples. We generate adversarial examples by attacking NetA and NetB using PGD method (run for 20 iterations with $\epsilon=0.1$) and then feed them to SODEF. Adv. examples gen. by: | SODEF | NetA | NetB
---|---|---|---
No Attack | 91.5 | 93.6 | 91.5
SODEF | 29.2 | - | -
NetA | 43.9 | 2.6 | -
NetB | 51.3 | - | 20.5

##  Appendix I Different ODE Solvers

We test different ODE solvers including Runge-Kutta of order 5, Runge-Kutta of
order 2, Euler method, Midpoint method, and Fourth-order Runge-Kutta with 3/8
rule. All the ODE solvers tested in Table 16 show similar performance. For
Euler method, Midpoint method, and Fourth-order Runge-Kutta with 3/8 rule, we
set step size to $0.05$. See https://github.com/rtqichen/torchdiffeq.

Table 16: Classification accuracy (%) using SODEF under FGSM attack on adversarial CIFAR10 examples with parameters $\epsilon=0.1$ and different ODE solvers. Model / Attack | clean | FGSM
---|---|---
Runge-Kutta of order 5 | 95.00 | 68.05
Runge-Kutta of order 2 | 95.12 | 69.95
Euler method | 95.20 | 70.22
Midpoint method | 95.10 | 70.86
Fourth-order Runge-Kutta with 3/8 rule | 95.07 | 69.50

##  Appendix J Proofs of Results in Paper

In this section, we provide detailed proofs of the results stated in the main
text.

###  J.1 Proof of Lemma 2

For simplicity, let $a=a(\mathbf{v}_{1},\ldots,\mathbf{v}_{k})$. Since
$0\leq\|\sum_{i=1}^{k}\mathbf{v}_{i}\|^{2}=\sum_{i}\|\mathbf{v}_{i}\|^{2}+2\sum_{i\leq
j}\mathbf{v}_{i}^{\intercal}\mathbf{v}_{j}\leq k+2\sum_{i\leq j}a=k+k(k-1)a$,
we have $a\geq 1/(1-k)$. We next prove $a=1/(1-k)$ is attainable. We can
construct a $k\times k$ matrix $\mathbf{B}$ with
$\mathbf{B}(i,i)=1,\forall\,i$ and $\mathbf{B}(i,j)=1/(1-k),\ \forall\,i\neq
j$. It can be verified that $\mathbf{B}$ is a diagonally dominant matrix. Due
to the fact that a symmetric diagonally dominant real matrix with nonnegative
diagonal entries is positive semi-definite and hence a Gram matrix w.r.t. the
Euclidean inner product [19], there exists a set of vectors $\mathbf{v}_{i}$
such that $\mathbf{v}_{i}^{\intercal}\mathbf{v}_{j}=\mathbf{B}(i,j)=1/(1-k)$
for $i\neq j$. In this case, $a(\mathbf{v}_{1},\ldots,\mathbf{v}_{k})=1/(1-k)$
and the proof is complete.

###  J.2 Proof of Corollary 1

We may decompose $\mathbf{B}$ as

$\displaystyle\mathbf{B}$
$\displaystyle=(\bm{\Sigma}^{1/2}\mathbf{U}^{\intercal})^{\intercal}(\bm{\Sigma}^{1/2}\mathbf{U}^{\intercal})$
$\displaystyle=\left(\mathbf{Q}\bm{\Sigma}^{1/2}\mathbf{U}^{\intercal}\right)^{\intercal}\left(\mathbf{Q}\bm{\Sigma}^{1/2}\mathbf{U}^{\intercal}\right),$

and hence $\mathbf{B}$ is the Gram matrix of
$\\{\mathbf{v}_{1},\ldots,\mathbf{v}_{k}\\}$. The result follows from the
construction of $\mathbf{B}$.

###  J.3 Proof of Lemma 1

The set of finite points $\\{\mathbf{z}_{1},...,\mathbf{z}_{k}\\}$ is closed
and this lemma is an immediate consequence of the Whitney extension theorem
[45].

###  J.4 Proof of Lemma 3

If a continuous function $f$ is such that $f(\mathbf{z})=0$
$\nu_{\bm{\phi}}$-almost surely, then $f(\mathbf{z})=0$ for all $\mathbf{z}\in
E^{\prime}$. Suppose otherwise, then there exists $\mathbf{z}_{0}$ such that
$\lVert f(\mathbf{z}_{0})\rVert_{2}>0$. From continuity, there exists an open
ball $B\in E^{\prime}$ containing $\mathbf{z}_{0}$ such that $\lVert
f(\mathbf{z})\rVert_{2}>0$ for all $\mathbf{z}\in B$. Then since $B$ has non-
zero Lebesgue measure, $\nu_{\bm{\phi}}(B)>0$, a contradiction. Therefore,
every $\mathbf{z}\in E^{\prime}$ has $\nabla f(\mathbf{z})=0$ and cannot be a
Lyapunov-stable equilibrium point.

###  J.5 Proof of Theorem 4

Consider
$f(\mathbf{z})=[f^{(1)}(\mathbf{z}^{(1)}),...,f^{(n)}(\mathbf{z}^{(n)})]$ with
each $f^{(i)}(\mathbf{z}^{(i)})\in C^{1}(\mathbb{R},\mathbb{R})$. Since
$f^{(i)}(\mathbf{z}^{(i)})$ only depends on $\mathbf{z}^{(i)}$, $\nabla
f_{\bm{\theta}}(\mathbf{z})$ is a diagonal matrix with all off-diagonal
elements being $0$. The constraint 4 is thus satisfied immediately and it
suffices to show that there exists such a $f$ satisfying the constraints 3 and
LABEL:eq:re2.

Select a $\mathbf{z}_{l}=(\mathbf{z}^{(1)}_{l},\ldots,\mathbf{z}^{(n)}_{l})$
from the interior of each $E_{l}$, $l=1,\dots,L$. Let
$f^{(i)}(\mathbf{z}^{(i)})=-\beta(\mathbf{z}^{(i)}-\mathbf{z}_{l}^{(i)})$ on
each $E_{l}$, where $\beta>0$. Then $f(\mathbf{z})$ satisfies LABEL:eq:re2 for
all $\beta>0$ and $\mathbf{z}_{l}$ is a Lyapunov-stable equilibrium point for
each $l$ since $\nabla f_{\bm{\theta}}(\mathbf{z}_{l})$ is a diagonal matrix
with negative diagonal values. Since each $E_{l}\subset\mathbb{R}^{n}$ is
compact, we have that $\forall\,\epsilon>0$, $\exists\>\beta>0$ sufficiently
small such that $|f^{(i)}(\mathbf{z}^{(i)})|<\epsilon$ for all
$\mathbf{z}\in\bigcup_{l}E_{l}$. The constraint 3 is therefore satisfied for
$f(\mathbf{z})$ with a sufficiently small $\beta$. Since $\bigcup_{l}E_{l}$ is
closed, the Whitney extension theorem [45] can be applied to extend
$f(\mathbf{z})$ to a function in $C^{1}(\mathbb{R}^{n},\mathbb{R}^{n})$.

## References

  * [1] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “ImageNet classification with deep convolutional neural networks,” in _Proc. Advances Neural Inf. Process. Syst._ , 2012.
  * [2] Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner, “Gradient-based learning applied to document recognition,” _Proceedings of the IEEE_ , vol. 86, no. 11, pp. 2278–2324, Nov. 1998.
  * [3] G. Hinton, L. Deng, D. Yu, G. E. Dahl, A. Mohamed, N. Jaitly, A. Senior, V. Vanhoucke, P. Nguyen, T. N. Sainath, and B. Kingsbury, “Deep neural networks for acoustic modeling in speech recognition: The shared views of four research groups,” _IEEE Signal Process. Mag._ , vol. 29, no. 6, pp. 82–97, Nov. 2012.
  * [4] D. Andor, C. Alberti, D. Weiss, A. Severyn, A. Presta, K. Ganchev, S. Petrov, and M. Collins, “Globally normalized transition-based neural networks,” in _Proc. Annu. Meeting Assoc, Comput. Linguistics_ , 2016.
  * [5] C. Szegedy, W. Zaremba, I. Sutskever, J. Bruna, D. Erhan, I. Goodfellow, and R. Fergus, “Intriguing properties of neural networks,” in _Proc. Int. Conf. Learning Representations_ , 2013.
  * [6] H. Yan, J. Du, V. Y. Tan, and J. Feng, “On robustness of neural ordinary differential equations,” in _Proc. Advances Neural Inf. Process. Syst._ , 2018, pp. 1–13.
  * [7] E. Haber and L. Ruthotto, “Stable architectures for deep neural networks,” _Inverse Problems_ , vol. 34, no. 1, pp. 1–23, Dec. 2017.
  * [8] X. Liu, S. Si, Q. Cao, S. Kumar, and C.-J. Hsieh, “How does noise help robustness? Explanation and exploration under the neural sde framework,” in _Proc. Conf. Comput. Vision Pattern Recognition_ , 2020, pp. 282–290.
  * [9] R. T. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud, “Neural ordinary differential equations,” _arXiv preprint arXiv:1806.07366_ , 2018.
  * [10] C.-T. Chen and B. Shafai, _Linear system theory and design_. New York: Oxford university press New York, 1999\.
  * [11] Y. LeCun, C. Corte, and C. Burges, “MNIST handwritten digit database,” _ATT Labs [Online]. Available: http://yann.lecun.com/exdb/mnist_ , vol. 2, 2010, (Last accessed: Dec 1, 2020).
  * [12] A. Krizhevsky and G. Hinton, “Learning multiple layers of features from tiny images,” _Master’s thesis, Department of Computer Science, University of Toronto_ , 2009.
  * [13] A. Mądry, A. Makelov, L. Schmidt, D. Tsipras, and A. Vladu, “Towards deep learning models resistant to adversarial attacks,” in _Proc. Int. Conf. Learning Representations_ , 2018.
  * [14] E. Haber and L. Ruthotto, “Stable architectures for deep neural networks,” _Inverse Problems_ , vol. 34, no. 1, 2018.
  * [15] Y. Huang, Y. Yu, H. Zhang, Y. Ma, and Y. Yao, “Adversarial robustness of stabilized neuralodes might be from obfuscated gradients,” _arXiv preprint arXiv:2009.13145_ , 2020.
  * [16] L. Mingjie, H. Lingshen, and L. Zhouchen, “Implicit Euler skip connections: Enhancing adversarial robustness via numerical stability,” in _Proc. Int. Conf. Machine Learning_ , 2020.
  * [17] I. J. Goodfellow, J. Shlens, and C. Szegedy, “Explaining and harnessing adversarial examples,” in _Proc. Int. Conf. Learning Representations_ , 2015\.
  * [18] D. Arrowsmith and C. M. Place, _Dynamical systems: differential equations, maps, and chaotic behaviour_. London: CRC Press, 1992.
  * [19] R. A. Horn and C. R. Johnson, _Matrix analysis_. New York: Cambridge university press, 2012.
  * [20] L. Van der Maaten and G. Hinton, “Visualizing data using t-sne,” _J. Mach. Learning Res._ , vol. 9, no. 11, pp. 2580–2605, Nov. 2008.
  * [21] K. Hornik, “Approximation capabilities of multilayer feedforward networks,” _Neural Networks_ , vol. 4, no. 2, pp. 251–257, Oct. 1991.
  * [22] Y. Luo, Y. Wong, S. M. Kankanhalli, and Q. Zhao, “Direction concentration learning: Enhancing congruency in machine learning,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 43, no. 1, pp. 1928 – 1946, Jun. 2021.
  * [23] M. Tan and Q. Le, “Efficientnet: Rethinking model scaling for convolutional neural networks,” in _Proc. Int. Conf. Mach. Learning_ , 2019, pp. 6105–6114.
  * [24] Github pytorch-cifar repository. Accessed: May 1, 2021. [Online]. Available: https://github.com/kuangliu/pytorch-cifar
  * [25] T. Pang, X. Yang, Y. Dong, H. Su, and J. Zhu, “Bag of tricks for adversarial training,” in _Proc. Int. Conf. Learning Representations_ , 2021.
  * [26] F. Croce and M. Hein, “Reliable evaluation of adversarial robustness with an ensemble of diverse parameter-free attacks,” in _Proc. Int. Conf. Mach. Learning_ , 2020, pp. 2206–2216.
  * [27] ——, “Minimally distorted adversarial examples with a fast adaptive boundary attack,” in _Proc. Int. Conf. Mach. Learning_ , 2020, pp. 2196–2205.
  * [28] M. Andriushchenko, F. Croce, N. Flammarion, and M. Hein, “Square attack: a query-efficient black-box adversarial attack via random search,” in _Proc. European Conf. Comput. Vision_. Springer, 2020, pp. 484–501.
  * [29] Y. Bengio, P. Simard, and P. Frasconi, “Learning long-term dependencies with gradient descent is difficult,” _IEEE Tran. Neural Networks_ , vol. 5, no. 2, pp. 157–166, 1994.
  * [30] M. Sharif, S. Bhagavatula, L. Bauer, and M. K. Reiter, “Accessorize to a crime: Real and stealthy attacks on state-of-the-art face recognition,” in _Proc. ACM SIGSAC Conf. Comput. Commun._ , 2016, pp. 1528–1540.
  * [31] S.-M. Moosavi-Dezfooli, A. Fawzi, and P. Frossard, “Deepfool: A simple and accurate method to fool deep neural networks,” in _Proc. Conf. Comput. Vision Pattern Recognition_ , 2016, pp. 2574–2582.
  * [32] N. Carlini and D. Wagner, “Audio adversarial examples: Targeted attacks on speech-to-text,” in _Proc. IEEE Security and Privacy Workshops_ , 2018, pp. 1–7.
  * [33] L. Schönherr, K. Kohls, S. Zeiler, T. Holz, and D. Kolossa, “Adversarial attacks against automatic speech recognition systems via psychoacoustic hiding,” _arXiv preprint arXiv:1808.05665_ , 2018.
  * [34] C. Xiang, C. R. Qi, and B. Li, “Generating 3d adversarial point clouds,” in _Proc. Conf. Comput. Vision Pattern Recognition_ , 2019, pp. 9136–9144.
  * [35] Y. Cao, C. Xiao, B. Cyr, Y. Zhou, W. Park, S. Rampazzi, Q. A. Chen, K. Fu, and Z. M. Mao, “Adversarial sensor attack on lidar-based perception in autonomous driving,” in _Proc. ACM SIGSAC Conf. Comput. Commun._ , 2019, pp. 2267–2281.
  * [36] A. Kurakin, I. J. Goodfellow, and S. Bengio, “Adversarial examples in the physical world,” in _Proc. Int. Conf. Learning Representations_ , 2017.
  * [37] N. Papernot, P. McDaniel, S. Jha, M. Fredrikson, Z. B. Celik, and A. Swami, “The limitations of deep learning in adversarial settings,” in _Proc. IEEE Eur. Symp. Security Privacy_ , 2016.
  * [38] N. Carlini and D. Wagner, “Towards evaluating the robustness of neural networks,” in _Proc. IEEE Symp. Security and Privacy_ , 2017.
  * [39] N. Papernot, P. McDaniel, X. Wu, S. Jha, and A. Swami, “Distillation as a defense to adversarial perturbations against deep neural networks,” in _Proc. IEEE Symp. Security and Privacy_ , 2016.
  * [40] G. Verma and A. Swami, “Error correcting output codes improve probability estimation and adversarial robustness of deep neural networks,” in _Proc. Advances Neural Inf. Process. Syst._ , 2019.
  * [41] Y. Song, Q. Kang, and W. P. Tay, “Error-correcting output codes with ensemble diversity for robust learning in neural networks,” in _Proc. AAAI Conf. Artificial Intelligence_ , Feb. 2021.
  * [42] D. Hendrycks and K. Gimpel, “A baseline for detecting misclassified and out-of-distribution examples in neural networks,” in _Proc. Int. Conf. Learning Representations_ , 2017.
  * [43] D. Meng and H. Chen, “Magnet: a two-pronged defense against adversarial examples,” in _Proc. SIGSAC Conf. Comput. Commun. Security_ , 2017.
  * [44] P. Samangouei, M. Kabkab, and R. Chellappa, “Defense-GAN: protecting classifiers against adversarial attacks using generative models,” in _Proc. Int. Conf. Learning Representations_ , 2018.
  * [45] E. J. McShane, “Extension of range of functions,” _Bulletin Amer. Math. Soc._ , vol. 40, no. 12, pp. 837–842, Jun. 1934.