## Transversal transport of magnons in a modified Lieb lattice

Pedro Gonçalves de Oliveiraa (P. G. de Oliveira)
E-mail address<EMAIL_ADDRESS>(corresponding author)

and

Antônio Sérgio Teixeira Piresa (A. S. T. Pires)
E-mail address<EMAIL_ADDRESS>
aDepartamento de Física, Universidade Federal de Minas Gerais, Belo Horizonte,
MG, CP702, 30123-970, Brazil.

### Abstract

We studied a two-band magnon insulating model whose geometry is that of a
modified Lieb lattice, on which one of the sites was removed. There is an
anisotropic ferromagnetic exchange interaction between the three nearest
neighbours, and the anisotropy opens a gap in the magnon energy band
structure. A non-vanishing Berry curvature is induced by a Dzyaloshinskii-
Moriya interaction (DMI). The topology of the bands is trivial (in the sense
of a null Chern number), but the finite Berry curvature induces Hall-like
transport effects whose coefficients were calculated. Their dependence on
temperature and applied magnetic field were studied and show resemblance with
other magnon insulating systems found on the literature.

Keywords: Spin waves; magnons; Dzyaloshinskii–Moriya interaction; transport;
Hall-like effects.

##  1 Introduction

Topological effects on condensed matter systems have been intensely studied
since the discovery of the quantum Hall effect by von Klitzing et al [1]. The
attention falls naturally on the so-called topological insulators (TIs), which
are electronic systems that have gapped bands on the bulk and robust
conducting edge (or surface) modes. These systems have different insulating
phases characterized by topological indices [2], and may show Hall-like
effects when subjected to a field or temperature gradient [3, 4, 5, 6, 7].
These effects arise in materials with strong spin-orbit coupling, and can be
related to the Berry phase and Berry curvature of the electronic bands [8].

In analogy to TIs in electronic systems, there has been a great interest in
topological magnon insulators (TMIs) lately. Magnons are spin-wave excitations
of the ground state of systems of localized spins, and when magnon bands have
non-trivial topology (in the form of finite Berry curvature) the same Hall-
like transport effects can arise [9, 10, 11, 12, 13, 14, 15, 16, 17]. Since
magnons are bosons, magnonic systems are intrinsically different from
electronic ones, which motivates their study. A notable fact is that, because
of their uncharged nature, magnons favour dissipationless transport, being of
great interest to spintronics [18]. While topological effects in magnon
systems were first discovered in a three-dimensional material with the
geometry of the pyrochlore lattice [19], the main theoretical interest
nowadays falls on two-dimensional lattices, where the most studied geometries
are the honeycomb lattice [20, 21, 22] and the kagome lattice [12, 23, 24, 25,
14]. The latter can be layered with triangular planes to form the pyrochlore
structure. Other lattices that were predicted to hold topological magnon
effects are the Shastry-Sutherland [26], square [27], checkerboard [28, 29,
30] and Lieb [31] lattices.

The latter is of particular interest because it’s the geometry that $CuO_{2}$
planes assume in high-$T_{c}$ cuprate superconductors [32]. In a tight-binding
approach it is a three-band model that shows a flat band [33] and a single
Dirac cone in the Brillouin zone. The energy gap can be opened with the
creation of a TI phase by an intrinsic spin-orbit interaction term [34, 35].
The Berry curvature and anomalous Hall effect of the electronic Lieb lattice
were studied by He et al [36]. A magnon model for the same lattice was studied
by Cao et al [31], where it was shown that topological insulating phases can
be induced by a complex hopping between next-near-neighbours. The insulating
system presents thermal spin Hall effect, and the temperature dependence of
the thermal Hall conductivity is related to its topological phases.

Based on the necessity of finding novel magnetic systems with non-trivial
topologies, we propose a study of a modified version of the ferromagnetic Lieb
lattice, where one of the inequivalent sites is removed. The topological
effects are induced by a Dzyaloshinskii-Moriya interaction (DMI) between the
next-next-near neighbours, which is equivalent to the above mentioned complex
hopping. The DMI breaks the time-reversal symmetry (TRS), and it is the most
common way of inducing non-vanishing Berry curvature and topological effects.
We present the system’s geometry and Hamiltonian in section 2, calculate its
magnon band structure in section 3, calculate and discuss its Hall-like
transport coefficients in section 4 and present our conclusions in section 5.

##  2 Model

We consider a lattice with two inequivalent sites (A and B) in each square
unit cell, with the following Hamiltonian:

$\displaystyle H=$ $\displaystyle-J_{1}\sum_{\left\langle
i,j\right\rangle}\left(S_{i}^{x}S_{j}^{x}+S_{i}^{y}S_{j}^{y}+\lambda
S_{i}^{z}S_{j}^{z}\right)-J_{2}\sum_{\left\langle\left\langle
i,j\right\rangle\right\rangle\in
A}\left(S_{i}^{x}S_{j}^{x}+S_{i}^{y}S_{j}^{y}+\lambda
S_{i}^{z}S_{j}^{z}\right)$ $\displaystyle-
J_{3}\sum_{\left\langle\left\langle\left\langle
i,j\right\rangle\right\rangle\right\rangle}\left(S_{i}^{x}S_{j}^{x}+S_{i}^{y}S_{j}^{y}+\lambda
S_{i}^{z}S_{j}^{z}\right)$
$\displaystyle-D\sum_{\left\langle\left\langle\left\langle
i,j\right\rangle\right\rangle\right\rangle}\nu_{ij}\left(S_{i}^{x}S_{j}^{y}-S_{i}^{y}S_{j}^{x}\right)-B\sum_{i}S_{i}^{z}$
(1)

The lattice can be seen as a modified Lieb lattice in which one of the three
inequivalent sites was removed. The unit cell is a square, which side was
taken as equal to one. A sketch of the lattice can be seen in Figure 1. There
are ferromagnetic exchange interactions (solid lines) between near-neighbours
$A$ and $B$ (strength $J_{1}$) and next-near-neighbours $A$ (strength
$J_{2}$). On the diagonals between $A$ and $B$ (dashed lines) there are a
$J_{3}$ exchange interaction and a Dzyaloshinskii-Moriya interaction (DMI)
[37, 38] of the form
$-\mathbf{D}_{ij}\cdot\left(\mathbf{S}_{i}\times\mathbf{S}_{j}\right)$. The
later is responsible for the finite Berry curvature and Hall-like transport
effects which are the focus of this paper.

Figure 1: Modified Lieb lattice studied on this paper. Circles (squares)
represent the A (B) sites. Solid lines are ferromagnetic exchange bonds
($J_{1}$ between A and B; $J_{2}$ between two A sites). Dashed lines represent
both $J_{3}$ and DM interaction. The unit cell is represented by the dashed
square.

The DM interaction is not forbidden by Moriya’s rules, for there is no center
of inversion in the midpoint of the bond [38]. Hence we can consider the
interaction as a regular DMI and there is no need for introducing and external
electric field to induce the interaction, as it is the case in some lattices
where the regular DMI is forbidden [39]. The same Moriya’s rules restrict the
DM vector to the $z$ direction, as a two-dimensional lattice is symmetric in
respect to a reflection upon its own plane. We take
$\mathbf{D}_{ij}=D\nu_{ij}\mathbf{\hat{z}}$ (fourth term in (1)), where,
$\nu_{ij}=\pm 1$ for different bond directions, following Figure 2. In all the
exchange interactions there is an anisotropy $\lambda>1$ on the $z$ direction,
which is responsible for the stabilization of the magnetic order in an easy-
axis configuration [40, 41] (otherwise the system would have no long-range
order according to the Mermim-Wagner theorem). The last term in (1) is a
Zeeman interaction with a constant magnetic field
$\mathbf{B=}B\mathbf{\hat{z}.}$

Figure 2: Configuration of the DMI vectors D on the diagonal bonds, departing
from an $A$ site.

##  3 Magnon bands

We are interested in the linear spin-wave regime, so we use the zeroth order
expansion of the Holstein-Primakoff representation:

$\displaystyle S_{i}^{+}$ $\displaystyle\approx\sqrt{2S}\ a_{i}\text{ ;
}S_{i}^{-}\approx\sqrt{2S}\ a_{i}^{\dagger}\text{ ;
}S_{i}^{z}=S-a_{i}^{\dagger}a_{i}$ (2) $\displaystyle S_{j}^{+}$
$\displaystyle\approx\sqrt{2S}\ b_{j}\text{ ; }S_{j}^{-}\approx\sqrt{2S}\
b_{j}^{\dagger}\text{ ; }S_{j}^{z}=S-b_{j}^{\dagger}b_{j}$

We are studying a ferromagnetic model, but since we have two inequivalent
sites we need two operators $a$ and $b$. Here, the index $i$ represents sites
in the sublattice $A$ and $j$ in sublattice $B$. This spin-wave approximation
works for low temperature, and the magnon-magnon interactions (higher order
therms on the expansion) are much smaller than the linear contribution and can
be neglected [42, 43].

After transforming, keeping only the quadratic terms and applying a Fourier
transform, we can write the harmonic momentum-space Hamiltonian as:

$H_{harm}=S\sum_{k}\psi_{k}^{\dagger}\left(h_{0}\hat{1}_{2}+\hat{M}_{k}\right)\psi_{k}$
(3)

where $\psi_{k}^{\dagger}=\left(a_{k}^{\dagger}\text{ \ \
}b_{k}^{\dagger}\right).$ Here, $h_{0}$ and the $2\times 2$ matrix
$\hat{M}_{k}$ are

$\displaystyle h_{0}$ $\displaystyle=J_{2}\left(\lambda-\cos
k_{x}\right)+2\lambda\left(J_{1}+2J_{3}\right)+\frac{B}{S}$ (4) $\displaystyle
M_{11}$ $\displaystyle=J_{2}\left(\lambda-\cos k_{x}\right)$ $\displaystyle
M_{12}$ $\displaystyle=M_{21}^{\ast}=-2J_{1}\cos\frac{k_{y}}{2}-4J_{3}\cos
k_{x}\cos\frac{k_{y}}{2}-4iDm_{k}$ $\displaystyle M_{22}$
$\displaystyle=-J_{2}\left(\lambda-\cos k_{x}\right)$

where $m_{k}=-\sin k_{x}\sin\frac{k_{y}}{2}$.

We can write the matrix $\hat{M}_{k}$ as a Pauli vector

$\hat{M}=h_{x}\left(\mathbf{k}\right)\hat{\sigma}_{x}+h_{y}\left(\mathbf{k}\right)\hat{\sigma}_{y}+h_{x}\left(\mathbf{k}\right)\hat{\sigma}_{z}$
(5)

