# A Quantitative Geometric Approach to Neural Network Smoothness

Zi Wang† Gautam Prakriya‡ Somesh Jha†
†University of Wisconsin-Madison ‡The Chinese University of Hong Kong
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Fast and precise Lipschitz constant estimation of neural networks is an
important task for deep learning. Researchers have recently found an intrinsic
trade-off between the accuracy and smoothness of neural networks, so training
a network with a loose Lipschitz constant estimation imposes a strong
regularization and can hurt the model accuracy significantly. In this work, we
provide a unified theoretical framework, a quantitative geometric approach, to
address the Lipschitz constant estimation. By adopting this framework, we can
immediately obtain several theoretical results, including the computational
hardness of Lipschitz constant estimation and its approximability.
Furthermore, the quantitative geometric perspective can also provide some
insights into recent empirical observations that techniques for one norm do
not usually transfer to another one.

We also implement the algorithms induced from this quantitative geometric
approach in a tool GeoLIP. These algorithms are based on semidefinite
programming (SDP). Our empirical evaluation demonstrates that GeoLIP is more
scalable and precise than existing tools on Lipschitz constant estimation for
$\ell_{\infty}$-perturbations. Furthermore, we also show its intricate
relations with other recent SDP-based techniques, both theoretically and
empirically. We believe that this unified quantitative geometric perspective
can bring new insights and theoretical tools to the investigation of neural-
network smoothness and robustness.

##  1 Introduction

The past decade has witnessed the unprecedented success of deep learning in
many machine learning tasks (Krizhevsky et al., 2012; Mikolov et al., 2013).
Despite the growing popularity of deep learning, researchers have also found
that neural networks are very vulnerable to adversarial attacks (Szegedy et
al., 2014; Goodfellow et al., 2015; Papernot et al., 2015). As a result, it is
important to train neural networks that are robust against those attacks Madry
et al. (2019). In recent years, the deep learning community starts to focus on
certifiably robust neural networks Albarghouthi (2021); Hein and
Andriushchenko (2017); Katz et al. (2017); Cohen et al. (2019); Raghunathan et
al. (2018); Wang et al. (2022); Leino et al. (2021). One way to achieve
certified robustness is to estimate the smoothness of neural networks, where
the smoothness is measured by the Lipschitz constant of the neural network.
Recent works have found that to achieve both high accuracy and low
Lipschitzness, the network has to significantly increase the model capacity
(Bubeck and Sellke, 2021). This implies that there is an intrinsic tension
between accuracy and smoothness of neural networks.

Commonly considered adversarial attacks are the $\ell_{\infty}$ and
$\ell_{2}$-perturbations in the input space. Leino et al. (2021); Cohen et al.
(2019) have successfully trained networks with low $\ell_{2}$-Lipschitz
constant, and Huang et al. (2021) trained networks with low local
Lipschitzness for $\ell_{\infty}$-perturbations. There are no known good
techniques to train neural networks with low global Lipschitzness for
$\ell_{\infty}$-perturbations. One critical step to train such networks is to
measure the Lipschitz constant more precisely and efficiently. Jordan and
Dimakis (2021) showed that it is $\mathsf{NP}$-hard to approximate the
Lipschitz constant for $\ell_{\infty}$-perturbations within a constant factor,
and proposed an exponential-time algorithm to compute the exact Lipschitz
constant. However, researchers are interested in more scalable approaches to
certify and train the networks.

In theory, activation functions in a neural network are correlated, and the
supremum of all feasible activation patterns induces the _true_ global
Lipschitz constant. In practice, we can treat all activations independently,
and consider the maximum gradient norm over all possible activation patterns.
This upper bound of Lipschitz constant, which we call the _Formal Global
Lipschitz_ (FGL) constant, is also beneficial for certifying the robustness of
neural networks Scaman and Virmaux (2018); Latorre et al. (2020); Chen et al.
(2020); Raghunathan et al. (2018); Fazlyab et al. (2019); Pauli et al. (2022).
Throughout the paper, we will mainly consider the estimation of FGL for
$\ell_{\infty}$ and $\ell_{2}$-perturbations.

We address the estimation from the quantitative geometric perspective.
Quantitative geometry aims to understand geometric structures via quantitative
parameters, and has connections with many mathematical disciplines, such as
functional analysis and combinatorics (Naor, 2013). In computer science,
quantitative geometry plays a central role in understanding the
approximability of discrete optimization problems. We approach those hard
discrete optimization problems by considering the efficiently solvable
continuous counterparts, and analyze the precision loss due to relaxation,
which is often the SDP relaxation (Goemans and Williamson, 1995; Nesterov,
1998; Alon and Naor, 2004). By adopting the quantitative geometric approach,
we can immediately understand the computational hardness and approximability
of FGL estimation. Our algorithms on two-layer networks are the natural SDP
relaxations using the quantitative geometric perspective.

Recent works on FGL estimation for $\ell_{\infty}$-perturbations employ
polynomial optimization methods (Latorre et al., 2020; Chen et al., 2020).
Polynomial optimization is a very general framework, and many problems can be
cast in this framework (Motzkin and Straus, 1965; Goemans and Williamson,
1995). Therefore, we argue that this is not a precise characterization of the
FGL estimation problem. On the other hand, there are also several SDP-based
techniques for FGL estimations. Raghunathan et al. (2018) proposed an SDP
algorithm to estimate the FGL for $\ell_{\infty}$-perturbations of two layer
networks, and Fazlyab et al. (2019) devised an SDP algorithm to estimate the
$\ell_{2}$-FGL. We will demonstrate the intricate relationships between our
algorithms and these existing SDP relaxations on two-layer networks, and their
relationships in return inspires us to design a complete framework and
algorithms for FGL estimations on multi-layer networks.

Several empirical studies have found that techniques on one
$\ell_{p}$-perturbations often do not transfer to another ones, even though
the authors claim that _in theory_ these techniques should transfer (Fazlyab
et al., 2019; Leino et al., 2021). This in-theory claim usually comes from a
_qualitative_ perspective. In finite-dimensional space, one can always bound
one $\ell_{p}$-norm from the other one, so techniques for one
$\ell_{p}$-perturbations can also provide another bound for a different
$\ell_{p}$-perturbations. However, this bound is loose and in practice not
useful (Latorre et al., 2020; Leino et al., 2021). Instead, we believe that
when transferring techniques from one norm to another one, we should consider
the quantitative geometric principle: we separate the geometry-dependent
component from the geometry-independent one in those techniques, and modify
the geometry-dependent component accordingly for a different normed space. As
we will demonstrate, our whole work is guided by this principle. We hope that
our unified quantitative geometric approach can bring insights to the
empirical hard-to-transfer observations, and new tools to address these
issues.

#### Contributions

To summarize, we have made the following contributions:

  1. 1.

We provide a unified theoretical framework to address FGL estimations, which
immediately provides the computational hardness and SDP-based approximation
algorithms on two-layer networks (Section 3). These computational hardness and
approximability results are the first for FGL estimations.

  2. 2.

We demonstrate the relations between our algorithms and other SDP-based tools,
which also inspires us to design the algorithms for multi-layer networks
(Sections 4 and 5).

  3. 3.

We implemented the algorithms and name the tool GeoLIP. We empirically
validate our theoretical claims, and compare GeoLIP with existing tools for
estimating FGL for $\ell_{\infty}$-perturbations. The result shows that GeoLIP
provides a tighter bound than existing tools on small networks, and much
better results than naive matrix-norm-product method for deep networks, which
existing tools cannot handle (Sections 6 and 7).

##  2 Preliminaries

#### Notations

Let $[n]=\\{1,\ldots,n\\}$. For two functions $f$ and $g$, $f\circ
g(x)=f(g(x))$. A 0-1-cube is $\\{0,1\\}^{n}$, and a norm-1 cube is
$\\{-1,1\\}^{n}$ for some $n>0$. $\mathbb{R}_{+}=[0,\infty)$. For any vector
$v\in\mathbb{R}^{n}$, $\text{diag}(v)$ is an $n\times n$ diagonal matrix, with
diagonal values $v$. Let $e_{n}=(1,\ldots,1)\in\mathbb{R}^{n}$ be an
$n$-dimensional vector of all $1$’s. $||v||_{p}$ denotes the $\ell_{p}$ norm
of $v$. We use $q$ to denote the Hölder dual of $p$ as a convention, i.e.,
$\frac{1}{p}+\frac{1}{q}=1$. If $v$ is an operator in the $\ell_{p}$-space,
the operator norm of $v$ is then $||v||_{q}$. Throughout the paper, we
consider the $\ell_{p}$-norm of the input’s perturbation, and therefore, the
$\ell_{q}$-norm of the gradient, which acts as an operator on the
perturbation.

A square matrix $X\succeq 0$ means that $X$ positive semidefinite (PSD).
$\text{tr}(X)$ denotes the trace of a square matrix $X$. If
$a,b\in\mathbb{R}^{n}$, $\langle a,b\rangle$ denotes the inner product of $a$
and $b$.

#### Lipschitz Function

Given two metric spaces $(X,d_{X})$ and $(Y,d_{Y})$, a function
$f:X\rightarrow Y$ is _Lipschitz_ continuous if there exist $K>0$ such that
for all $x_{1},x_{2}\in X$,

$d_{Y}(f(x_{2}),f(x_{1}))\leq Kd_{X}(x_{2},x_{1}).$ (1)