so that the dispersion relation is simply

$\frac{E_{\pm}}{\hbar}=\omega_{\pm}\left(\mathbf{k}\right)=S\left(h_{0}\left(\mathbf{k}\right)\pm
h\left(\mathbf{k}\right)\right)$ (6)

where
$h\left(\mathbf{k}\right)=\left\|\mathbf{h}\left(\mathbf{k}\right)\right\|=\sqrt{h_{x}^{2}\left(\mathbf{k}\right)+h_{y}^{2}\left(\mathbf{k}\right)+h_{z}^{2}\left(\mathbf{k}\right)},$
and

$\displaystyle h_{x}\left(\mathbf{k}\right)$
$\displaystyle=-2\cos\frac{k_{y}}{2}\left(J_{1}+2J_{3}\cos k_{x}\right)$ (7)
$\displaystyle h_{y}\left(\mathbf{k}\right)$ $\displaystyle=4Dm_{k}$
$\displaystyle h_{z}\left(\mathbf{k}\right)$
$\displaystyle=J_{2}\left(\lambda-\cos k_{x}\right)$

Sometimes it is useful to perform the substitution $\tan\phi=\frac{D}{J_{3}},$
arriving at

$h_{x}\left(\mathbf{k}\right)=-2J_{1}\cos\frac{k_{y}}{2}-4J_{D}\gamma_{k}\cos\phi\text{
\ , \ }h_{y}\left(\mathbf{k}\right)=4J_{D}m_{k}\sin\phi$ (8)

where $\gamma_{k}=\cos k_{x}\cos\frac{k_{y}}{2}$ and
$J_{D}=\sqrt{J_{3}^{2}+D^{2}}$. The phase $\phi$ can be seen as a magnetic
flux generated by the DM term. In the case of a pure DM interaction in the
diagonals (no exchange $J_{3}$), the phase is $\phi=\pi/2$.

The explicit dispersion relation is

$\displaystyle\omega_{\pm}\left(\mathbf{k}\right)$
$\displaystyle=SJ_{2}\left(\lambda-\cos
k_{x}\right)+2S\lambda\left(J_{1}+2J_{3}\right)+B+$ (9) $\displaystyle\pm
S\sqrt{4\cos^{2}\frac{k_{y}}{2}\left(J_{1}+2J_{3}\cos
k_{x}\right)^{2}+J_{2}^{2}\left(\lambda-\cos
k_{x}\right)^{2}+16D^{2}m_{k}^{2}}$

The band structure of the system is plotted in Figure 3 for $S=1/2$, $B=0$,
$J_{1}=1$, $J_{2}=0.5$, $J_{3}=0.2$, $D=0.1$ and $\lambda=1.5$. There is a gap
of width $2SJ_{2}\left(\lambda-1\right)$ at the high-symmetry point
$X^{\prime}=\left(0,\pm\pi\right)$. The gap vanishes into a Dirac point in the
isotropic limit ($\lambda=1$) independently of the value of other parameters.
Therefore, the anisotropy is responsible for the gap, and not the DM
interaction, as it is the case of other magnon insulating systems [20, 22, 31,
30]. For a zero applied field $B=0$ and small enough values of $D$ the valence
band has a single global minimum at the $\Gamma$ point, which has zero value
if and only if $\lambda=1$. Hence, in the isotropic limit we have a gapless
system with a Goldstone mode. For our purposes this possibility is absent,
since we defined $\lambda>1$ for an easy axis configuration. Therefore, we
have an insulating system with no Goldstone mode.

Figure 3: Band structure of our system for $S=1/2$, $B=0$, $J_{1}=1$,
$J_{2}=0.5$, $J_{3}=0.2$, $D=0.1$ and $\lambda=1.5$.

##  4 Transport coefficents

For magnon systems it is known that an external in-plane magnetic field
gradient generates not only a parallel spin current, but also a spin response
on the transverse direction [9]:

$j_{x}^{S}=-\sigma_{xy}\left(\partial_{y}B\right)$ (10)

This is the so-called spin Hall effect of magnons, and it is responsible for a
protected spin current on the edges of a 2D magnet [16, 17]. In a
semiclassical picture, this current can be explained by the effect of the
borders, which act as an effective field gradient, confining the magnons
inside the magnet and generating a spin current perpendicular to the gradient
(along the edge)$\cdot$

The transverse conductivity $\sigma_{xy}$ can be obtained from the Berry
curvature of the system [10, 11]:

$\sigma_{xy}=-\frac{1}{\hbar V_{BZ}}\sum_{\lambda}\int_{BZ}dk_{x}dk_{y}\text{
}n_{\lambda}\left(\mathbf{k}\right)\Omega_{\lambda}\left(\mathbf{k}\right)$
(11)

where $\lambda$ sweeps the magnon energy bands. Here,
$n_{\lambda}\left(\mathbf{k}\right)=\left(e^{\hbar\omega_{\lambda}\left(\mathbf{k}\right)/k_{B}T}-1\right)^{-1}$
is the Bose-Einstein distribution and
$\Omega_{\lambda}\left(\mathbf{k}\right)$ is the off-plane component of the
(vector) Berry curvature of the band defined as [10, 11, 8]:

$\mathbf{\Omega_{\lambda}}\left(\mathbf{k}\right)=i\left\langle\nabla_{\mathbf{k}}u_{\lambda}\left(\mathbf{k}\right)\right|\times\left|\nabla_{\mathbf{k}}u_{\lambda}\left(\mathbf{k}\right)\right\rangle$
(12)

where $\left|u_{\lambda}\left(\mathbf{k}\right)\right\rangle$ is the Bloch
wave function of the $\lambda$ band. For a two-band system parametrized by
$\mathbf{h}\left(\mathbf{k}\right)=\left(h_{x}\left(\mathbf{k}\right),h_{y}\left(\mathbf{k}\right),h_{z}\left(\mathbf{k}\right)\right)$
as ours, the Berry curvatures can be written as [22, 8]

$\Omega_{+}\left(\mathbf{k}\right)=-\frac{1}{2h^{3}}\mathbf{h}\cdot\left(\partial_{k_{x}}\mathbf{h}\times\partial_{k_{y}}\mathbf{h}\right)$
(13)

and $\Omega_{-}\left(\mathbf{k}\right)=-\Omega_{+}\left(\mathbf{k}\right)$.
For our system the definition above gives:

$\displaystyle\Omega_{+}\left(\mathbf{k}\right)=$
$\displaystyle-2\frac{DJ_{2}}{h\left(\mathbf{k}\right)^{3}}\Bigg{\\{}J_{1}\left[\sin^{2}k_{x}-\left(\lambda-\cos
k_{x}\right)\cos k_{x}\sin^{2}\frac{k_{y}}{2}\right]+$ (14)
$\displaystyle+2J_{3}\left[\left(\lambda-\cos
k_{x}\right)\left(\sin^{2}k_{x}\cos^{2}\frac{k_{y}}{2}-\cos^{2}k_{x}\sin^{2}\frac{k_{y}}{2}\right)+\sin^{2}k_{x}\cos
k_{x}\right]\Bigg{\\}}$

The Berry curvature of the valence band through the Brillouin zone is plotted
in Figure 4. As we could expect, the curvature is highly concentrated around
the point in the Brillouin zone where the energy gap occurs [8].

Figure 4: Berry curvature of the valence band of the system for $S=1/2$,
$J_{1}=1$, $J_{2}=0.5$, $J_{3}=0.2$, $D=0.1$ and $\lambda=1.5$.

The Chern number is defined as proportional to the integral of the Berry
curvature in the Brillouin zone and is an integer number which labels the
inequivalent topological phases in Chern insulators [2]. Despite of the finite
Berry curvature, the Chern numbers of the bands are null for any combination
of parameters values (provided that $\lambda>0$), which means that the system
doesn’t present a non-trivial topological insulating phase. Nevertheless, the
non-vanishing Berry curvature gives rise to Hall-like effects like the spin
Hall effect shown above.

The Berry curvature is related to transverse spin and heat currents in
response to an applied temperature gradient [15]. These are the spin Nernst
effect [44],

$j_{x}^{N}=-\alpha_{xy}\left(\partial_{y}T\right)$ (15)

and the thermal Hall effect [10, 11, 12]

$j_{x}^{Q}=-\kappa_{xy}\left(\partial_{y}T\right).$ (16)

Here, $j_{x}^{N}$ is the spin current and $j_{x}^{Q}\,$ the heat current. The
transport coefficients are defined as [10, 11, 44, 21]:

$\displaystyle\alpha_{xy}$ $\displaystyle=-\frac{k_{B}}{\hbar
V_{BZ}}\sum_{\lambda}\int_{BZ}dk_{x}dk_{y}\text{ \
}c_{1}\left(n_{\lambda}\left(\mathbf{k}\right)\right)\Omega_{\lambda}\left(\mathbf{k}\right)$
$\displaystyle\kappa_{xy}$ $\displaystyle=-\frac{k_{B}^{2}T}{\hbar
V_{BZ}}\sum_{\lambda}\int_{BZ}dk_{x}dk_{y}\text{ \
}c_{2}\left(n_{\lambda}\left(\mathbf{k}\right)\right)\Omega_{\lambda}\left(\mathbf{k}\right)$
(17)

with

$\displaystyle c_{1}\left(x\right)$
$\displaystyle=\left(1+x\right)\ln\left(1+x\right)-x\ln x$ $\displaystyle
c_{2}\left(x\right)$
$\displaystyle=\left(1+x\right)\ln\left(\frac{1+x}{x}\right)^{2}-\left(\ln
x\right)^{2}-2Li_{2}\left(-x\right)$ (18)

where $Li_{2}\left(x\right)$ is the Spence’s dilogarithm function.

The spin Hall conductivity $\sigma_{xy}$ is plotted as a function of
temperature for three values of $D$ in Figure 5 (in all transport coefficients
plots the temperature and applied magnetic field are shown in units of the
exchange energy $J_{1}$). The plot shows a monotonically rising behavior of
$\sigma_{xy}$, similar to what could be observed in the checkerboard lattice
[30, 29]. At zero temperature $\sigma_{xy}$ is zero due to the absence of
magnon excitations. That is a consequence of the fact that boson numbers are
not conserved and vanish in the zero temperature limit. However, as the
temperature increases from zero, magnons are thermally excited and
$\sigma_{xy}$ becomes finite. At low temperatures, the lower band dominates.

In all temperature plots the range of the temperature axis was chosen to best
show the character of the curve, but we wouldn’t expect the model to work at
such high temperatures. In fact, the linear spin-wave approximation deals with
perturbations of the ordered ground state and only works at low temperatures.
On way to quantify the validity of the model is calculating the expected value
of the total boson number:

$\Delta=\frac{1}{V_{BZ}}\sum_{\lambda}\int_{BZ}dk_{x}dk_{y}\ n_{\lambda}$ (19)

The spin-wave approach works for $\Delta\ll S$. Some values of $\Delta$ for
different spins and temperatures are represented in Table 1. As we can see,
the approximation works better for low temperatures and high values of spin.

Figure 5: Spin Hall conductivity versus temperature for $S=1/2$, $B=0$,
$J_{1}=1$, $J_{2}=0.5$, $J_{3}=0$, $\lambda=1.2$ and different values of $D$.
Table 1: Expectation value of the boson number for different temperatures and
spin. The theory parameters are $B=0$, $J_{1}=1$, $J_{2}=0.5$, $J_{3}=0$,
$\lambda=1.2$ and $D=0.2$. Temperature is in units of $J_{1}$. Figure 6: Spin
Nernst coefficient versus temperature for $S=1/2$, $B=0$, $J_{1}=1$,
$J_{2}=0.5$, $J_{3}=0$, $\lambda=1.2$ and different values of $D$.

In Figure 6 we present the spin Nernst coefficient as a function of
temperature. The coefficient $c_{1}(x)$ decreases with $x$ and this leads to a
flattening of $\alpha_{xy}$ for high temperature. We see a monotonic response
to the temperature, similar to other magnon systems like the AFM checkerboard
and FM Kagome lattices [28, 25], but in contrast to the FM checkerboard and
both FM and AFM honeycomb lattices [30, 22, 44].

The thermal Hall conductivity versus temperature is shown in Figure 7, and is
monotonically increasing with no sign change. This behaviour is similar to the
FM honeycomb, AFM checkerboard and AFM Kagome lattices [21, 29, 24]. The FM
Lieb and Kagome lattices also show the same behaviour for some choices of
interaction parameters; for other combinations, $\kappa_{xy}(T)$ has a local
minimum/maximum and can even change sign with the increase of the temperature
[31, 23]. This heterogeneous character of the $\kappa_{xy}(T)$ for different
parameters may be related to different topological phases of the insulating
system, indexed by the Chern numbers of the bands. As mentioned above, our
system is an insulator with a single trivial phase (the Chern number is zero
for any combination of parameters), so we wouldn’t expect any change in the
character of the curve for different parameters. The same can be said about
the other transport coefficients. In fact, the only change we detected in the
transport coefficients was from a quantitative nature (see Figure 8 for an
example).

Figure 7: Thermal Hall conductivity versus temperature for $S=1/2$, $B=0$,
$J_{1}=1$, $J_{2}=0.5$, $J_{3}=0$, $\lambda=1.2$ and different values of $D$.
Figure 8: Spin Hall conductivity versus temperature for $S=1/2$, $B=0$,
$J_{1}=J_{2}=J_{3}=1$, $D=0.2$ and $\lambda=1.2$. This plot exemplifies the
fact that the transport coefficients curves change only quantitatively for
different combination of parameters.

Figure 9 shows $\sigma_{xy}$ as a function of the applied magnetic field $B$
perpendicular to the lattice plane, for $D=0.2$ and $T=0.1$. From Eq. 14 we
see that a magnetic field does not affect the Berry curvature. However,
increasing $B$, $\omega_{+}$ and $\omega_{-}$ increase, and so, for a given
$T$, a smaller number of magnons are excited and $\sigma_{xy}$ decreases. A
strong magnetic field diminishes the thermal population difference between the
bands leading to a suppression of $\sigma_{xy}$, and so affecting
substantially the conductivity. The same behaviour is noted in the other
transport coefficients (Figures 10 and 11). This behaviour was predicted for
generic ferromagnetic 2D films in the dipolar regime [13], and was also
observed in theorical calculations on the checkerboard and Kagome lattices
[30, 14, 24].

Figure 9: Spin Hall conductivity versus magnetic field for $S=1/2$, $T=0.1$,
$J_{1}=1$, $J_{2}=0.5$, $J_{3}=0$, $\lambda=1.2$ and $D=0.2$. Figure 10: Spin
Nernst coefficient versus magnetic field for $S=1/2$, $T=0.1$, $J_{1}=1$,
$J_{2}=0.5$, $J_{3}=0$, $\lambda=1.2$ and $D=0.2$. Figure 11: Thermal Hall
conductivity versus magnetic field for $S=1/2$, $T=0.1$, $J_{1}=1$,
$J_{2}=0.5$, $J_{3}=0$, $\lambda=1.2$ and $D=0.2$.

##  5 Conclusions

We studied a two-band ferromagnetic magnon model with the geometry of a
modified Lieb lattice. An easy-axis anisotropy induces a gap in the
$X^{\prime}=(0,\pm\pi)$ point. When a uniform Dzyaloshinsky-Moriya interaction
is present between next-next-near neighbours, we find a non-vanishing Berry
curvature due to time-reversal symmetry breaking.

The Chern numbers of the bands are null, so the insulating system has only one
(trivial) phase. Nevertheless, the finite Berry curvature induces three magnon
Hall-like effects whose transport coefficients were studied: the spin Hall
effect, the thermal spin Hall effect and spin Nernst effect. The response of
the coefficients to the temperature is monotonic without sign change, and
resembles other topological magnon systems found on the literature. The
presence of an external off-plane magnetic field through a Zeeman interaction
minimizes the transport effects, as we expected by thermodynamical
considerations.

As far as we know, up to now there is no material described by the lattice
studied here. Nevertheless, there is a variety of compounds described by the
Lieb lattice. Thus, we believe that through a modification of the Lieb lattice
one could synthesize a compound where our model could be used. Another
possibility is in the field of optical lattices, where advances in
synthesizing techniques make it possible to mimic DM interactions using laser
beams.

## Acknowledgments

This work was supported by CAPES (Coordenação de Aperfeiçoamento de Pessoal de
Nível Superior) and CNPq (Conselho Nacional de Desenvolvimento Científico e
Tecnológico).

## References

  * [1] K. v. Klitzing, G. Dorda, and M. Pepper, “New method for high-accuracy determination of the fine-structure constant based on quantized hall resistance,” Phys. Rev. Lett., vol. 45, pp. 494–497, 1980.
  * [2] M. Fruchart and D. Carpentier, “An introduction to topological insulators,” C.R. Phys., vol. 14, no. 9-10, pp. 779–815, 2013.
  * [3] J. E. Hirsch, “Spin hall effect,” Phys. Rev. Lett., vol. 83, pp. 1834–1837, 1999.
  * [4] J. Sinova, D. Culcer, Q. Niu, N. A. Sinitsyn, T. Jungwirth, and A. H. MacDonald, “Universal intrinsic spin hall effect,” Phys. Rev. Lett., vol. 92, p. 126603, 2004.
  * [5] S. Murakami, N. Nagaosa, and S.-C. Zhang, “Dissipationless quantum spin current at room temperature,” Science, vol. 301, no. 5638, pp. 1348–1351, 2003.
  * [6] N. Nagaosa, “Anomalous hall effect –a new perspective–,” J. Phys. Soc. Jpn., vol. 75, no. 4, p. 042001, 2006.
  * [7] D. Xiao, Y. Yao, Z. Fang, and Q. Niu, “Berry-phase effect in anomalous thermoelectric transport,” Phys. Rev. Lett., vol. 97, p. 026603, 2006.
  * [8] D. Xiao, M.-C. Chang, and Q. Niu, “Berry phase effects on electronic properties,” Rev. Mod. Phys., vol. 82, no. 3, pp. 1959–2007, 2010.
  * [9] J. H. Han and H. Lee, “Spin chirality and hall-like transport phenomena of spin excitations,” J. Phys. Soc. Jpn., vol. 86, no. 1, p. 011007, 2017\.
  * [10] R. Matsumoto and S. Murakami, “Theoretical prediction of a rotating magnon wave packet in ferromagnets,” Phys. Rev. Lett., vol. 106, no. 19, p. 197202, 2011.
  * [11] R. Matsumoto and S. Murakami, “Rotational motion of magnons and the thermal hall effect,” Phys. Rev. B, vol. 84, no. 18, p. 184406, 2011.
  * [12] H. Katsura, N. Nagaosa, and P. A. Lee, “Theory of the thermal hall effect in quantum magnets,” Phys. Rev. Lett., vol. 104, p. 066403, 2010.
  * [13] R. Matsumoto, R. Shindou, and S. Murakami, “Thermal hall effect of magnons in magnets with dipolar interaction,” Phys. Rev. B, vol. 89, p. 054420, 2014\.
  * [14] H. Lee, J. H. Han, and P. A. Lee, “Thermal hall effect of spins in a paramagnet,” Phys. Rev. B, vol. 91, p. 125413, 2015.
  * [15] S. Park, N. Nagaosa, and B.-J. Yang, “Thermal hall effect, spin nernst effect, and spin density induced by a thermal gradient in collinear ferrimagnets from magnon–phonon interaction,” Nano Lett., vol. 20, no. 4, pp. 2741–2746, 2020.
  * [16] F. Meier and D. Loss, “Magnetization transport and quantized spin conductance,” Physical Review Letters, vol. 90, no. 16, 2003.
  * [17] S. Fujimoto, “Hall effect of spin waves in frustrated magnets,” Phys. Rev. Lett., vol. 103, p. 047203, 2009.
  * [18] A. V. Chumak, V. I. Vasyuchka, A. A. Serga, and B. Hillebrands, “Magnon spintronics,” Nat. Phys., vol. 11, no. 6, pp. 453–461, 2015.
  * [19] Y. Onose, T. Ideue, H. Katsura, Y. Shiomi, N. Nagaosa, and Y. Tokura, “Observation of the magnon hall effect,” Science, vol. 329, no. 5989, pp. 297–299, 2010.
  * [20] S. A. Owerre, “A first theoretical realization of honeycomb topological magnon insulator,” J. Phys.: Condens. Matter, vol. 28, no. 38, p. 386001, 2016\.
  * [21] S. A. Owerre, “Topological honeycomb magnon hall effect: A calculation of thermal hall conductivity of magnetic spin excitations,” J. Appl. Phys., vol. 120, no. 4, p. 043903, 2016.
  * [22] S. K. Kim, H. Ochoa, R. Zarzuela, and Y. Tserkovnyak, “Realization of the haldane-kane-mele model in a system of localized spins,” Phys. Rev. Lett., vol. 117, no. 22, p. 227201, 2016.
  * [23] A. Mook, J. Henk, and I. Mertig, “Magnon hall effect and topology in kagome lattices: A theoretical investigation,” Phys. Rev. B, vol. 89, p. 134409, 2014.
  * [24] P. Laurell and G. A. Fiete, “Magnon thermal hall effect in kagome antiferromagnets with dzyaloshinskii-moriya interactions,” Phys. Rev. B, vol. 98, p. 094419, 2018.
  * [25] A. A. Kovalev and V. Zyuzin, “Spin torque and nernst effects in dzyaloshinskii-moriya ferromagnets,” Phys. Rev. B, vol. 93, p. 161106, 2016\.
  * [26] M. Malki and G. S. Uhrig, “Topological magnon bands for magnonics,” Phys. Rev. B, vol. 99, p. 174412, 2019.
  * [27] M. Kawano and C. Hotta, “Thermal hall effect and topological edge states in a square-lattice antiferromagnet,” Phys. Rev. B, vol. 99, no. 5, p. 054422, 2019.
  * [28] A. S. T. Pires, “Magnon spin nernst effect on the antiferromagnetic checkerboard lattice,” Phys. Lett. A, vol. 383, no. 32, p. 125887, 2019\.
  * [29] A. S. T. Pires, “Topological magnons in the antiferromagnetic checkerboard lattice,” Physica E, vol. 118, p. 113899, 2020.
  * [30] A. S. T. Pires, “Topological magnons on the checkerboard lattice,” Physica B, vol. 602, p. 412490, 2021.
  * [31] X. Cao, K. Chen, and D. He, “Magnon hall effect on the lieb lattice,” J. Phys.: Condens. Matter, vol. 27, no. 16, p. 166003, 2015.
  * [32] S. Kambe and O. Ishii, “Correlation between tc, in-plane cuocu bond length, and buckling of the cuo2 plane in cuprate superconductors,” Physica C, vol. 341-348, pp. 555–556, 2000.
  * [33] M. Niţă, B. Ostahie, and A. Aldea, “Spectral and transport properties of the two-dimensional lieb lattice,” Phys. Rev. B, vol. 87, p. 125428, 2013.
  * [34] C. Weeks and M. Franz, “Topological insulators on the lieb and perovskite lattices,” Phys. Rev. B, vol. 82, p. 085310, 2010.
  * [35] W. Beugeling, J. C. Everts, and C. Morais Smith, “Topological phase transitions driven by next-nearest-neighbor hopping in two-dimensional lattices,” Phys. Rev. B, vol. 86, p. 195129, 2012.
  * [36] Y. He, J. Moore, and C. M. Varma, “Berry phase and anomalous hall effect in a three-orbital tight-binding hamiltonian,” Phys. Rev. B, vol. 85, p. 155106, 2012.
  * [37] I. Dzyaloshinsky, “A thermodynamic theory of “weak” ferromagnetism of antiferromagnetics,” J. Phys. Chem. Solids, vol. 4, no. 4, pp. 241–255, 1958.
  * [38] T. Moriya, “Anisotropic superexchange interaction and weak ferromagnetism,” Phys. Rev., vol. 120, pp. 91–98, Oct. 1960.
  * [39] H. Katsura, N. Nagaosa, and A. V. Balatsky, “Spin current and magnetoelectric effect in noncollinear magnets,” Phys. Rev. Lett., vol. 95, no. 5, p. 057205, 2005.
  * [40] C. Herring and C. Kittel, “On the theory of spin waves in ferromagnetic media,” Phys. Rev., vol. 81, no. 5, pp. 869–880, 1951.
  * [41] P. Bruno, “Spin-wave theory of two-dimensional ferromagnets in the presence of dipolar interactions and magnetocrystalline anisotropy,” Phys. Rev. B, vol. 43, no. 7, pp. 6015–6021, 1991.
  * [42] F. J. Dyson, “General theory of spin-wave interactions,” Phys. Rev., vol. 102, no. 5, pp. 1217–1230, 1956.
  * [43] T. Oguchi, “Theory of spin-wave interactions in ferro- and antiferromagnetism,” Phys. Rev., vol. 117, no. 1, pp. 117–123, 1960.
  * [44] R. Cheng, S. Okamoto, and D. Xiao, “Spin nernst effect of magnons in collinear antiferromagnets,” Phys. Rev. Lett., vol. 117, no. 21, p. 217202, 2016\.