The smallest such $K$ satisfying Equation 1, denoted by $K_{f}$, is called the
Lipschitz constant of $f$. For neural networks, $X$ is in general
$\mathbb{R}^{m}$ equipped with the $\ell_{p}$-norm. We will only consider the
case when $Y=\mathbb{R}$. In actual applications such as a classification
task, a neural network has multiple outputs. The prediction is the class with
the maximum score. One can then use the margin between each pair of class
predictions and its Lipschitz constant to certify the robustness of a given
prediction (Raghunathan et al., 2018; Leino et al., 2021). From Rademacher’s
theorem, if $f$ is Lipschitz continuous, then $f$ is is almost everywhere
differentiable, and $K_{f}=\sup_{x}||\nabla f(x)||_{q}$.

#### Neural Network as Function

A neural network $f:\mathbb{R}^{m}\rightarrow\mathbb{R}$ is characterized by
the following computation:

$f_{1}(x)=W_{1}x+b_{1};\;f_{i}(x)=W_{i}\text{diag}(\sigma(f_{i-1}(x)))+b_{i},i=2,\ldots,d.$

where $W_{i}\in\mathbb{R}^{n_{i+1}\times n_{i}}$ is the weight matrix between
the layers, $n_{1}=m$, $d$ is the depth of the neural network, $\sigma$
denotes an activation function, $b_{i}\in\mathbb{R}^{n_{i+1}}$ is the bias
term, and $f=f_{d}\circ\cdots\circ f_{1}$. Because we only consider the
$\mathbb{R}$ as the codomain of $f$, $W_{d}=u\in\mathbb{R}^{1\times n_{d}}$ is
a vector. From chain rule, the gradient of this function is

$\nabla
f(x)=W_{1}^{T}\Pi_{i=2}^{d}\text{diag}(\sigma^{\prime}(f_{i-1}(x))W_{i}^{T}.$
(2)

Common activation functions, include ReLU, sigmoid functions, are (almost)
everywhere differentiable. As a result, we are interested in the supremum
operator norm of Equation 2.

However, checking all possible inputs $x$ is infeasible, and common activation
functions have bounded derivative values, saying $[a,b]$. We are then
interested in the following value instead:

$\max_{v_{i}\in[a,b]^{n_{i}}}||W_{1}^{T}\cdot\text{diag}(v_{1})\cdot\cdots\cdot\text{diag}(v_{d})W_{d}^{T}||_{q},$
(3)

where $n_{i}$ is the dimension of each $\text{diag}(v_{i})$. We call this
value the _formal_ global Lipschitz (FGL) constant because we treat all
activation functions independent but in reality not all activation patterns
are feasible. Therefore, this is only an upper bound of the _true_ global
Lipschitz constant of the neural network. However, it is the value studied in
most global Lipschitzness literature (Scaman and Virmaux, 2018; Fazlyab et
al., 2019; Latorre et al., 2020), and also turns out useful in certifying the
robustness of neural networks (Raghunathan et al., 2018; Leino et al., 2021;
Pauli et al., 2022).

In this paper, we focus on the $\ell_{p}$-perturbation on the input, where
$p=\infty$ or $p=2$. Because $q$ is the Hölder dual of $p$, we are interested
in the value of Equation 3, when $q=1$ (for $p=\infty$), and $q=2$ (for
$p=2$).

Notice that in the ReLU-network case, $[a,b]$ is $\\{0,1\\}$, because the
derivative of ReLU is either $0$ or $1$ almost everywhere. The objective in
Equation 3 becomes an integer optimization problem. We will use ReLU-networks
as the illustration for the rest of the paper because of the popularity of
ReLU in practice and the easy presentation of the 0-1-cube. However, the
algorithms presented in this work can be adapted with minor adjustments to
other activation functions.

##  3 Two Layer Neural Networks

In this section, we will consider the single-hidden-layer neural network case.
We will show the equivalence between the FGL estimation and the matrix mixed-
norm problems. This immediately yields the complexity of FGL estimation, and
approximation algorithms for this estimation.

#### Problem Description

The two-layer neural networks can be defined in the following way:

$f(x)=u\text{diag}(y)Wx,$

where $W\in\mathbb{R}^{n\times m}$, $\text{diag}(y)\in\mathbb{R}^{n\times n}$,
and $u\in\mathbb{R}^{1\times n}$. The FGL of this network is then

$\max_{y\in\\{0,1\\}^{n}}||W^{T}\text{diag}(y)u^{T}||_{q}.$

If we expand the matrix multiplication, it is easy to check that this equals
to the following equation:

$\max_{y\in\\{0,1\\}^{n}}||W^{T}\text{diag}(u)y||_{q}.$

Let $A=W^{T}\text{diag}(u)$, then the above equation is

$\max_{y\in\\{0,1\\}^{n}}||Ay||_{q}.$ (4)

###  3.1 $\ell_{\infty}$-FGL Estimation

We will consider a natural SDP relaxation to Equation 4 when $q=1$, and
analyze the result using the celebrated _Grothendieck Inequality_ , which is a
fundamental tool in functional analysis.

#### Mixed Norm Problem

The $\infty\rightarrow 1$ mixed-norm of a matrix is defined as

$||A||_{\infty\rightarrow 1}=\max_{||x||_{\infty}=1}||Ax||_{1}.$

The mixed-norm problem appears similar to Equation 4 when $q=1$, except for
that the instead of a norm-$1$-cube, the cube in Equation 4 is a 0-1-cube.

#### SDP for the Mixed-Norm Problem

Recall that for $v\in\mathbb{R}^{m}$,
$||v||_{1}=\max_{||u||_{\infty}=1}\langle u,v\rangle$. We can reformulate the
mixed-norm problem as follows:

$\max_{x\in\\{-1,1\\}^{n}}||Ax||_{1}=\max_{(x,y)\in\\{-1,1\\}^{n+m}}\langle
Ax,y\rangle.$ (5)

If we let $z=(x^{T},y^{T})$, we can have

$\max_{(x,y)\in\\{-1,1\\}^{n+m}}\langle Ax,y\rangle=z\cdot B\cdot z^{T},$

where $B$ is a $(m+n)\times(m+n)$ matrix. The last $m$ rows and first $n$
columns of $B$ is $A$, and the rest are $0$: $B=\begin{pmatrix}0&0\\\
A&0\end{pmatrix}$.

The natural SDP relaxation of the mixed-norm problem is:

$\displaystyle\begin{split}\max\text{tr}(BX)\\\ s.t.\;X\succeq
0,X_{ii}=1,&i\in[n+m],\end{split}$ (6)

In other words, we treat $X$ as the SDP matrix relaxed from the rank-1 matrix
$z^{T}\cdot z$.

In fact, the SDP relaxation for $\infty\rightarrow 1$ mixed-norm problem
admits a constant approximation factor. Grothendieck (1956) developed the
local theory of Banach spaces, and showed that there exists an absolute value
$K_{G}$ such that following theorem holds:

######  Theorem 3.1.

For any $m,n\geq 1$, $A\in\mathbb{R}^{n\times m}$, and any Hilbert space $H$

$\max_{u_{i},v_{j}\in B(H)}\sum_{i,j}A_{ij}\langle u_{i},v_{j}\rangle_{H}\leq
K_{G}||A||_{\infty\rightarrow 1},$

where $B(H)$ denotes the unit ball of the Hilbert space. The precise value of
$K_{G}$ is still an outstanding open problem (Krivine, 1979; Braverman et al.,
2011). It is known that $K_{G}<1.783$. The approximation factor of the SDP
relaxation in Equation 6 is $K_{G}$.

#### Complexity of FGL Estimation

Alon and Naor (2004) showed that both the cut-norm and the mixed-norm of a
matrix $A$ is $\mathsf{NP}$-hard, or specifically $\mathsf{MAXSNP}$-hard, via
a reduction to the Max-Cut problem. The cut-norm of a matrix is essentially
the same as defined in Equation 5 except for that the constraint is a
0-1-cube. As a result, one can view the FGL estimation as a mixture of the
cut-norm and the mixed-norm problems. We will similarly show that the FGL
estimation problem is $\mathsf{MAXSNP}$-hard, via a cube rescaling technique
we introduce later.

######  Theorem 3.2.

FGL estimation for ReLU-networks is $\mathsf{MAXSNP}$-hard.

#### From 0-1 Cube to Norm-1 Cube

If we can transform the 0-1 cube in Equation 4 to a norm-1 cube, and formulate
an equivalent optimization problem, then one can apply the SDP program in
Equation 6 to compute an upper bound of FGL.

Let $x_{i}=(t_{i}+1)/2$, where $t_{i}\in\\{-1,1\\}$. We have

$\displaystyle\begin{split}&\max_{x\in\\{0,1\\}^{n}}||Ax||_{1}\\\
=&\max_{x\in\\{0,1\\}^{n},y\in\\{-1,1\\}^{m}}y^{T}Ax\\\
=&\max_{t\in\\{-1,1\\}^{n},y\in\\{-1,1\\}^{m}}\frac{1}{2}y^{T}A(t+e_{n})\end{split}$
(7)

Let $OPT_{1}$ be the optimal value of

$\max_{(t,y)\in\\{-1,1\\}^{n+m}}y^{T}A(t+e_{n}).$

Introduce another variable $\tau\in\\{-1,1\\}$, and let $OPT_{2}$ be the
optimal value of

$\max_{(t,y,\tau)\in\\{-1,1\\}^{n+m+1}}y^{T}A(t+\tau e_{n}).$ (8)

######  Lemma 3.3.

$OPT_{1}=OPT_{2}$.

###### Proof.

Clearly $OPT_{2}\geq OPT_{1}$.

Now if $(\hat{t},\hat{y},\tau=-1)$ is an optimal solution to Equation 8, then
$(-\hat{t},-\hat{y},\tau=1)$ is also an optimal solution, so $OPT_{2}\leq
OPT_{1}$. ∎

Now let $z=(t,\tau)$, and we can verify that $y^{T}A(t+\tau e_{n})=y^{T}Bz$,
where $B=\begin{pmatrix}A&Ae_{n}\end{pmatrix}$.

Now the semidefinite program to the FGL constant is

$\displaystyle\begin{split}\max\;&\text{tr}(BX)\\\ s.t.\;X\succeq
0,X_{ii}=1&,i\in[n+m+1],\end{split}$ (9)

where $B$ is a $(n+1+m)\times(n+1+m)$ matrix, and $B=\begin{pmatrix}0&0&0\\\
A&Ae_{n}&0\end{pmatrix}$

As a result, we have the following theorem:

######  Theorem 3.4.

There exists a polynomial-time approximation algorithm to estimate the FGL of
one-hidden-layer neural network, moreover, the approximation ratio is $K_{G}$.

###### Proof.

Let $B=\begin{pmatrix}0&0&0\\\ A&Ae_{n}&0\end{pmatrix}$. Combing Equations 7,
8 and 9, the approximation algorithm for Equation 4 where $q=1$ is induced by
the following SDP program:

$\displaystyle\begin{split}\max\;&\frac{1}{2}\text{tr}(BX)\\\ s.t.\;X\succeq
0,X_{ii}=1&,i\in[n+m+1],\end{split}$ (10)

∎

#### Proof of Theorem 3.2

###### Proof.

We will use the cube rescaling techniques introduced above. Alon and Naor
(2004) showed that matrix cut-norm is $\mathsf{MAXSNP}$-hard. We will show
that if one can solve the FGL estimation problem, then one can find the cut
norm of a matrix.

Given a matrix $A$, the cut norm of a matrix $A\in\mathbb{R}^{m\times n}$ is
defined as

$CN(A)=\max_{x\in\\{0,1\\}^{n},y\in\\{0,1\\}^{m}}\langle Ax,y\rangle.$

We need to transform $y$ from 0-1 cube to norm-1 cube, so similarly let
$y_{i}=(t_{i}+1)/2$, where $t_{i}\in\\{-1,1\\}$. The we will have

$CN(A)=\max_{x\in\\{0,1\\}^{n},y\in\\{0,1\\}^{m}}\langle
Ax,y\rangle=\frac{1}{2}\max_{x\in\\{0,1\\}^{n},t\in\\{-1,1\\}^{m}}\langle
Ax,(t+e_{m})\rangle.$

Let $B=\begin{pmatrix}A\\\ e_{m}^{T}A\end{pmatrix}$. From above we know that

$\max_{x\in\\{0,1\\}^{n},t\in\\{-1,1\\}^{m}}\langle
Ax,(t+e_{m})\rangle=\max_{x\in\\{0,1\\}^{n},(t,\tau)\in\\{-1,1\\}^{m+1}}\langle
Bx,(t,\tau)\rangle.$

One can then construct a two layer neural network, where the first weight
matrix is $B^{T}$, and the second weight matrix is
$(1,\ldots,1)\in\mathbb{R}^{n}$. Because the network we consider has only one
output, the second weight matrix is only a vector. The FGL of this network is
exactly twice of the cut-norm of $A$. ∎

###  3.2 $\ell_{2}$-FGL estimation

Scaman and Virmaux (2018) showed that the $\ell_{2}$-FGL estimation is
$\mathsf{NP}$-hard. If $q=2$ in Equation 4, the objective becomes similar to a
$\infty\rightarrow 2$ mixed-norm problem. As we shall see, this is a quadratic
optimization problem with a PSD weight matrix over a cube. It can be viewed as
a generalization of the graph Max-Cut problem, because in the quadratic
optimization formulation of Max-Cut, the matrix is the Laplacian of the graph,
a special PSD matrix (Goemans and Williamson, 1995). Nesterov (1998)
generalized Goemans-Williamson’s technique and analyzed the case when $A$ is
PSD, showing that the natural SDP relaxation in this case produces an
$\frac{\pi}{2}$-approximation ratio. This provides a
$\sqrt{\frac{\pi}{2}}$-approximation algorithm for the $\ell_{2}$-FGL
estimation.

Now let $q=2$ in Equation 4, we will have:

$\max_{y\in\\{0,1\\}^{n}}||Ay||_{2}.$

In other words, we only need to solve the following program:

$\max_{y\in\\{0,1\\}^{n}}y^{T}(A^{T}A)y.$ (11)

Let $M=A^{T}A$, then $M$ is a PSD matrix. We have demonstrated the scaling
techniques in Section 3.1. Let $x\in\\{-1,1\\}^{n+1}$, one can verify that

$\max_{y\in\\{0,1\\}^{n}}y^{T}My=\frac{1}{4}\max_{x\in\\{-1,1\\}^{n+1}}x^{T}\hat{M}x,$
(12)

where $\hat{M}=\begin{pmatrix}M&Me_{n}\\\
e_{n}^{T}M&e_{n}^{T}Me_{n}\end{pmatrix}$.

It is easy to verify that if $M$ is PSD, $\hat{M}$ is also PSD. Because
$M=A^{T}A$, $\hat{M}=(A,Ae_{n})^{T}\cdot(A,Ae_{n})$. Now we can consider the
following natural SDP relaxation to
$\max_{x\in\\{-1,1\\}^{n+1}}x^{T}\hat{M}x$:

$\displaystyle\begin{split}\max\text{tr}(\hat{M}X)\\\ s.t.\;X\succeq
0,X_{ii}=1,&i=1,\ldots,n+1.\end{split}$ (13)

Nesterov (1998) proved that this SDP relaxation admits a
$\frac{\pi}{2}$-approximation factor. This implies the following theorem:

######  Theorem 3.5.

There exists a polynomial-time approximation algorithm to estimate the
$\ell_{2}$-FGL of one-hidden-layer neural network with approximation factor
$\sqrt{\frac{\pi}{2}}$.

###### Proof.

Let $\hat{M}=\begin{pmatrix}M&Me_{n}\\\
e_{n}^{T}M&e_{n}^{T}Me_{n}\end{pmatrix}$, where $M=A^{T}A$. Combining
Equations 11, 12 and 13, the approximation algorithm for Equation 4 where
$q=2$ is induced by the following SDP program:

$\displaystyle\begin{split}\max\frac{1}{2}\sqrt{\text{tr}(\hat{M}X)}\\\
s.t.\;X\succeq 0,X_{ii}=1,&i=1,\ldots,n+1.\end{split}$ (14)

∎

######  Remark 3.6.

As we have discussed, for two layer networks, the $\ell_{p}$-FGL estimation is
essentially the $\infty\rightarrow q$ mixed-norm problem. Indeed the mixed-
norm problem is an outstanding topic in theoretical computer science. As
discussed in Bhattiprolu et al. (2018), the $\infty\rightarrow q$ mixed norm
problem has constant approximation algorithms if $q\leq 2$, and is hard to
approximate within almost polynomial factors when $q>2$. Because when $q>2$,
its Hölder dual $p<2$. This implies that for two layer networks, the FGL
estimation is much harder for $\ell_{p}$-perturbations when $p<2$.

##  4 Relations to Existing SDP Works

There are other existing SDP works on $FGL$ estimations, and here we discuss
the relationship between these works and ours. Table 1 presents existing SDP
relaxations and their scope of application.

###  4.1 $\ell_{\infty}$-FGL Estimation

Raghunathan et al. (2018) also formulated an SDP to the one-hidden-layer
neural network. If we compare their formulation with ours, theirs is
essentially

$\displaystyle\begin{split}\max\frac{1}{4}\text{tr}(CX)\\\ s.t.\;X\succeq
0,X_{ii}=1,&i\in[n+m+1],\end{split}$

where $C$ is a $(m+n+1)\times(m+n+1)$ matrix, and
$C=\begin{pmatrix}0&0&A^{T}\\\ 0&0&e_{n}^{T}A^{T}\\\ A&Ae_{n}&0\end{pmatrix}$.
This is essentially the same as ours in Equation 10. However, we provide a
rigorous derivation and simpler formulation, and also a sound theoretical
analysis of the bound, which illustrate more insights to this problem.
Raghunathan et al. (2018) treated the SDP relaxation as a heuristic for a hard
quadratic programming problem. We proved that this relaxation is not only a
heuristic, but in fact an approximation algorithm with a tight bound.

Table 1: Existing SDP relaxations for $\ell_{p}$-FGL
|  $\ell_{2}$-FGL |  $\ell_{\infty}$-FGL
---|---|---
Two-Layer Network | LipSDP-Neuron | Raghunathan et al. (2018)
Multi-Layer Network | LipSDP-Neuron | Not Available

###  4.2 $\ell_{2}$-FGL Estimation

Fazlyab et al. (2019) proposed LipSDP, another SDP algorithm for the
$\ell_{2}$-FGL estimation problem. Fazlyab et al. (2019) provided several
versions of LipSDP to balance the precision and scalability. Pauli et al.
(2022) demonstrated that the most precise version of LiPSDP, LiPSDP-Network,
fails to produce an upper bound for $\ell_{2}$-FGL. Here we will study the
less precise version, LiPSDP-Neuron, and all the subsequent references of
LiPSDP are to this version.

Surprisingly, even though the approach in LipSDP appears quite different from
the mixed-norm formulation, the empirical evaluation shows that the results
are exactly the same. We now show that LipSDP is a dual program of our
formulation for estimating $\ell_{2}$-FGL on a single-hidden-layer network. In
fact, the duality provides another evidence of LipSDP-Neuron’s correctness,
and also hints that LipSDP-Network is likely wrong, which was confirmed by
Pauli et al. (2022).

LipSDP (Fazlyab et al., 2019) for the two-layer network is

$\min_{\zeta,\lambda}\Big{\\{}\sqrt{\zeta}:\begin{pmatrix}-2abW^{T}WT-\zeta
I_{m}&(a+b)W^{T}T\\\ (a+b)TW&-2T+u^{T}u\end{pmatrix}\preceq 0,\lambda_{i}\geq
0\Big{\\}},$

where $T=\text{diag}(\lambda)$, for $\lambda\in\mathbb{R}_{+}^{n}$. This SDP
program appears very different from Equation 14, but their empirical results
are exactly the same as we will demonstrate in Section 6.

#### A New Optimization Program

We will construct a new optimization program, which we show is equivalent to
Equation 11, and LipSDP is the SDP relaxation for its dual program.

Let $\Delta x\in\mathbb{R}^{m}$ be a perturbation on $x$, and
$\Delta\sigma(x)\in\mathbb{R}^{n}$ denote the induced perturbation of $\Delta
x$ on $\text{diag}(y)$, and $w_{i}$ be the row vector of $W$. Let the input of
the $i$-th activation node on $\text{diag}(y)$ be $\gamma_{i}$. Hence,
$\gamma_{i}=w_{i}x$, and $\Delta\gamma_{i}=w_{i}\Delta x$.

The constraint from the activation function is
$\frac{\Delta\sigma(x)_{i}}{\Delta\gamma_{i}}\in[a,b]$, in other words,
$\frac{\Delta\sigma(x)_{i}}{w_{i}\Delta x}\in[a,b]$. One can write the range
constraint as

$(\Delta\sigma(x)_{i}-a\cdot w_{i}\Delta x)(\Delta\sigma(x)_{i}-b\cdot
w_{i}\Delta x)\leq 0.$

This can be written in the quadratic form:

$\begin{pmatrix}w_{i}\Delta x\\\
\Delta\sigma(x)_{i}\end{pmatrix}^{T}\begin{pmatrix}-2ab&a+b\\\
a+b&-2\end{pmatrix}\begin{pmatrix}w_{i}\Delta x\\\
\Delta\sigma(x)_{i}\end{pmatrix}\geq 0,\;\;\;\;\forall i\in[n].$ (15)

Since for the two layer network, $f(x)=u\sigma(x)$, then $\Delta
f(x)=u\Delta\sigma(x)$. The objective for $\ell_{2}$-FGL estimation is

$\max_{\Delta
x,\Delta\sigma(x)}\sqrt{\frac{(u\Delta\sigma(x))^{2}}{\langle\Delta x,\Delta
x\rangle}}.$ (16)

Notice that because $u\in\mathbb{R}^{1\times n}$, $u\Delta\sigma(x)$ is a
scalar. We can view each $y_{i}$ in Equation 11 as
$\frac{\Delta\sigma(x)_{i}}{\Delta\gamma_{i}}$, the derivative of
$\sigma(x)_{i}$ without the limit. Therefore,
$\Delta\sigma(x)_{i}=y_{i}\Delta\gamma_{i}=w_{i}y_{i}\Delta x$.
$u\Delta\sigma(x)=\Delta x\sum_{i}^{n}u_{i}y_{i}w_{i}=\Delta x(Ay)$, where
$A=W^{T}\text{diag}(u)$ as defined in Equation 4.

As a result, from Cauchy–Schwarz inequality, the above objective is

$\displaystyle\begin{split}\max_{\Delta
x,\Delta\sigma(x)}\frac{(u\Delta\sigma(x))^{2}}{\langle\Delta x,\Delta
x\rangle}&=\max_{y}(Ay)^{2},\\\ s.t.\;y\in&[a,b]^{n},\end{split}$

This demonstrates the equivalence between the new optimization program and
Equation 11 when $[a,b]=\\{0,1\\}$ for $\sigma=\text{ReLU}$.

#### Duality in LipSDP

Now we will show that LipSDP is the dual program to the program formulated
above. Let us introduce a variable $\zeta$ such that

$\zeta-\frac{(u\Delta\sigma(x))^{2}}{(\Delta x)^{2}}\geq 0.$

In other words,

$\zeta(\Delta x)^{2}-(u\Delta\sigma(x))^{2}\geq 0.$ (17)

For each constraint in Equation 15, let us introduce a dual variable
$\lambda_{i}\geq 0$. Multiply each constraint with $\lambda_{i}$, and we have

$\begin{pmatrix}\Delta x\\\
\Delta\sigma(x)_{i}\end{pmatrix}^{T}\begin{pmatrix}-2ab\lambda_{i}w_{i}^{T}w_{i}&(a+b)\lambda_{i}w_{i}^{T}\\\
(a+b)\lambda_{i}w_{i}&-2\lambda_{i}\end{pmatrix}\begin{pmatrix}\Delta x\\\
\Delta\sigma(x)_{i}\end{pmatrix}\geq 0,\;\;\;\;\forall i\in[n].$

Sum all of them, then we have

$\begin{pmatrix}\Delta x\\\
\Delta\sigma(x)\end{pmatrix}^{T}\begin{pmatrix}-2abW^{T}WT&(a+b)W^{T}T\\\
(a+b)TW&-2T\end{pmatrix}\begin{pmatrix}\Delta x\\\
\Delta\sigma(x)\end{pmatrix}\geq 0,$

where $T=\text{diag}(\lambda)$ is the $n\times n$ diagonal matrix of dual
variables.

On the other hand, Equation 17 can be rewritten in the quadratic form:

$\begin{pmatrix}\Delta x\\\
\Delta\sigma(x)\end{pmatrix}^{T}\begin{pmatrix}\zeta I_{m}&0\\\
0&-u^{T}u\end{pmatrix}\begin{pmatrix}\Delta x\\\
\Delta\sigma(x)\end{pmatrix}\geq 0.$

As a result, the dual program for the new optimization program is

$\min_{\zeta,\lambda}\Big{\\{}\sqrt{\zeta}:\begin{pmatrix}-2abW^{T}WT-\zeta
I_{m}&(a+b)W^{T}T\\\ (a+b)TW&-2T+u^{T}u\end{pmatrix}\preceq 0,\lambda_{i}\geq
0\Big{\\}},$

which is exactly LipSDP.

This dual SDP derivation is of similar form in Ben-Tal and Nemirovski (2001,
Ch.4.3.1).

##  5 Multiple Layer Neural Networks

For a multi-layer neural network, the gradient becomes a high-degree
polynomial and its $\ell_{q}$-norm estimation becomes a high-degree polynomial
optimization problem over a cube, which is in general a hard problem
(Lasserre, 2015). However, Fazlyab et al. (2019) devised a compact matrix
representation of multi-layer neural networks so that there is no need to
handle high-degree polynomials. We will show that this compact matrix
representation of neural network has a natural dual program interpretation,
which allows us to transfer the technique to the $\ell_{\infty}$-FGL
estimation.

We will first give an SDP dual program of $\ell_{\infty}$-FGL estimation on
two-layer networks. The multi-layer case follows easily. In the end, we will
briefly discuss the gradient-polynomial approach to multi-layer network FGL
estimations.

###  5.1 Dual for $\ell_{\infty}$-FGL estimation

We first formulate the $\ell_{\infty}$-FGL estimation with a similar program
as in Section 4.2, and then provide a dual program solution to it similar to
LipSDP. The difference here is that now we consider
$\ell_{\infty}$-perturbations to the input space instead of $\ell_{2}$. Hence,
the objective becomes

$\max_{\Delta x,\Delta\sigma(x)}\frac{|u\Delta\sigma(x)|}{||\Delta
x||_{\infty}}.$

If we add an extra constraint $||\Delta x||_{\infty}=1$, the above objective
becomes

$\max_{\Delta
x,\Delta\sigma(x)}\frac{1}{2}(u\Delta\sigma(x)+u\Delta\sigma(x)).$ (18)

We have $u\Delta\sigma(x)+u\Delta\sigma(x)$, because $u\Delta\sigma(x)$ is a
scalar, and we need a symmetric matrix representation in the SDP program.

The constraints are

$\begin{pmatrix}w_{i}\Delta x\\\
\Delta\sigma(x)_{i}\end{pmatrix}^{T}\begin{pmatrix}2ab&-(a+b)\\\
-(a+b)&2\end{pmatrix}\begin{pmatrix}w_{i}\Delta x\\\
\Delta\sigma(x)_{i}\end{pmatrix}\leq 0,\;\Delta(x)_{j}^{2}\leq 1,\;\forall
i\in[n],j\in[m]$

We can write $\Delta(x)_{j}^{2}\leq 1$ as

$\begin{pmatrix}1\\\ \Delta x_{j}\end{pmatrix}^{T}\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}\begin{pmatrix}1\\\ \Delta x_{j}\end{pmatrix}\geq 0$