11institutetext: Leiden University, The Netherlands 11email:
<EMAIL_ADDRESS>22institutetext: University of
Oxford, United Kingdom 22email<EMAIL_ADDRESS>33institutetext: RWTH Aachen University, Germany
33email<EMAIL_ADDRESS>
# Automated Design of Linear Bounding Functions
for Sigmoidal Nonlinearities in Neural Networks

Matthias König Work done while at University of Oxford.11 Xiyue Zhang (✉) 22
Holger H. Hoos 1133 Marta Kwiatkowska 22 Jan N. van Rijn 11

###### Abstract

The ubiquity of deep learning algorithms in various applications has amplified
the need for assuring their robustness against small input perturbations such
as those occurring in adversarial attacks. Existing complete verification
techniques offer provable guarantees for all robustness queries but struggle
to scale beyond small neural networks. To overcome this computational
intractability, incomplete verification methods often rely on convex
relaxation to over-approximate the nonlinearities in neural networks. Progress
in tighter approximations has been achieved for piecewise linear functions.
However, robustness verification of neural networks for general activation
functions (_e.g._ , Sigmoid, Tanh) remains under-explored and poses new
challenges. Typically, these networks are verified using convex relaxation
techniques, which involve computing linear upper and lower bounds of the
nonlinear activation functions. In this work, we propose a novel parameter
search method to improve the quality of these linear approximations.
Specifically, we show that using a simple search method, carefully adapted to
the given verification problem through state-of-the-art algorithm
configuration techniques, improves the average global lower bound by 25% on
average over the current state of the art on several commonly used local
robustness verification benchmarks.

###### Keywords:

Neural Network Verification Automated Algorithm Configuration Convex
Relaxation

##  1 Introduction

(a) $\mathcal{S}^{+}$

(b) $\mathcal{S}^{-}$

(c) $\mathcal{S}^{\pm}$

Figure 1: Linear bounding rules for different cases of the Sigmoid activation
function. The x-axis shows the pre-activation bounds, while the y-axis
indicates the output of the activation function.

Over the last decade, deep learning algorithms have gained increasing
significance as essential tools across diverse application domains and usage
scenarios. Their applications range from manoeuvre advisory systems in
unmanned aircraft to face recognition in mobile phones (see, _e.g._ , [18]).
Simultaneously, it is now widely acknowledged that neural networks are
susceptible to adversarial attacks, where a given input is manipulated to
cause misclassification [31]. Remarkably, in image recognition tasks, the
perturbation required can be so subtle that it remains virtually imperceptible
to the human eye.

Numerous methods have been proposed to evaluate the robustness of neural
networks against adversarial attacks. Early methods involve empirical attacks
[12, 22, 7]; however, these approaches do not provide a comprehensive
assessment of neural network robustness, as a defence mechanism against one
type of attack might still be vulnerable to another, potentially novel, class
of attacks. Consequently, there have been efforts to develop approaches that
formally verify the robustness of neural networks against adversarial attacks
[19, 11, 37, 6, 33, 3, 13]. These formal verification methods enable a
principled assessment of neural network robustness, providing provable
guarantees on desirable properties, typically in the form of a pair of input-
output specifications. However, complete verification of neural networks,
where the verifier is theoretically guaranteed to provide a definite answer to
the property under verification, is a challenging NP-complete problem [19].
Solving this problem usually requires computationally expensive methods,
_e.g._ , SMT solvers [19] or mixed integer linear programming systems [34],
limiting scalability and efficiency.

The aforementioned computational complexity is mainly due to the nonlinear
activation functions in a neural network, which results in the neural network
verification problem becoming a non-convex optimisation problem. In light of
this, incomplete methods have been proposed that exploit _convex relaxation_
techniques, which approximate the nonlinearities using linear symbolic bounds,
to provide sound and efficient verification [1, 30, 39, 11, 28]; completeness
can be achieved through branch and bound (see, _e.g._ [5]). Most formal
verification methods are limited to ReLU-based networks (see [23]), which
satisfy the piecewise linear property. However, convex relaxation techniques
are applicable to _more general_ commonly-used activation functions, such as
Sigmoid or Tanh, by approximating them in terms of piecewise linear functions.
This enables an extension of formal verification methods based on convex
relaxation to general activation functions, though these remain under-
explored.

An important application of neural network verification tools is _robustness
certification_ , which computes guarantees that the prediction of the network
is stable (invariant) around a given input point. $\mathrm{CROWN}$ [39] is the
first generic framework leveraging adaptive linear bounds to efficiently
certify robustness for general activation functions. Another series of
techniques, _e.g._ , $\mathrm{DeepZ}$ [28], $\mathrm{DeepPoly}$ [30] employs
abstract interpretation and abstract transformers for commonly-used activation
functions, such as Sigmoid and Tanh. $\mathrm{CROWN}$ computes linear upper
and lower bounds of nonlinear functions, while $\mathrm{DeepZ}$ and
$\mathrm{DeepPoly}$ incorporate convex relaxation into the abstract
transformer.

Convex relaxation methods typically over-approximate the activation functions
in all nonlinear neurons in the neural network, which inevitably introduces
imprecision. Consequently, the corresponding verification algorithms may fail
to prove the robustness of a neural network when the original network
satisfies the specification and thus weaken certification guarantees. A
globally, _i.e._ , network-wise, tighter over-approximation is crucial for
strong certification guarantees. However, for Sigmoidal activations, it
remains unclear how to design or configure the linear approximation of
nonlinear neurons to achieve globally tight output bounds, which are directly
used to determine the robustness of neural networks with respect to a specific
property under investigation.

Motivated by the need to reduce imprecision of the bounding functions (as an
over-approximation of the activation function), in this work we introduce an
automated and systematic method to compute tighter bounding functions to
improve certification guarantees. These bounding functions are defined by the
tangent point, where they touch the activation function. To this end, we
propose a novel parameter search method for identifying the tangent points to
find better-suited linear bounding functions by considering different cases of
Sigmoidal activation functions. To tackle the infinite search space (of
tangent points) for the bounding functions, our approach leverages state-of-
the-art algorithm configuration techniques.

Concretely, we use automated algorithm configuration techniques to find
optimal hyper-parameters of the search method used for obtaining the tangent
points of the linear bounding functions, which has previously been done using
binary search [39]. These hyper-parameters control the initial tangent point
per neuron as well as the rate at which these initial points are updated (we
will refer to the latter as a multiplier for the remainder of this work) until
a feasible bounding function has been obtained. Notice that these hyper-
parameters are set for the entire network, _i.e._ , all neurons share the same
starting point and multiplier; however, they eventually result in different
tangent points, as the search process only ends once a feasible bounding
function has been found.

Moreover, we show that, by using our proposed method, we improve the average
lower bound on the network output by 25% on average across several
verification benchmarks, and can certify robustness for instances that were
previously unsolved.

##  2 Related Work

In the following, we give some background on using convex relaxation for
neural network robustness verification and on automated algorithm
configuration.

###  2.1 Convex Relaxation for Neural Network Verification

We use $f:\mathbb{R}^{n}\to\mathbb{R}^{m}$ to denote a neural network trained
to make predictions on an $m$-class classification problem. Let
$\textbf{W}^{(i)}$ denote the weight matrix and $\textbf{b}^{(i)}$ denote the
bias for the $i$-th layer. We use $\hat{z}^{(i)}$ to denote the pre-activation
neuron values, and $z^{(i)}$ to denote the post-activation ones, such that we
have $\hat{z}^{(i)}=\textbf{W}^{(i)}z^{(i-1)}+\textbf{b}^{(i)}$ (_i.e._ , the
post-activation values of layer $i-1$ are linearly weighted to become the pre-
activation values of layer $i$). We use $\sigma$ to denote the activation
functions in intermediate layers, such that $z=\sigma(\hat{z})$. In this work,
we focus on neural networks with sigmoidal activation functions, including
Sigmoid $\sigma(x)=1/(1+e^{-x})$ and Tanh
$\sigma(x)=(e^{x}-e^{-x})/(e^{x}+e^{-x})$. We use $f(x)$ to denote the neural
network output for all classes and $f_{j}(x)$ to denote the output associated
with class $j$ (with $1\leq j\leq m$). The final decision is computed by
taking the label with the greatest output value, _i.e._ ,
$\operatorname*{arg\,max}_{j}f_{j}(x)$.

Neural network verification. The robustness to adversarial perturbations is
one of the most important properties of neural networks, requiring that the
predictions of a neural network should be preserved for a local input region
$\mathcal{C}$, typically defined as a $l_{p}$ norm ball with a radius of
$\epsilon$ around the original input $x_{0}$, in the following way:

$\forall
x\in\mathcal{C},\,\operatorname*{arg\,max}_{j}f_{j}(x)=\operatorname*{arg\,max}_{j}f_{j}(x_{0}).$
(1)

For an input instance $x_{0}$ with ground-truth label $y_{0}$, verifying the
robustness property can then be transformed into proving that $\forall
x\in\mathcal{C},f_{y_{0}}(x)-f_{j}(x)\geq 0$ for all $j\neq y_{0}$ where
$j\in[1,m]$.

Given a verification problem, for example, to check whether the output
constraint $f_{y_{0}}(x)-f_{j}(x)\geq 0$ is satisfied, we can append an
additional layer at the end of the neural network, such that the output
property can be merged into the new neural network function
$g:\mathbb{R}^{n}\to\mathbb{R}$ where $g(x)=f_{y_{0}}(x)-f_{j}(x)$. While this
should formally be done for each possible network output $j\neq y_{0}$, to
simplify the notation we use a single $g$.

In this way, the verification problem can be formulated canonically as
follows, i.e., to prove or falsify:

$\forall x\in\mathcal{C},g(x)\geq 0$ (2)

One way to verify the property is to solve the optimisation problem
$\min_{x\in\mathcal{C}}g(x)$. However, due to the nonlinearity of the
activation function $\sigma$, the neural network verification problem is NP-
complete [19]. To address such intractability, state-of-the-art verification
algorithms leverage convex relaxation to transform the verification problem
into a convex optimisation problem.

######  Definition 1 (Convex relaxation of activation functions)

For a nonlinear activation function $\sigma(\hat{z})$ with pre-activation
bounds $\hat{z}\in[l,u]$, convex relaxation relaxes the non-convex equality
constraint $z=\sigma(\hat{z})$ to two convex inequality constraints by
computing the linear lower and upper bounding functions
$h_{L}(\hat{z})=\alpha_{L}\hat{z}+\beta_{L}$ and
$h_{U}(\hat{z})=\alpha_{U}\hat{z}+\beta_{U}$, such that
$h_{L}(\hat{z})\leq\sigma(\hat{z})\leq h_{U}(\hat{z})$.

In Figure 1, the pre-activation bounds are indicated by the grey-shaded area,
$h_{L}$ is illustrated by the red dotted line, whereas $h_{U}$ is indicated by
the green dotted line. Notice that the parameters
$\alpha_{L},\beta_{L},\alpha_{U},\beta_{U}$ of the linear functions depend on
the pre-activation bounds $l$ and $u$. With such convex inequality
constraints, we can propagate the relaxations through layers using an
efficient back-substitution procedure to compute the linear lower bounds
(denoted by $\underline{g}$) and upper bounds (denoted by $\overline{g}$) for
the neural network $g$.

######  Definition 2 (Convex relaxation of neural networks)

A convex relaxation of the neural network $g:\mathbb{R}^{n}\to\mathbb{R}$ over
an input region $\mathcal{C}$, are two linear functions $\underline{g}$ and
$\overline{g}$ such that $\underline{g}(x)\leq g(x)\leq\overline{g}(x)$ for
all $x\in\mathcal{C}$.

Using convex relaxation, the verification problem is then reduced to the
following convex optimisation problem:

$g^{\ast}=\min_{x\in\mathcal{C}}\underline{g}(x)$ (3)

The robustness property is proved if the optimal solution $g^{\ast}\geq 0$, as
we have $\forall x\in\mathcal{C},\,g(x)\geq\underline{g}(x)$. Meanwhile,
convex relaxation of the nonlinear constraints inevitably introduces
approximation. As a consequence, the computed minimum might fail to satisfy
$\min_{x\in\mathcal{C}}\underline{g}(x)\geq 0$ even in cases in which the
network is robust, _i.e._ , $\min_{x\in\mathcal{C}}g(x)\geq 0$. Instead, the
true verification remains unknown. Therefore, improving the quality of the
approximation is crucial to reducing false failures in robustness
verification, and thus strengthening certification guarantees, which is the
aim of our contribution.

###  2.2 Automated Algorithm Configuration

In general, the algorithm configuration problem can be described as follows:
Given an algorithm $A$ (also referred to as the target algorithm) with
parameter configuration space $\Theta$, a set of problem instances $\Pi$, and
a cost metric $c:\Theta\times\Pi\rightarrow\mathbb{R}$, find a configuration
$\theta^{\ast}\in\Theta$ that minimises cost $c$ across the instances in
$\Pi$:

$\theta^{\ast}\in\underset{\theta\in\Theta}{\mathrm{arg\,min}}\underset{\pi\in\Pi}{\sum}c(\theta,\pi)$
(4)

The general workflow of the algorithm configuration procedure starts with
selecting a configuration $\theta\in\Theta$ and an instance $\pi\in\Pi$. Next,
the configurator initialises a run of algorithm $A$ with configuration
$\theta$ on instance $\pi$, and measures the resulting cost $c(\theta,\pi)$.
The configurator uses this information about the target algorithm’s
performance to find a configuration that performs well on the training
instances. This is enabled by a surrogate model, which provides a posterior
probability distribution that characterises the potential cost for
$c(\theta,\pi)$ at a configuration $\theta$. At every time point $t$, we have
a new observation of the cost $c_{t}$ at a new configuration point
$\theta_{t}$. The posterior distribution is updated based on the augmented
observation set $S_{t}=S_{t-1}\cup\\{(\theta_{t},c_{t})\\}$, and the next
configuration point is selected by maximising the acquisition function in the
following form:

$\theta_{t+1}=\underset{\theta\in\Theta}{\mathrm{arg\,max}}\,a_{t}(\theta,S_{t})$
(5)

where $a_{t}$ denotes the acquisition function. A typical choice for the
acquisition functions is the expected improvement [16].

Once the configuration budget (e.g., time budget or number of trials) is
exhausted, the procedure returns the current incumbent $\theta^{\ast}$, which
represents the best configuration found so far. Finally, when running the
target algorithm with configuration $\theta^{\ast}$, it should result in lower
cost (such as average running time) or improved solution quality across the
benchmark set.

Automated algorithm configuration has already been shown to work effectively
in the context of formal neural network verification [21], but also in other,
related domains, such as SAT solving [14, 17], scheduling [8], mixed-integer
programming [15, 24], evolutionary algorithms [2], answer set solving [10], AI
planning [35], and machine learning [32, 9].

##  3 Method

We consider the task of local robustness verification with Sigmoidal
activation functions. To this end, we focus on convex relaxation based
perturbation analysis as employed in the $\mathrm{CROWN}$ framework [39], and
use automated algorithm configuration techniques to improve the linear bounds
of the nonlinear activation functions. Notice that, by default,
$\mathrm{CROWN}$ employs binary search to obtain the points at which the
linear bounding functions are tangent to the activation function [39].

In this study, on the other hand, we used $\mathrm{SMAC}$ [16] to guide the
search for suitable tangent points. $\mathrm{SMAC}$ is a widely known, freely
available, state-of-the-art configurator based on sequential model-based
optimisation (also known as Bayesian optimisation). The main idea of
$\mathrm{SMAC}$ is to construct and iteratively update a statistical model of
target algorithm performance to guide the search for promising configurations;
$\mathrm{SMAC}$ uses a Random Forest regressor [4].

###  3.1 Configuration Objective

The objective of the configuration procedure is to maximise the minimum of the
global lower bound $g^{\ast}$ as defined in Equation (3) for each instance
under verification. Notice that the optimisation procedure that solves
Equation (3) is performed separately for each instance, _i.e._ , $|\Pi|=1$.
Recall that an instance is verified to be robust when $g^{\ast}>0$, as
outlined in Section 2.1.

By design, $\mathrm{SMAC}$ solves a minimisation problem; see Equation (4).
Since we are interested in maximising the lower bound of the network output,
we apply appropriate sign changes and define the cost metric $c$ as the
negative of the global lower bound.

###  3.2 Configuration Space

The $\mathrm{CROWN}$ framework [39] proposed a general certification solution
for neural networks with nonlinear activation functions, which relies on
convex relaxation to compute the output bounds. Since Sigmoidal activation
functions (Sigmoid/Tanh/Arctan) share the same features, that is, convex on
the negative side ($x<0$) and concave on the other side ($x>0$), the authors
of [39] leveraged this feature and proposed a general method to compute the
parameters (_i.e._ , the tangent points) of linear upper and lower functions
$h_{U}$, $h_{L}$.