Now let us introduce $n+m$ non-negative dual variables $(\lambda,\tau)$, where
$\lambda\in\mathbb{R}^{n}_{+}$, and $\tau\in\mathbb{R}^{m}_{+}$. If we
multiply each dual variable with the constraint and add all the constraints
together, we will have

$\begin{pmatrix}1\\\ \Delta x\\\
\Delta\sigma(x)\end{pmatrix}^{T}\begin{pmatrix}\sum_{j=1}^{m}\tau_{j}&0&0\\\
0&-2abW^{T}WT_{1}-T_{2}&(a+b)W^{T}T_{1}\\\
0&(a+b)T_{1}W&-2T_{1}\end{pmatrix}\begin{pmatrix}1\\\ \Delta x\\\
\Delta\sigma(x)\end{pmatrix}\geq 0,$

where $T_{1}=\text{diag}(\lambda)$ and $T_{2}=\text{diag}(\tau)$.

As a result, we can incorporate the objective Equation 18, and construct the
following dual program:

$\min_{\zeta,\lambda,\tau}\Big{\\{}\frac{\zeta}{2}:\begin{pmatrix}\sum_{j=1}^{m}\tau_{j}-\zeta&0&u\\\
0&-2abW^{T}WT_{1}-T_{2}&(a+b)W^{T}T_{1}\\\
u^{T}&(a+b)T_{1}W&-2T_{1}\end{pmatrix}\preceq 0,\lambda_{i},\tau_{j}\geq
0\Big{\\}}.$ (19)

###  5.2 Multi-layer Extension

We can simply extend the dual program to multiple-layer network. We vectorize
all the units in the input layer and all hidden layers, and then constrain
them using layer-wise inequalities to formulate an optimization problem.
Following the dual program solution, we will have a matrix representation of
the SDP program. Coincidentally, this dual program solution is exactly the
same as the compact neural network representation proposed by Fazlyab et al.
(2019) for $\ell_{2}$-FGL estimation. Let us use a three layer network as an
example. Consider

$f(x)=u\text{diag}(z)V\text{diag}(y)Wx,$ (20)

where $x\in\mathbb{R}^{l\times 1}$, $W\in\mathbb{R}^{n\times l}$,
$\text{diag}(y)\in\mathbb{R}^{n\times n}$, $V\in\mathbb{R}^{m\times n}$,
$\text{diag}(z)\in\mathbb{R}^{m\times m}$, $u\in\mathbb{R}^{1\times m}$.

Let $\Delta x$ denote the perturbation on the input layer, $\Delta y$ be the
perturbation on $\text{diag}(y)$ layer, and $\Delta z$ be the perturbation on
$\text{diag}(z)$ layer. Let $w_{i}$ be the row vector of $W$, and $v_{j}$ the
row vector of $V$. As shown in Section 4.2, we have the following constraints:

$\begin{pmatrix}\Delta x\\\ \Delta
y_{i}\end{pmatrix}^{T}\begin{pmatrix}-2abw_{i}^{T}w_{i}&(a+b)w_{i}^{T}\\\
(a+b)w_{i}&-2\end{pmatrix}\begin{pmatrix}\Delta x\\\ \Delta
y_{i}\end{pmatrix}\geq 0,$ (21)

$\begin{pmatrix}\Delta y\\\ \Delta
z_{j}\end{pmatrix}^{T}\begin{pmatrix}-2abv_{j}^{T}v_{j}&(a+b)v_{j}^{T}\\\
(a+b)v_{j}&-2\end{pmatrix}\begin{pmatrix}\Delta y\\\ \Delta
z_{j}\end{pmatrix}\geq 0,$ (22)

and also the box constraint on $\Delta x$:

$\begin{pmatrix}1\\\ \Delta x_{k}\end{pmatrix}^{T}\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}\begin{pmatrix}1\\\ \Delta x_{k}\end{pmatrix}\geq 0,$ (23)

for $i\in[n]$, $j\in[m]$, $k\in[l]$. The objective is:

$\max\frac{1}{2}(u\cdot\Delta z+u\cdot\Delta z).$ (24)

Following the same dual program approach, let
$\Lambda_{1}\in\mathbb{R}_{+}^{l}$ be the dual variables for the box
constraints (Equation 23), $\Lambda_{2}\in\mathbb{R}_{+}^{n}$ be the dual
variables for constraints on the first hidden layer (Equation 21),
$\Lambda_{3}\in\mathbb{R}_{+}^{m}$ be the dual variables for constraints on
the second hidden layer (Equation 22), we can formulate the following SDP
program:

$\min_{\zeta,\Lambda_{i}}\Big{\\{}\frac{\zeta}{2}:M\preceq
0,i\in[3]\Big{\\}},$

where

$M=\begin{pmatrix}\sum_{i=1}^{l}\Lambda_{1i}-\zeta&0&0&u\\\
0&-2abW^{T}WT_{2}-T_{1}&(a+b)W^{T}T_{2}&0\\\
0&(a+b)T_{2}W&-2T_{2}-2abV^{T}VT_{3}&(a+b)V^{T}T_{3}\\\
u^{T}&0&(a+b)T_{3}V&-2T_{3}\end{pmatrix},$

and $T_{1}=\text{diag}(\Lambda_{1})$, $T_{2}=\text{diag}(\Lambda_{2})$ and
$T_{3}=\text{diag}(\Lambda_{3})$. This is the SDP dual program of Equations
21, 22, 23 and 24.

For a general $d$-layer multi-layer network,

$f(x)=u\text{diag}(y_{d-1})W_{d-1}\cdots\text{diag}(y_{1})W_{1}x,$

where $x\in\mathbb{R}^{n_{1}\times 1}$, $W_{i}\in\mathbb{R}^{n_{i+1}\times
n_{i}}$, $\text{diag}(y_{i})\in\mathbb{R}^{n_{i+1}\times n_{i+1}}$, for
$i\in[d-1]$, $u\in\mathbb{R}^{1\times n_{d}}$. Let
$\Lambda_{i}\in\mathbb{R}^{n_{i}}_{+}$ for $i\in[d]$. We can then formulate
the following SDP program:

$\min_{\zeta,\Lambda_{i}}\Big{\\{}\frac{\zeta}{2}:(L+N)\preceq
0,i\in[d]\Big{\\}},$

where

$L=\begin{pmatrix}0&0&0&\dots&0&0\\\
0&-2abW_{1}^{T}W_{1}T_{2}&(a+b)W_{1}^{T}T_{2}&\dots&0&0\\\
0&(a+b)T_{2}W_{1}&-2T_{2}-2abW_{2}^{T}W_{2}T_{3}&\dots&0&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&0&\dots&(a+b)T_{d}W_{d}&-2T_{d}\\\ \end{pmatrix},$ (25)
$N=\begin{pmatrix}\sum_{i=1}^{n_{1}}\Lambda_{1i}-\zeta&0&\dots&u\\\
0&-T_{1}&\dots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
u^{T}&0&\cdots&0\end{pmatrix},$

and $T_{i}=\text{diag}(\Lambda_{i})$ for $i\in[d]$. $A$ encodes the layer-wise
constraints from each hidden unit, and $B$ captures the box constraint on each
input perturbation, and the optimization objective.

######  Remark 5.1.

If we expand the matrix inequality derived from the compact neural-network
representation , i.e., $\begin{pmatrix}A\\\
B\end{pmatrix}^{T}\begin{pmatrix}-2abT&(a+b)T\\\
(a+b)T&-2T\end{pmatrix}\begin{pmatrix}A\\\ B\end{pmatrix}$ in Fazlyab et al.
(2019, Theorem 2), where $A$, $B$ are giant sparse matrices that encode a
multi-layer network (see Fazlyab et al. (2019) for the definitions of $A$ and
$B$), we will have exactly the same matrix for network constraints as $L$
(Equation 25) in the dual program formulation. In other words, we provide a
dual optimization interpretation to the compact neural-network representation
in LipSDP.

###  5.3 Polynomial Optimization Approach to FGL estimation

We now briefly discuss the gradient approach to estimate the FGL. We still use
the three layer network of Equation 20 as an example, and the multi-layer
version follows easily.

The gradient vector of this network is

$W^{T}\text{diag}(y)V^{T}\text{diag}(z)u^{T}.$

The $i$-th component of this vector is then

$\sum_{k=1}^{m}\sum_{j=1}^{n}(u_{k}V_{kj}W_{ji})\cdot(y_{j}z_{k})$

Therefore, the $\ell_{p}$-norm of the gradient ends up being a polynomial
optimization problem over a cube. For example, the $\ell_{1}$-norm
(corresponding to $\ell_{\infty}$-perturnations) of the gradient is

$\max_{x_{i}\in\\{-1,1\\},y_{j}\in\\{0,1\\},z_{k}\in\\{0,1\\}}\sum_{i,j,k=1}^{l,n,m}T_{ijk}\cdot
x_{i}y_{j}z_{k},$ (26)

where $T_{ijk}=W_{ji}V_{kj}u_{k}$.

This is essentially a tensor cut-norm problem, and it is an open problem
whether there exists an approximation algorithm within a constant factor to
the general tensor cut-norm problem (Kannan, 2010). Notice that Equation 26 is
not a general tensor-cut-norm problem, because the tensor is generated from
the weight matrices. For example, if we fix $j$, the projected matrices of $T$
are of rank-1. Each vector in $T_{:,j,:}$ is the product of $V_{kj}u_{k}$ with
the vector $W_{j:}$:

$\forall k:T_{:,j,k}=W_{j:}V_{kj}u_{k}$

However, we do not have the technique to exploit the low-rank structure of
these special polynomial optimization problems.

##  6 Evaluation

The primary goal of our work is to provide a theoretical framework, and also
algorithms for $\ell_{\infty}$-FGL estimations on practically used networks.
The $\ell_{2}$-FGL can be computed using LipSDP (Fazlyab et al., 2019). We
have implemented the algorithms using MATLAB (MATLAB, 2021), the CVX toolbox
(CVX Research, 2020) and MOSEK solver (ApS, 2019), and name the tool GeoLIP.
To validate our theory and the applicability of our algorithms, we want to
empirically answer the following three research questions:

RQ1: How precise is GeoLIP compared to existing methods? RQ2: How scalable is
GeoLIP compared to existing methods? RQ3: Are the dual programs devised in
GeoLIP valid?

#### Experimental Design

To answer the first two questions, we will run GeoLIP and existing tools that
measure $\ell_{\infty}$-FGL on various feed-forward neural networks trained
with the MNIST dataset (LeCun and Cortes, 2010). The tools and network
configurations will be specified. We will record the computed
$\ell_{\infty}$-FGL to compare the precision, and the computation time to
compare the scalability.