Based on the curvature of Sigmoidal activation functions, each nonlinear
neuron of the $i$-th layer (denoted as $[n_{i}]$) is categorised into one of
the three cases: $\mathcal{S}^{+}$, $\mathcal{S}^{-}$, and $\mathcal{S}^{\pm}$
where $\mathcal{S}^{+}=\\{j\in[n_{i}]\mid 0\leq l_{j}^{(i)}\leq
u_{j}^{(i)}\\}$, $\mathcal{S}^{-}=\\{j\in[n_{i}]\mid l_{j}^{(i)}\leq
u_{j}^{(i)}\leq 0\\}$, and $\mathcal{S}^{\pm}=\\{j\in[n_{i}]\mid
l_{j}^{(i)}\leq 0\leq u_{j}^{(i)}\\}$. Intuitively, $\mathcal{S}^{+}$
represents the case in which the pre-activation bounds are both positive,
$\mathcal{S}^{-}$ represents the case in which they are both negative and
$\mathcal{S}^{\pm}$ represents the case in which $l$ is negative and $u$ is
positive. Different bounding rules are then proposed for these three cases;
these are illustrated in Figure 1.

####  3.2.1 Bounding rules for $S^{+}$ domain

For any node $j\in\mathcal{S}^{+}$, the activation function
$\sigma(\hat{z}_{j})$ is concave; hence, a tangent line of
$\sigma(\hat{z}_{j})$ (represented in green in Figure 1) at any tangent point
$\hat{z}^{\ast}_{j}\in[l_{j}^{(i)},u_{j}^{(i)}]$ is a valid upper bounding
function $h_{U,j}^{(i)}$, and the linear function passing the two endpoints,
$(l_{j}^{(i)},\sigma(l_{j}^{(i)}))$ and $(u_{j}^{(i)},\sigma(u_{j}^{(i)}))$,
serves as a valid lower bounding function $h_{L,j}^{(i)}$.

Algorithm 1 Our proposed tangent point search method for sigmoidal functions.
It details only the search method in the $S^{+}$ domain (for space reasons).
The hyper-parameter $s$ determines the initial tangent point and the hyper-
parameter $\psi$ determines the rate by which $s$ is changed. The values of
$s$ and the update rate $\psi$ are not exposed as hyper-parameters in vanilla
$\mathrm{CROWN}$. In our case, these are optimised by $\mathrm{SMAC}$, where
$\mathrm{SMAC}$ uses a fixed budget of $n_{\mathtt{max}}$ evaluation calls.
Notice that in vanilla $\mathrm{CROWN}$, binary search is performed once
initial bounds have been obtained to move tangent points closer to 0; this
step is not performed in our method.

1:procedure OptimiseSearchParameters($x$, $s$, $\psi$)

2: Initialize cost $c_{0}=\inf$ $\triangleright$ first iteration of the
optimisation loop can overwrite it

3: Initialise observation set $S_{0}=\\{\emptyset\\}$

4: for $n=1$, $\cdots$, $n_{\mathtt{max}}$ do

5: Select starting point and multiplier rate based on observation set

6:
$(s_{n},\psi_{n})\leftarrow\underset{(s,\psi)\in\Theta}{\mathrm{arg\,max}}\,a_{n}((s,\psi),S_{n-1})$

7: for all node $j$ in network $f$ do

8: if $j\in S^{+}$ then

9: Initialise upper bound tangent point for current node
$\hat{z}^{\ast}_{j}\leftarrow s_{n}$

10: while Tangent point $\hat{z}^{\ast}_{j}$ leads to invalid upper bound
$h_{U,j}^{(i)}$ do

11: Update $\hat{z}^{\ast}_{j}$:
$\hat{z}^{\ast}_{j}\leftarrow\hat{z}^{\ast}_{j}\cdot\psi_{n}$

12: end while

13: else if $j\in S^{-}$ then

14: Determine tangent point for lower bounding function, in a similar fashion
as above

15: else$\triangleright$ In this case, $j\in S^{\pm}$;

16: Determine tangent point for both lower bounding and upper bounding
function, in a similar fashion as above

17: end if

18: end for

19: Evaluate cost $c_{n}=-1\cdot g^{\ast}(x)$ $\triangleright$ Linear bounding
layer-by-layer via convex relaxation as per Def. 2 and Eq. 3

20: Augment observation set
$S_{n}=S_{n-1}\cup\bigl{\\{}\bigl{(}(s_{n},\psi_{n}),c_{n}\bigr{)}\bigr{\\}}$

21: Train surrogate model using $S_{n}$

22: if $c_{n}<c^{\ast}$ then

23: $s^{\ast}\leftarrow s_{n}$, $\psi^{\ast}\leftarrow\psi_{n}$,
$c^{\ast}\leftarrow c_{n}$

24: end if

25: end for

26: return $c^{\ast}$, ($s^{\ast},\psi^{\ast}$)

27:end procedure

To select good values of $\hat{z}^{\ast}_{j}$, we propose a search method
whose behaviour depends on two hyper-parameters, which are optimised by SMAC:

  * •

starting point $s$ at which each tangent point $\hat{z}^{\ast}_{j}$ is
initialised;

  * •

multiplier $\psi$ to change the value of $\hat{z}^{\ast}_{j}$ if the linear
bound is found to be invalid.

This search method is illustrated in Algorithm 1. The algorithm initiates and
maintains the best configuration $\theta^{\ast}=(s^{\ast},\psi^{\ast})$ found
so far, which leads to a tighter bounding function (defined by the tangent
points $\hat{z}^{\ast}_{j}$ for each node $j$) and, thus, a tighter global
lower bound $g^{\ast}$. To find optimal settings for the hyper-parameters $s$
and $\psi$, SMAC samples values of $s$ and $\psi$ from a pre-defined
configuration space $\Theta$. Configurations $\theta$, as introduced in
Equation (4), are sampled from $s\in[0.01,2]$ and $\psi\in[1.01,3]$. Note that
these ranges are based on empirical observations; values outside this range
typically result in extremely loose bounds. In principle, other ranges could
be used to sample values of $s$ and $\psi$, respectively, as long as $s>0$ and
$\psi>1$. Also note that the search method will be the same for the entire
network; however, it will result in different tangent points per node.

Based on a given hyperparameter configuration, we can identify the tangent
points $\hat{z}^{\ast}_{j}$ for each node $j$, and compute the global lower
bound $g^{\ast}$ as well as the cost value $c=-1\cdot g^{\ast}(x)$. At each
iteration, we can acquire the parameter configuration based on the updated
surrogate model (trained with the new observation set $S_{n}$). With the
selected candidate configuration, we perform a sanity check to ensure the
validity of the linear upper bounding function incurred by the tangent points
$\hat{z}^{\ast}_{j}$ at each node. An upper bound is considered invalid if, at
any given point $\hat{z}^{\ast}_{j}$, the value of
$h_{U,j}^{(i)}(\hat{z}^{\ast}_{j})$ is smaller than the value of the nonlinear
sigmoidal $\sigma(\hat{z}^{\ast}_{j})$. The global lower bound $g_{n}$ (also
the cost value $c_{n}$) is then updated with the new upper bounding function.
The loop terminates and returns the best-achieved cost value $c^{\ast}$ as
well as the corresponding configuration $(s^{\ast},\psi^{\ast})$ when the
maximum iteration limit $n_{\mathtt{max}}$ is reached. Notice that SMAC has
several additional options that can slightly deviate from the description
[16].

While Algorithm 1 elaborates on tangent point search for the upper bounding
function, it can be similarly configured to search for lower bounding
functions with two key modifications. Firstly, the parameters are initialised
to different value ranges. Secondly, the validity evaluation is to check the
bounding function will always return a lower value than the sigmoidal
function. We introduce the details for $S^{-}$ and $S^{\pm}$ domains in the
following.

####  3.2.2 Bounding rules for $S^{-}$ domain

Symmetrically, for any node $j\in\mathcal{S}^{-}$, $h_{U,j}^{(i)}$ is defined
as the linear function passing the two endpoints and $h_{L,j}^{(i)}$ can be a
tangent line of $\sigma(\hat{z}_{j})$ at any point
$\hat{z}^{\ast}_{j}\in[l_{j}^{(i)},u_{j}^{(i)}]$. In this case, we employ a
similar search method, except that a configuration $\theta^{*}$ is now sampled
from $s\in[-0.01,-2]$ and $\psi\in[-1.01,-3]$. Again, other ranges could be
used to sample values of $s$ and $\psi$, respectively, as long as $s<0$ and
$\psi<-1$. Furthermore, a bound is considered invalid if, at any given point
$\hat{z}^{\ast}_{j}$, the value of $h_{L,j}^{(i)}(\hat{z}^{\ast}_{j})$ is
larger than the value of $\sigma(\hat{z}^{\ast}_{j})$.

####  3.2.3 Bounding rules for $S^{\pm}$ domain

Lastly, for any node $j\in\mathcal{S}^{\pm}$, $h_{U,j}^{(i)}$ is a tangent
line passing $(l_{j}^{(i)},\sigma(l_{j}^{(i)}))$ and a tangent point
$(\hat{z}^{\ast}_{1},\sigma(\hat{z}^{\ast}_{1}))$, where
$\hat{z}^{\ast}_{1}\geq 0$, and $h_{L,j}^{(i)}$ is a tangent line that passes
$(u_{j}^{(i)},\sigma(u_{j}^{(i)}))$ and a tangent point
$(\hat{z}^{\ast}_{2},\sigma(\hat{z}^{\ast}_{2}))$, where
$\hat{z}^{\ast}_{2}\leq 0$. Again, we employ the same search method to obtain
$\hat{z}^{\ast}_{j}$ for each bound. Moreover, for any $h_{U,j}^{(i)}$, a
configuration $\theta^{*}$ is sampled from $s\in[0.01,2]$ and
$\psi\in[1.01,3]$, while for any $h_{L,j}^{(i)}$, a configuration $\theta^{*}$
is sampled from $s\in[-0.01,-2]$ and $\psi\in[-1.01,-3]$. Notice that we can
use similar rules for bounding Tanh functions, which share the Sigmoidal shape
but range from -1 to 1 instead of 0 to 1.

Table 1: Experimental results obtained from our automated bound configuration
method, compared to the vanilla $\mathrm{CROWN}$ algorithm, which relies on
binary search to obtain the tangent points (with $\alpha$-optimisation for
Sigmoid networks). Boldfaced values indicate superior performance.