To answer RQ3, we will run the natural SDP relaxations for $\ell_{q}$-FGL
estimations proposed in Section 3, LipSDP for $\ell_{2}$-FGL estimation, and
the dual program devised in Section 5.1 for $\ell_{\infty}$-FGL on two-layer
neural networks, and compare their computed FGL.

#### SETUP

Our main baseline tool is LiPopt (Latorre et al., 2020), which is an
$\ell_{\infty}$-FGL estimation tool.111Another method is proposed by Chen et
al. (2020), however, the code is not available and we are not able to compare
it with GeoLIP. Notice that LiPopt is based on the Python Gurobi solver
(Gurobi Optimization, LLC, 2022), while we use a MATLAB CVX and MOSEK solver.

_MP_ stands for the weight-matrix norm product method. This is the naive upper
bound of FGL.

_BruF_ stands for an brute-force search for the maximum activation pattern,
via an exhaustive enumeration of all possible activation patterns. This is the
ground truth for FGL estimations, however, because this search take
exponential time, we are only able to run it on networks with a few hidden
units.

_Sample_ means that we randomly sample $200,000$ points in the input space and
compute the gradient norm at those points. Notice that this is a lower bound
of the true Lipschitz constant, and so a lower bound of FGL.

All the experiments are run on a workstation with forty-eight Intel® Xeon®
Silver 4214 CPUs running at 2.20GHz, and 258 GB of memory. We use N/A in the
result tables to denote the computation takes too much time ($>10$ hours).

#### Network Setting

We run the experiments on fully-connected feed-forward neural networks,
trained with the MNIST dataset for $10$ epochs using the ADAM optimizer
(Kingma and Ba, 2017). All the trained networks have accuracy greater than
$92\%$ on the test data. For two-layer networks, we use $8$, $16$, $64$,
$128$, $256$ hidden nodes. For multiple-layer networks, we consider $3$, $7$,
$8$-layer networks, and each hidden layer has $64$ ReLU units.

Because MNIST has $10$ classes, we will report the estimated FGL with respect
to label $8$ as in Latorre et al. (2020), and the average running time for the
computation: we record the total computation time for all 10 classes from each
tool, and report average time per class.

LiPopt relies on a linear programming (LP) hierarchy for the polynomial
optimization problem. We will use LiPopt-k to denote the $k$-th degree of the
LP hierarchy. As Latorre et al. (2020) pointed out, the degree has to be at
least the depth of the network to compute a valid bound, so we have to use at
least LiPopt-k for $k$-layer networks.

#### Single Hidden Layer

We consider two layer neural networks with different numbers of hidden units.
The results are summarized in Tables 2 and 3.

Table 2: $\ell_{\infty}$-FGL Estimation of Different Methods for Two Layer
Network
#Units | DGeoLIP | NGeoLIP | LiPopt-2 | MP | Sample | BruF
---|---|---|---|---|---|---
8 | 142.19 | 142.19 | 180.38 | 411.90 | 134.76 | 134.76
16 | 185.18 | 185.18 | 259.44 | 578.54 | 175.24 | 175.24
64 | 287.60 | 287.60 | 510.00 | 1207.70 | 253.89 | N/A
128 | 346.27 | 346.27 | 780.46 | 2004.34 | 266.22 | N/A
256 | 425.04 | 425.04 | 1011.65 | 2697.38 | 306.98 | N/A
Table 3: Average Running Time of Different Methods for Two Layer Network
$\ell_{\infty}$-FGL Estimation (in seconds)
# Hidden Units | DGeoLIP | NGeoLIP | LiPopt-2
---|---|---|---
8 | 23.1 | 21.5 | 1533
16 | 28.1 | 22.3 | 1572
64 | 93.4 | 31.7 | 1831
128 | 292.5 | 42.2 | 2055
256 | 976.0 | 70.9 | 2690

#### Multiple Hidden Layers

We consider $3$, $7$, $8$-layer neural networks. Each hidden layer in the
network has $64$ ReLU units. The results are summarized in Tables 4 and 5.

Table 4: $\ell_{\infty}$-FGL Estimation of Different Methods for Multi-layer
Network
# Layers | GeoLIP | Matrix Product | Sample | liptopt
---|---|---|---|---
3 | 529.42 | 9023.65 | 311.88 | N/A
7 | 5156.5 | $1.423*10^{7}$ | 1168.8 | N/A
8 | 8327.2 | $8.237*10^{7}$ | 1130.6 | N/A
Table 5: Average Running Time of GeoLIP for Multi-layer Network
$\ell_{\infty}$-FGL Estimation (in seconds)
3-Layer Net | 7-Layer Net | 8-Layer Net
---|---|---
120.9 | 284.3 | 329.5

#### Duality

We measure the $\ell_{2}$ and $\ell_{\infty}$-FGL on two-layer networks. The
$\ell_{2}$-FGL is measured by LipSDP and the natural SDP relaxation as in
Equation 14, which we call NGeoLIP. The $\ell_{2}$-FGL is measured by the
natural SDP relaxation as in Equation 10, which is NGeoLIP, and the dual form
as in Equation 19, which we call DGeoLIP. Still, we consider networks with
$8$, $16$, $64$, $128$, $256$ hidden nodes. The results are summarized in
Tables 2 and 6. In addition, we also record their running time in Tables 3 and
7.

Table 6: $\ell_{2}$-FGL Estimation of Different Methods for Two Layer Network
#Units | NGeoLIP | LipSDP | MP | Sample | BruF
---|---|---|---|---|---
8 | 6.531 | 6.531 | 11.035 | 6.527 | 6.527
16 | 8.801 | 8.801 | 13.936 | 8.795 | 8.799
64 | 12.573 | 12.573 | 22.501 | 11.901 | N/A
128 | 15.205 | 15.205 | 30.972 | 13.030 | N/A
256 | 18.590 | 18.590 | 35.716 | 14.610 | N/A
Table 7: Average Running Time of LipSDP and NGeoLIP for Two Layer Network
$\ell_{2}$-FGL Estimation (in seconds)
# Hidden Units | DGeoLIP | NGeoLIP
---|---|---
8 | 11.5 | 1.2
16 | 15.7 | 1.2
64 | 64.2 | 1.3
128 | 216.1 | 1.7
256 | 758.1 | 4.1

##  7 Discussion

In the following, we will discuss how our evaluation addresses the research
questions, and their implications. We first address scalability, and then
precision. In the end, we discuss duality and the related quantitative
geometric perspective .

###  7.1 Scalability

#### Scalability of GeoLIP

In the experiments of LiPopt, we only used LiPopt-2. In theory, if one can go
higher in the LP hierarchy in LiPopt, the result becomes more precise.
However, in the case of fully connected neural networks, using degree-3 in
LiPopt is already impractical. For example, for the simplest network we use,
i.e., the single hidden layer neural network with $8$ hidden units, using
LiPopt-3, one $\ell_{\infty}$-FGL computation needs at least $200$ hours
projected by LiPopt. As a result, for all the LiPopt-related experiments, we
are only able to run LiPopt-2. LiPopt is unable to handle neural networks with
more than two layers because this requires LiPopt with degree beyond 2.

Even if we only consider LiPopt-2 on two layer networks, the running time is
still much higher compared to GeoLIP. This demonstrates the great advantage of
GeoLIP in terms of scalability compared with LiPopt.

#### Sparsity might be a bug, not a feature

Latorre et al. (2020) argued that LiPopt can exploit the sparsity of neural
network, and works better on architectures with higher sparsity. Empirically,
the authors examined LiPopt on networks that are trained to encourage
sparsity. However, improving sparsity of a neural network may hurt its
robustness. Bubeck and Sellke (2021) showed that for a network to achieve both
high accuracy and smoothness, the neural network has to increase its capacity
significantly. Increasing sparsity effectively decreases the number of
neurons, and therefore may hurt the model’s robustness. A thorough
investigation of sparsity is necessary when it comes to the robustness of
neural networks, and this can be a future work.

###  7.2 Precison

If we compare LiPopt-2 with GeoLIP on two layer networks from Table 2, it is
clear that GeoLIP produces more precise results. For networks with depth
greater than $2$, we can only compare GeoLIP with the naive matrix-norm-
product method. As we can see from all experiments, GeoLIP’s estimation on FGL
is always much lower than the matrix-norm product, for example, for the
$8$-layer network, the FGL measured by the matrix-norm product is around
$10000$ times greater than GeoLIP’s result. Given that GeoLIP is the only tool
that can scale to deep networks other than the matrix-norm product method,
this demonstrates the importance of GeoLIP.

We have also shown that the two-layer network $\ell_{\infty}$-FGL estimation
from GeoLIP has a theoretical guarantee with an approximation factor
$K_{G}<1.783$. If we compare the results from GeoLIP and Sampling in Table 2,
which is a lower bound of true Lipschitz constant, the ratio is still within
$1.783$. This validates our theoretical claim.

More interestingly, we showed that GeoLIP’s approximation factor for
$\ell_{2}$-FGL estimation on two layer networks is
$\sqrt{\frac{\pi}{2}}\approx 1.253$. The $\ell_{2}$-FGL is very close to the
sampled lower bound of true Lipschitz constant in Table 6. On the other hand,
because the result from GroLIP is an upper bound of FGL, and this result is
not much greater than the sampled lower bound of true Lipschitz constant, this
empirically demonstrates that the true Lipschitz constant is not very
different from the FGL on two-layer networks.

###  7.3 Duality

In Section 4.2, we have demonstrated the duality between NGeoLIP and LipSDP
for $\ell_{2}$-FGL estimation on two layer networks, even though the
approaches appears drastically different. The experiments show that the
results of LipSDP and GeoLIP on two-layer-network $\ell_{2}$-FGL estimation
are exactly the same, which empirically demonstrates the duality between
LipSDP and GeoLIP. Similarly, the dual program proposed in Section 5.1 also
matches the natural SDP relaxation of $\ell_{\infty}$-FGL estimation on two-
layer networks. These results also provide empirical evidence for the
correctness of all the SDP programs. Though Pauli et al. (2022) showed that
the most precise version of LipSDP is invalid for estimating an upper bound of
$\ell_{2}$-FGL, our dual-program argument shows that the less precise version
of LipSDP is correct.