|  |  |  |  Avg. Global Lower Bound $\mathbf{g^{\ast}}$ | # Certified Instances
---|---|---|---|---|---
Dataset | Network | Activation | Epsilon | Baseline | Configured | Improvement | Baseline | Configured
CIFAR10 | ConvMed | Sigmoid | 0.0313 | -15.611 | -5.506 | 184% | 0 | 1
CIFAR10 | ConvMed | Sigmoid | 0.0157 | -1.727 | -1.633 | 6% | 13 | 13
CIFAR10 | ConvMed | Sigmoid | 0.0078 | -0.709 | -0.702 | 1% | 29 | 31
CIFAR10 | ConvMed | Sigmoid | 0.0039 | -0.440 | -0.440 | - | 39 | 39
CIFAR10 | FNN_6x500 | Sigmoid | 0.0313 | -18.448 | -16.433 | 12% | 0 | 0
CIFAR10 | FNN_6x500 | Sigmoid | 0.0157 | -17.819 | -16.234 | 10% | 0 | 0
CIFAR10 | FNN_6x500 | Sigmoid | 0.0078 | -15.075 | -13.478 | 12% | 0 | 0
CIFAR10 | FNN_6x500 | Sigmoid | 0.0039 | -6.731 | -5.915 | 14% | 5 | 8
MNIST | ConvMed | Sigmoid | 0.3 | -5.153 | -4.175 | 23% | 3 | 4
MNIST | ConvMed | Sigmoid | 0.12 | -9.484 | -8.563 | 10% | 0 | 0
MNIST | ConvMed | Sigmoid | 0.06 | -0.110 | -0.070 | 57% | 54 | 54
MNIST | ConvMed | Sigmoid | 0.03 | 2.753 | 2.756 | - | 92 | 92
CIFAR10 | ConvMed | Tanh | 0.0313 | -65.849 | -55.565 | 19% | 0 | 0
CIFAR10 | ConvMed | Tanh | 0.0157 | -44.504 | -33.723 | 32% | 0 | 0
CIFAR10 | ConvMed | Tanh | 0.0078 | -10.835 | -9.348 | 16% | 2 | 2
CIFAR10 | ConvMed | Tanh | 0.0039 | -2.299 | -2.259 | 2% | 15 | 16

##  4 Setup for Empirical Evaluation

We will empirically investigate the effectiveness of the proposed search
procedure. We consider two types of neural network architecture: convolutional
neural networks (CNNs) and fully connected neural networks (FNNs). We evaluate
the effectiveness of our approach on the $\mathrm{ERAN}$ benchmark [25, 26,
30, 27, 28]. Following the naming conventions, we refer to these networks as
ConvMed and FNN_6x500. These networks are adversarially trained on the
CIFAR-10 dataset with $\epsilon=0.0313$. In addition, to investigate whether
our findings hold for additional datasets, we consider a CNN with similar
architecture that is adversarially trained on the MNIST dataset with
$\epsilon=0.3$. Furthermore, to demonstrate the generality of our method to
other Sigmoidal activation functions, we perform the evaluation on a neural
network trained on CIFAR-10 with similar architecture, adversarially trained
with $\epsilon=0.0313$, but using Tanh activation functions. Further details
about the considered networks can be found in the ERAN repository.

We verify each network for local robustness with respect to the first 100
instances in the test set of the MNIST and CIFAR-10 datasets, respectively,
initially with perturbation radii equivalent to the values used during
training. We also evaluate our method on a broader range of perturbation
radii; specifically, we considered $\epsilon\in\\{0.0157,0.0078,0.0039\\}$ for
CIFAR-10 and $\epsilon\in\\{0.12,0.06,0.03\\}$ for MNIST. Notice that these
values of $\epsilon$ are in line with commonly chosen values from the
verification literature [20, 36, 29].

For verification, we employ $\mathrm{CROWN}$ with default settings as provided
by the authors. For Sigmoid-based networks, $\mathrm{CROWN}$ also performs
$\alpha$-optimisation [38]; however, this is not implemented for Tanh
activation functions, where we run $\mathrm{CROWN}$ without further
optimisation. As we are purely interested in the global lower bound on the
network output, we skip the PGD attack used for upper bound computation. For
the configuration procedure, we set the number of evaluation calls performed
by SMAC to 150, _i.e._ , the configuration process terminates after 150
trials.

(a) CIFAR10, CNN,
$\epsilon=0.0313$

(b) CIFAR10, CNN,
$\epsilon=0.0078$

(c) CIFAR10, FNN,
$\epsilon=0.0039$

(d) MNIST, CNN,
$\epsilon=0.3$

Figure 2: Experimental results obtained for Sigmoid-based networks. Each dot
represents a problem instance and the global lower bound, _i.e._ , the value
of $g^{\ast}$, for that instance achieved by the baseline approach (x-axis)
_vs_ our method (y-axis).

All experiments are performed on a cluster of machines equipped with NVIDIA
GeForce GTX 1080 Ti GPUs with 11 GB video memory.

##  5 Experimental Results and Discussion

(a) Layer 1

(b) Layer 2

(c) Layer 3

Figure 3: Probability density functions of the values of $\hat{z}^{\ast}_{j}$
obtained by our method as well as vanilla $\mathrm{CROWN}$ for lower and upper
bounding functions per activation layer. Remember that $\hat{z}^{\ast}_{j}$
determines the tangent point of the bounding function of a given node $j$.

Table 1 shows the average global lower bound over the given test instances
achieved by the vanilla $\mathrm{CROWN}$ algorithm, which represents our
baseline, and those achieved by $\mathrm{CROWN}$ in combination with our
configured bounding method. In addition, we report the absolute number of
instances for which robustness certification could be obtained by each
approach.

###  5.1 Sigmoid-based Networks

We first report the results obtained for the Sigmoid-based CNN trained on the
CIFAR-10 dataset. These can be found in Table 1. When certifying this network
with $\epsilon=0.03$, we obtained an improvement in the average global lower
bound of 184% (-15.611 _vs_ -5.506). This is also visualised in Figure 2(a).
As shown there, the global lower bound is improved by almost an order of
magnitude for some of the instances. Furthermore, using our configured bounds,
we could certify robustness for an instance that could not be solved by the
baseline approach.

For $\epsilon=0.0157$ and $\epsilon=0.0078$, we achieved improvements of 6%
and 1%, respectively. Interestingly, although the latter is a rather small
improvement in the average global lower bound, our method could verify two
additional instances, which the baseline method was unable to solve. These
results are visualised in Figure 2(b). Although instances generally lie very
close to the equality line, our method could improve on instances with lower
bounds very close to 0, for which even a very small increase can lead to
certified robustness.

Lastly, when $\epsilon=0.0039$, we achieved a similar performance as the
baseline method. This indicates that the effectiveness of configured bounds
decreases as the perturbation radius becomes very small.

Next, we investigate whether our approach extends to fully connected neural
networks; experimental results are shown in Table 1. In general, bounds
obtained for this network type are looser than those obtained for the CNN,
irrespective of the perturbation radius. Nonetheless, our method achieved
improvements in the average global lower bound between 10 and 14 per cent
across all considered perturbation radii.

For FNNs, we achieved the greatest improvement when the perturbation radius is
smallest, _i.e._ $\epsilon=0.0039$. In this scenario, we were able to certify
robustness for 3 additional instances, which were previously unsolved; see
also Figure 2(c) for a visualisation of these results.

Table 1 also shows the results from our experiments on the CNN trained on the
MNIST dataset. When verifying this network with $\epsilon=0.3$, we achieved an
improvement in the average global lower bound of 23%. Furthermore, using our
configured bounds, we could again certify robustness for an instance that
could not be solved previously. This is also visualised in Figure 2(d).

For $\epsilon=0.12$ and $\epsilon=0.06$, we achieved improvements in the
average global lower bound of 10% and 57%, respectively. Lastly, when
$\epsilon=0.03$, our method did not improve over the baseline. This again
shows that, for CNNs, configuring the bounds is less effective if the
perturbation radius is minimal.

###  5.2 Tanh-based Networks

Next, we report the results obtained for the Tanh-based CNN trained on the
CIFAR-10 dataset; these are also presented in Table 1. Notably, we achieved
consistent improvement for any given value of $\epsilon$. Furthermore, we
found that the global lower bounds are generally much lower than those
obtained for the Sigmoid-based CNN, although verifying the same properties.
Nevertheless, when $\epsilon=0.0039$, our method enables the verification of
an additional instance, which was previously unsolved. Overall, our results
demonstrate the strength of our approach, and its potential to improve
verification performance on networks with non-piecewise linear activation
functions in general.

###  5.3 Distribution of Tangent Points

To gain a better understanding of the difference between bounding functions
obtained by our method and vanilla $\mathrm{CROWN}$, we perform an empirical
analysis of tangent point distributions obtained for the linear bounding
functions of the Sigmoid-based CNN network trained on CIFAR-10 when verified
for local robustness with $\epsilon=0.0313$. Notice that this benchmark shows
the greatest improvement in the average global lower bound.

Figure 3 shows probability density functions of tangent points for lower and
upper bounding functions per activation layer of the considered network.
Notably, they show that the distribution of tangent points found by the
configured bounding algorithm is much more centred around a specific value
than those obtained by the baseline approach. Furthermore, the difference
between the empirical distribution functions of tangent points obtained by our
search method and the baseline approach was determined as statistically
significant by means of a Kolmogorov–Smirnov test with a standard significance
threshold of 0.05. Moreover, these empirical observations hint towards the
existence of an optimal region for bounding parameters of Sigmoidal activation
functions, which might be difficult to identify using the baseline approach
without automated configuration of the tangent point search method.

##  6 Conclusions and Future Work

In this work, we have shown that automated algorithm configuration can provide
a systematic and effective way for designing bounding functions of
nonlinearities beyond the commonly studied piecewise linear activation
functions (_e.g._ , ReLU). Specifically, our new method achieved consistent
improvements in average global lower bound across several robustness
verification benchmarks and perturbation radii.

At the same time, we see several fruitful avenues for future work. First of
all, the proposed search method for obtaining the tangent points is controlled
by only two hyper-parameters. A more sophisticated method, allowing for a more
fine-grained configuration of the search method, could improve performance
further. In addition, one could configure the hyper-parameters of the search
method as well as those of the $\alpha$-optimisation method, _e.g._ , the
learning rate of the projected gradient algorithm, jointly. Overall, we see
this study as a promising step towards the automated design of versatile and
efficient neural network verification algorithms.