If we compare the running time in Tables 3 and 7, the dual program takes more
time to solve than the natural relaxation. This is particularly true when the
number of hidden neurons increases. From the reported numbers of variables and
equality constraints by CVX, the dual program and natural relaxation have
similar numbers. It is also observed that the CPU usage is higher when the
natural relaxation is being solved. We want to point out that the running time
and optimization algorithm are solver-dependent, and efficiently solving SDP
is beyond the scope of this work. It is an interesting future direction to
exploit the block structure of the dual programs, and develop algorithms that
are compatible with those programs, because training low Lipschitz network is
a critical task, and it is promising to incorporate the SDP programs.

####  $\ell_{2}$ versus $\ell_{\infty}$ FGLs

If we compare results from Tables 2 and 6, we can also find that the
discrepancy between matrix product method and sampled lower bound is much
smaller in the $\ell_{2}$ case. This could also explain why Gloro works for
$\ell_{2}$-perturbations but not the $\ell_{\infty}$ case in practice, where
Leino et al. (2021) used matrix-norm product to compute the Lipschitz constant
of the network in Gloro.

#### SDP Relaxation

Applying SDP on intractable combinatorial optimization problem was pioneered
by the seminal Goemans-Williamson algorithm for the Max-Cut problem (Goemans
and Williamson, 1995). Our $\ell_{2}$-FGL estimation algorithm can be viewed
as a direct extension from Goemans and Williamson (1995). The
$\sqrt{\pi/2}$-approximation factor is from Nesterov (1998). The only
requirement is that the weight matrix in the quadratic polynomial is PSD.
Since for Max-Cut, the weight matrix is the graph Laplacian matrix, a special
PSD matrix, so the analysis produces a tighter approximation factor than
$\pi/2$. Under the unique games conjecture (Khot, 2002), the approximation
factor in Goemans-Williamson algorithm is optimal. On one hand, our dual-
program interpretation of LipSDP-Neuron and the empirical evaluation provide
evidence that LipSDP-Neuron is correct. On the other hand, a more precise
estimation of $\ell_{2}$-FGL than the natural SDP relaxation might also
improve the Max-Cut’s algorithm, which is a very hard task. This hints that
LipSDP-Network might be incorrect, and indeed Pauli et al. (2022) demonstrated
the flaw in LipSDP-Network.

#### Limitation of Polynomial Optimization

The polynomial optimization method is a very general framework, and therefore
solving it is also hard (Lasserre, 2015). For example, finding the maximum
independent set of a graph can be encoded as a quadratic polynomial
optimization problem with variables on a cube (Motzkin and Straus, 1965). We
argue that approaching FGL-estimation from the perspective of polynomial
optimization loses the accurate characterization of this problem. For example,
for two-layer network, we have provided constant approximation algorithms to
estimate FGL in both $\ell_{\infty}$ and $\ell_{2}$ cases. However, for a
general polynomial optimization problem on a cube, we cannot achieve constant
approximation, because the maximum independent set problem cannot be
approximated to a constant factor in polynomial time unless
$\mathsf{P}=\mathsf{NP}$ (Trevisan, 2004).

#### Quantitative Geometric Perspective on Smoothness

Recent robustness studies find that techniques for one perturbations cannot
transfer to other perturbations. For example, Fazlyab et al. (2019) proposed
an SDP for the $\ell_{2}$-FGL estimation, and argue that this can generalize
to other $\ell_{p}$-FGL estimation using Hölder’s inequality. It turns out
this scaling is too loose (Latorre et al., 2020). We argue that this is a
qualitative perspective. From Hölder’s inequality, in the finite-dimensional
space, one can always bound one $\ell_{p}$ norm from another one. All the
$\ell_{p}$-spaces appear similar qualitatively, because once we have one
$\ell_{p}$-norm, we can derive an upper bound of a different $\ell_{p}$-norm,
but quantitatively, this bound is too loose and not useful in practice.

In this work, we adopt the the quantitative geometric perspective for
measuring the smoothness. For $\ell_{2}$-FGL, the perturbation set is a ball;
and for $\ell_{\infty}$-FGL, the perturbation set is a cube. The high
dimensional quantitative property of cubes and balls are quite counter-
intuitive compared to the low dimensional world, so we need to rigorously
incorporate the geometry to the techniques when considering different norms.

When we study the FGL estimation for two-layer networks, we consider the
precise mixed-norm problems for different perturbations. Similarly, when we
design the optimization program, instead of bounding the operator’s $\ell_{p}$
norm from another one, we encode the geometric constraint from the
perturbation set to the SDP. In fact, the dual optimization program is only
constrained by the underlying perturbation geometry and the layer-wise
restriction from each hidden unit, so the constraints and objective exactly
encode the FGL-estimation problem without additional assumptions. Because
often the SDP relaxation for intractable problems give the optimal known
algorithms, we conjecture that GeoLIP and LipSDP are (nearly) optimal for
$\ell_{\infty}$ and $\ell_{2}$-FGL estimations. More specifically,

Conjecture: If $\ell_{\infty}$ and $\ell_{2}$-FGL estimation problems have
constant-factor approximation algorithms (dependent or independent of network
depth), then GeoLIP and LipSDP also provide constant approximation guarantee.

One might be able to study the approximability of FGL estimations from the
polynomial optimization problems described in Section 5.3.

##  8 Related Work

We have discussed extensively recent works on FGL estimations (Fazlyab et al.,
2019; Latorre et al., 2020; Chen et al., 2020; Pauli et al., 2022; Raghunathan
et al., 2018), and compare GeoLIP with them. In particular, Raghunathan et al.
(2018) proposed to incorporate their SDP to train certifiably robust neural
networks via the dual program. However, their SDP formulation only works for
two-layer networks. Because the SDP programs in our work are already in the
dual form, and work for multi-layer networks, it is an immediate future work
to train general certifiably robust networks.

The global Lipschitz constant of a neural network is data-independent. In
practice, the Lipschitz constant restricted on a specific dataset might have a
lower Lipschitz constant, and further improves the certifiable robustness on
the dataset. Latorre et al. (2020) also proposed to use LiPopt to estimate the
local Lipschitz constant. Estimating the local Lipschitz constant is not the
problem studied in our work, and there are tools specifically designed for
local perturbations and the Lipschitz constant (Laurel et al., 2022; Zhang et
al., 2019).

##  9 Conclusion

In this work, we have provided a quantitative geometric framework for FGL
estimations, and also algorithms for the $\ell_{\infty}$-FGL estimation. We
hope that this framework can bring new insights and theoretical tools to the
study of neural-network smoothness and robustness. One important lesson is
that when transferring techniques from one perturbations to another ones, we
should also transfer the underlying geometry. In the meantime, we believe that
the algorithms from our work can be a crucial step to train smooth neural
networks.