####  6.0.1 Acknowledgment

This research was partially supported by TAILOR, a project funded by EU
Horizon 2020 research and innovation program under GA No. 952215. It further
received funding from the ERC under the European Union’s Horizon 2020 research
and innovation program (FUN2MODEL, grant agreement No. 834115) and ELSA:
European Lighthouse on Secure and Safe AI project (grant agreement No.
101070617 under UK guarantee).

## References

  * [1] Bak, S., Tran, H.D., Hobbs, K., Johnson, T.T.: Improved Geometric Path Enumeration for Verifying ReLU Neural Networks. In: Proceedings of the 32nd International Conference on Computer Aided Verification (CAV 2020). pp. 66–96 (2020)
  * [2] Bezerra, L.C., López-Ibánez, M., Stützle, T.: Automatic Component-Wise Design of Multiobjective Evolutionary Algorithms. IEEE Transactions on Evolutionary Computation 20(3), 403–417 (2015)
  * [3] Botoeva, E., Kouvaros, P., Kronqvist, J., Lomuscio, A., Misener, R.: Efficient Verification of ReLU-based Neural Networks via Dependency Analysis. In: Proceedings of the 34th AAAI Conference on Artificial Intelligence (AAAI-20). pp. 3291–3299 (2020)
  * [4] Breiman, L.: Random Forests. Machine Learning 45(1), 5–32 (2001)
  * [5] Bunel, R., Lu, J., Turkaslan, I., Torr, P.H.S., Kohli, P., Kumar, M.P.: Branch and Bound for Piecewise Linear Neural Network Verification. Journal of Machine Learning Research 21, 42:1–42:39 (2020)
  * [6] Bunel, R., Turkaslan, I., Torr, P., Kohli, P., Mudigonda, P.K.: A Unified View of Piecewise Linear Neural Network Verification. In: Advances in Neural Information Processing Systems 31 (NeurIPS 2018). pp. 1–10 (2018)
  * [7] Carlini, N., Wagner, D.: Towards Evaluating the Robustness of Neural Networks. In: Proceedings of the 38th IEEE Symposium on Security and Privacy (IEEE S&P 2017). pp. 39–57 (2017)
  * [8] Chiarandini, M., Fawcett, C., Hoos, H.H.: A Modular Multiphase Heuristic Solver for Post Enrolment Course Timetabling. In: Proceedings of the 7th International Conference on the Practice and Theory of Automated Timetabling (PATAT 2008) (2008)
  * [9] Feurer, M., Springenberg, J.T., Hutter, F.: Initializing Bayesian Hyperparameter Optimization via Meta-Learning. In: Proceedings of the 29th AAAI Conference on Artificial Intelligence (AAAI–15). pp. 1128–1135 (2015)
  * [10] Gebser, M., Kaminski, R., Kaufmann, B., Schaub, T., Schneider, M.T., Ziller, S.: A Portfolio Solver for Answer Set Programming: Preliminary Report. In: Proceedings of the 10th International Conference on Logic Programming and Nonmonotonic Reasoning (LPNMR2019). pp. 1–6 (2011)
  * [11] Gehr, T., Mirman, M., Drachsler-Cohen, D., Tsankov, P., Chaudhuri, S., Vechev, M.: AI2: Safety and Robustness Certification of Neural Networks with Abstract Interpretation. In: Proceedings of the 39th IEEE Symposium on Security and Privacy (IEEE S&P 2018). pp. 3–18 (2018)
  * [12] Goodfellow, I.J., Shlens, J., Szegedy, C.: Explaining and Harnessing Adversarial Examples. In: Proceedings of the 3rd International Conference on Learning Representations (ICLR 2015). pp. 1–11 (2015)
  * [13] Henriksen, P., Lomuscio, A.: Efficient Neural Network Verification via Adaptive Refinement and Adversarial Search. In: Proceedings of the 24th European Conference on Artificial Intelligence (ECAI 2020). pp. 2513–2520 (2020)
  * [14] Hutter, F., Babic, D., Hoos, H.H., Hu, A.J.: Boosting Verification by Automatic Tuning of Decision Procedures. In: Proceedings of Formal Methods in Computer Aided Design (FMCAD’07). pp. 27–34 (2007)
  * [15] Hutter, F., Hoos, H.H., Leyton-Brown, K.: Automated Configuration of Mixed Integer Programming Solvers. In: Proceedings of the 7th International Conference on Integration of Artificial Intelligence (AI) and Operations Research (OR) Techniques in Constraint Programming (CPAIOR 2010). pp. 186–202 (2010)
  * [16] Hutter, F., Hoos, H.H., Leyton-Brown, K.: Sequential Model-Based Optimization for General Algorithm Configuration. In: Proceedings of the 5th International Conference on Learning and Intelligent Optimization (LION 5). pp. 507–523 (2011)
  * [17] Hutter, F., Lindauer, M., Balint, A., Bayless, S., Hoos, H., Leyton-Brown, K.: the Configurable SAT Solver Challenge (CSSC). Artificial Intelligence 243, 1–25 (2017)
  * [18] Julian, K.D., Kochenderfer, M.J., Owen, M.P.: Deep neural network compression for aircraft collision avoidance systems. Journal of Guidance, Control, and Dynamics 42(3), 598–608 (2019)
  * [19] Katz, G., Barrett, C., Dill, D.L., Julian, K., Kochenderfer, M.J.: Reluplex: An Efficient SMT Solver for Verifying Deep Neural Networks. In: Proceedings of the 29th International Conference on Computer Aided Verification (CAV 2017). pp. 97–117 (2017)
  * [20] König, M., Bosman, A.W., Hoos, H.H., van Rijn, J.N.: Critically assessing the state of the art in neural network verification. Journal of Machine Learning Research 25(12), 1–53 (2024)
  * [21] König, M., Hoos, H.H., Rijn, J.N.v.: Speeding up neural network robustness verification via algorithm configuration and an optimised mixed integer linear programming solver portfolio. Machine Learning 111(12), 4565–4584 (2022)
  * [22] Kurakin, A., Goodfellow, I., Bengio, S.: Adversarial examples in the physical world. arXiv preprint arXiv:1607.02533 (2016)
  * [23] Li, L., Xie, T., Li, B.: Sok: Certified robustness for deep neural networks. In: Proceedings of the 44th IEEE Symposium on Security and Privacy (SP2023). pp. 1289–1310 (2023)
  * [24] Lopez-Ibanez, M., Stützle, T.: Automatically improving the anytime behaviour of optimisation algorithms. European Journal of Operational Research 235(3), 569–582 (2014)
  * [25] Müller, C., Serre, F., Singh, G., Püschel, M., Vechev, M.: Scaling Polyhedral Neural Network Verification on gpus. In: Proceedings of Machine Learning and Systems 3 (MLSys 2021). pp. 1–14 (2021)
  * [26] Singh, G., Ganvir, R., Püschel, M., Vechev, M.: Beyond the Single Neuron Convex Barrier for Neural Network Certification. In: Advances in Neural Information Processing Systems 32 (NeurIPS 2019). pp. 1–12 (2019)
  * [27] Singh, G., Gehr, T.: Boosting Robustness Certification of Neural networks. In: Proceedings of the 7th International Conference on Learning Representations (ICLR 2019). pp. 1–12 (2019)
  * [28] Singh, G., Gehr, T., Mirman, M., Püschel, M., Vechev, M.: Fast and Effective Robustness Certification. In: Advances in Neural Information Processing Systems 31 (NeurIPS 2018). pp. 1–12 (2018)
  * [29] Singh, G., Gehr, T., Püschel, M., Vechev, M.: An abstract domain for certifying neural networks. In: Proceedings of the 3rd ACM on Programming Languages (POPL 2019). pp. 1–30 (2019)
  * [30] Singh, G., Gehr, T., Püschel, M., Vechev, M.: An Abstract Domain for Certifying Neural Networks. In: Proceedings of the 46th ACM SIGPLAN Symposium on Principles of Programming Languages (ACMPOPL 2019). pp. 1–30 (2019)
  * [31] Szegedy, C., Zaremba, W., Sutskever, I., Bruna, J., Erhan, D., Goodfellow, I.J., Fergus, R.: Intriguing properties of neural networks. In: Proceedings of the 2nd International Conference on Learning Representations (ICLR 2014). pp. 1–10 (2014)
  * [32] Thornton, C., Hutter, F., Hoos, H.H., Leyton-Brown, K.: Auto-WEKA: Combined Selection and Hyperparameter Optimization of Classification Algorithms. In: Proceedings of the 19th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD2013). pp. 847–855 (2013)
  * [33] Tjeng, V., Xiao, K., Tedrake, R.: Evaluating Robustness of Neural Networks with Mixed Integer Programming. In: Proceedings of the 7th International Conference on Learning Representations (ICLR 2019). pp. 1–21 (2019)
  * [34] Tjeng, V., Xiao, K.Y., Tedrake, R.: Evaluating robustness of neural networks with mixed integer programming. In: 7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019. OpenReview.net (2019), https://openreview.net/forum?id=HyGIdiRqtm
  * [35] Vallati, M., Fawcett, C., Gerevini, A.E., Hoos, H., Saetti, A.: Automatic Generation of Efficient Domain-Specific Planners from Generic Parametrized Planners. In: Proceedings of the 6th Annual Symposium on Combinatorial Search (SOCS). pp. 184–192 (2013)
  * [36] Wang, S., Zhang, H., Xu, K., Lin, X., Jana, S., Hsieh, C.J., Kolter, J.Z.: Beta-crown: Efficient bound propagation with per-neuron split constraints for neural network robustness verification. In: Advances in Neural Information Processing Systems 34 (NeurIPS 2021). pp. 29909–29921 (2021)
  * [37] Xiang, W., Tran, H.D., Johnson, T.T.: Output Reachable Set Estimation and Verification for Multilayer Neural Networks. IEEE Transactions on Neural Networks and Learning Systems 29(11), 5777–5783 (2018)
  * [38] Xu, K., Zhang, H., Wang, S., Wang, Y., Jana, S., Lin, X., Hsieh, C.: Fast and complete: Enabling complete neural network verification with rapid and massively parallel incomplete verifiers. In: Proceedings of the 9th International Conference on Learning Representations (ICLR 2021) (2021)
  * [39] Zhang, H., Weng, T.W., Chen, P.Y., Hsieh, C.J., Daniel, L.: Efficient Neural Network Robustness Certification with General Activation Functions. Advances in Neural Information Processing Systems 31 (NeurIPS 2018) 31, 4944–4953 (2018)