## References

  * (1)
  * Albarghouthi (2021) Aws Albarghouthi. 2021\.  Introduction to Neural Network Verification.  arXiv:2109.10317 [cs.LG]
  * Alon and Naor (2004) Noga Alon and Assaf Naor. 2004.  Approximating the Cut-Norm via Grothendieck’s Inequality. In _Proceedings of the Thirty-Sixth Annual ACM Symposium on Theory of Computing_ (Chicago, IL, USA) _(STOC ’04)_. Association for Computing Machinery, New York, NY, USA, 72–80.  https://doi.org/10.1145/1007352.1007371
  * ApS (2019) MOSEK ApS. 2019\.  _The MOSEK optimization toolbox for MATLAB manual. Version 9.0._ http://docs.mosek.com/9.0/toolbox/index.html
  * Ben-Tal and Nemirovski (2001) Aharon Ben-Tal and Arkadi Nemirovski. 2001.  _Lectures on Modern Convex Optimization_.  Society for Industrial and Applied Mathematics.  https://doi.org/10.1137/1.9780898718829 arXiv:https://epubs.siam.org/doi/pdf/10.1137/1.9780898718829
  * Bhattiprolu et al. (2018) Vijay Bhattiprolu, Mrinalkanti Ghosh, Venkatesan Guruswami, Euiwoong Lee, and Madhur Tulsiani. 2018\.  Inapproximability of Matrix $p\rightarrow q$ Norms.  arXiv:1802.07425 [cs.CC]
  * Braverman et al. (2011) Mark Braverman, Konstantin Makarychev, Yury Makarychev, and Assaf Naor. 2011.  The Grothendieck Constant is Strictly Smaller than Krivine’s Bound. In _2011 IEEE 52nd Annual Symposium on Foundations of Computer Science_. 453–462.  https://doi.org/10.1109/FOCS.2011.77
  * Bubeck and Sellke (2021) Sebastien Bubeck and Mark Sellke. 2021.  A Universal Law of Robustness via Isoperimetry. In _Advances in Neural Information Processing Systems_ , A. Beygelzimer, Y. Dauphin, P. Liang, and J. Wortman Vaughan (Eds.).  https://openreview.net/forum?id=z71OSKqTFh7
  * Chen et al. (2020) Tong Chen, Jean-Bernard Lasserre, Victor Magron, and Edouard Pauwels. 2020.  Semialgebraic Optimization for Lipschitz Constants of ReLU Networks.  arXiv:2002.03657 [math.OC]
  * Cohen et al. (2019) Jeremy M Cohen, Elan Rosenfeld, and J. Zico Kolter. 2019\.  Certified Adversarial Robustness via Randomized Smoothing.  arXiv:1902.02918 [cs.LG]
  * CVX Research (2020) Inc. CVX Research. 2020\.  CVX: Software for Disciplined Convex Programming, Version 2.2, Build 1148.  http://cvxr.com/cvx.
  * Fazlyab et al. (2019) Mahyar Fazlyab, Alexander Robey, Hamed Hassani, Manfred Morari, and George Pappas. 2019.  Efficient and Accurate Estimation of Lipschitz Constants for Deep Neural Networks. In _Advances in Neural Information Processing Systems_ , H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Eds.), Vol. 32. Curran Associates, Inc.  https://proceedings.neurips.cc/paper/2019/file/95e1533eb1b20a97777749fb94fdb944-Paper.pdf
  * Goemans and Williamson (1995) Michel X. Goemans and David P. Williamson. 1995.  Improved Approximation Algorithms for Maximum Cut and Satisfiability Problems Using Semidefinite Programming.  _J. ACM_ 42, 6 (nov 1995), 1115–1145.  https://doi.org/10.1145/227683.227684
  * Goodfellow et al. (2015) Ian J. Goodfellow, Jonathon Shlens, and Christian Szegedy. 2015.  Explaining and Harnessing Adversarial Examples.  arXiv:1412.6572 [stat.ML]
  * Grothendieck (1956) A. Grothendieck. 1956\.  Résumé de la théorie métrique des produits tensoriels topologiques.  Bol. Soc. Mat. São Paulo 8, 1-79 (1956)..
  * Gurobi Optimization, LLC (2022) Gurobi Optimization, LLC. 2022\.  Gurobi Optimizer Reference Manual.  https://www.gurobi.com
  * Hein and Andriushchenko (2017) Matthias Hein and Maksym Andriushchenko. 2017.  Formal Guarantees on the Robustness of a Classifier against Adversarial Manipulation. In _Proceedings of the 31st International Conference on Neural Information Processing Systems_ (Long Beach, California, USA) _(NIPS’17)_. Curran Associates Inc., Red Hook, NY, USA, 2263–2273.
  * Huang et al. (2021) Yujia Huang, Huan Zhang, Yuanyuan Shi, J Zico Kolter, and Anima Anandkumar. 2021.  Training Certifiably Robust Neural Networks with Efficient Local Lipschitz Bounds.  _NeurIPS_ (2021).
  * Jordan and Dimakis (2021) Matt Jordan and Alexandros G. Dimakis. 2021.  Exactly Computing the Local Lipschitz Constant of ReLU Networks.  arXiv:2003.01219 [stat.ML]
  * Kannan (2010) Ravindran Kannan. 2010\.  Spectral Methods for Matrices and Tensors. In _Proceedings of the Forty-Second ACM Symposium on Theory of Computing_ (Cambridge, Massachusetts, USA) _(STOC ’10)_. Association for Computing Machinery, New York, NY, USA, 1–12.  https://doi.org/10.1145/1806689.1806691
  * Katz et al. (2017) Guy Katz, Clark Barrett, David L. Dill, Kyle Julian, and Mykel J. Kochenderfer. 2017.  Reluplex: An Efficient SMT Solver for Verifying Deep Neural Networks. In _Computer Aided Verification_ , Rupak Majumdar and Viktor Kunčak (Eds.). Springer International Publishing, Cham, 97–117.
  * Khot (2002) Subhash Khot. 2002\.  On the Power of Unique 2-Prover 1-Round Games. In _Proceedings of the Thiry-Fourth Annual ACM Symposium on Theory of Computing_ (Montreal, Quebec, Canada) _(STOC ’02)_. Association for Computing Machinery, New York, NY, USA, 767–775.  https://doi.org/10.1145/509907.510017
  * Kingma and Ba (2017) Diederik P. Kingma and Jimmy Ba. 2017.  Adam: A Method for Stochastic Optimization.  arXiv:1412.6980 [cs.LG]
  * Krivine (1979) J.L Krivine. 1979\.  Constantes de Grothendieck et fonctions de type positif sur les sphères.  _Advances in Mathematics_ 31, 1 (1979), 16–30.  https://doi.org/10.1016/0001-8708(79)90017-3
  * Krizhevsky et al. (2012) Alex Krizhevsky, Ilya Sutskever, and Geoffrey E. Hinton. 2012\.  Imagenet classification with deep convolutional neural networks. In _Advances in Neural Information Processing Systems_.
  * Lasserre (2015) Jean Bernard Lasserre. 2015\.  _An Introduction to Polynomial and Semi-Algebraic Optimization_.  Cambridge University Press.  https://doi.org/10.1017/CBO9781107447226
  * Latorre et al. (2020) Fabian Latorre, Paul Rolland, and Volkan Cevher. 2020\.  Lipschitz constant estimation of Neural Networks via sparse polynomial optimization. In _International Conference on Learning Representations_.  https://openreview.net/forum?id=rJe4_xSFDB
  * Laurel et al. (2022) Jacob Laurel, Rem Yang, Gagandeep Singh, and Sasa Misailovic. 2022\.  A Dual Number Abstraction for Static Analysis of Clarke Jacobians.  _Proc. ACM Program. Lang._ 6, POPL, Article 56 (jan 2022), 30 pages.  https://doi.org/10.1145/3498718
  * LeCun and Cortes (2010) Yann LeCun and Corinna Cortes. 2010.  MNIST handwritten digit database.  http://yann.lecun.com/exdb/mnist/.  (2010).  http://yann.lecun.com/exdb/mnist/
  * Leino et al. (2021) Klas Leino, Zifan Wang, and Matt Fredrikson. 2021.  Globally-Robust Neural Networks. In _International Conference on Machine Learning (ICML)_.
  * Madry et al. (2019) Aleksander Madry, Aleksandar Makelov, Ludwig Schmidt, Dimitris Tsipras, and Adrian Vladu. 2019\.  Towards Deep Learning Models Resistant to Adversarial Attacks.  arXiv:1706.06083 [stat.ML]
  * MATLAB (2021) MATLAB. 2021.  _9.11.0.1837725 (R2021b) Update 2_.  The MathWorks Inc., Natick, Massachusetts.
  * Mikolov et al. (2013) Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg Corrado, and Jeffrey Dean. 2013.  Distributed Representations of Words and Phrases and Their Compositionality. In _Proceedings of the 26th International Conference on Neural Information Processing Systems - Volume 2_ (Lake Tahoe, Nevada) _(NIPS’13)_. Curran Associates Inc., Red Hook, NY, USA, 3111–3119.
  * Motzkin and Straus (1965) T. S. Motzkin and E. G. Straus. 1965.  Maxima for Graphs and a New Proof of a Theorem of Turán.  _Canadian Journal of Mathematics_ 17 (1965), 533–540.  https://doi.org/10.4153/CJM-1965-053-6
  * Naor (2013) Assaf Naor. 2013\.  Quantitative geometry.  _Proceedings of the National Academy of Sciences_ 110, 48 (2013), 19202–19205.  https://doi.org/10.1073/pnas.1320388110 arXiv:https://www.pnas.org/content/110/48/19202.full.pdf
  * Nesterov (1998) Yu Nesterov. 1998\.  Semidefinite relaxation and nonconvex quadratic optimization.  _Optimization Methods and Software_ 9, 1-3 (1998), 141–160.  https://doi.org/10.1080/10556789808805690 arXiv:https://doi.org/10.1080/10556789808805690
  * Papernot et al. (2015) Nicolas Papernot, Patrick McDaniel, Somesh Jha, Matt Fredrikson, Z. Berkay Celik, and Ananthram Swami. 2015\.  The Limitations of Deep Learning in Adversarial Settings.  arXiv:1511.07528 [cs.CR]
  * Pauli et al. (2022) Patricia Pauli, Anne Koch, Julian Berberich, Paul Kohler, and Frank Allgöwer. 2022.  Training Robust Neural Networks Using Lipschitz Bounds.  _IEEE Control Systems Letters_ 6 (2022), 121–126.  https://doi.org/10.1109/LCSYS.2021.3050444
  * Raghunathan et al. (2018) Aditi Raghunathan, Jacob Steinhardt, and Percy Liang. 2018\.  Certified Defenses against Adversarial Examples. In _International Conference on Learning Representations_.  https://openreview.net/forum?id=Bys4ob-Rb
  * Scaman and Virmaux (2018) Kevin Scaman and Aladin Virmaux. 2018.  Lipschitz Regularity of Deep Neural Networks: Analysis and Efficient Estimation. In _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ (Montréal, Canada) _(NIPS’18)_. Curran Associates Inc., Red Hook, NY, USA, 3839–3848.
  * Szegedy et al. (2014) Christian Szegedy, Wojciech Zaremba, Ilya Sutskever, Joan Bruna, Dumitru Erhan, Ian Goodfellow, and Rob Fergus. 2014\.  Intriguing properties of neural networks.  arXiv:1312.6199 [cs.CV]
  * Trevisan (2004) Luca Trevisan. 2004\.  Inapproximability of Combinatorial Optimization Problems.  arXiv:cs/0409043 [cs.CC]
  * Wang et al. (2022) Zi Wang, Aws Albarghouthi, Gautam Prakriya, and Somesh Jha. 2022.  Interval Universal Approximation for Neural Networks.  _Proc. ACM Program. Lang._ 6, POPL, Article 14 (jan 2022), 29 pages.  https://doi.org/10.1145/3498675
  * Zhang et al. (2019) Huan Zhang, Pengchuan Zhang, and Cho-Jui Hsieh. 2019\.  RecurJac: An Efficient Recursive Algorithm for Bounding Jacobian Matrix of Neural Networks and Its Applications.  arXiv:1810.11783 [stat.ML]