Semialgebraic Representation of Monotone Deep Equilibrium Models and Applications to Certification
Year: 2022

Abstract:
Deep equilibrium models are based on implicitly defined functional relations and have shown competitive performance compared with the traditional deep networks. Monotone operator equilibrium networks (monDEQ) retain interesting performance with additional theoretical guaranties. Existing certification tools for classical deep networks cannot directly be applied to monDEQs for which much fewer tools exist. We introduce a semialgebraic representation for ReLU based monDEQs which allows to approximate the corresponding input output relation by semidefinite programming (SDP). We present several applications to network certification and obtain SDP models for the following problems : robustness certification, Lipschitz constant estimation, ellipsoidal uncertainty propagation. We use these models to certify robustness of monDEQs w.r.t. a generalLq norm. Experimental results show that the proposed models outperform existing approaches for monDEQ certification. Furthermore, our investigations suggest that monDEQs are much more robust to L2 perturbations than L∞ perturbations.

heading: 1 Introduction
text: With the increasing success of Deep Neural Networks (DNN) (e.g. computer vision, natural language processing), one witnesses a significant increase in size and complexity (topology and activation functions). This generates difficulties for theoretical analysis and a posteriori performance evaluation. This is problematic for applications where robustness issues are crucial, for example inverse problems (IP) in scientific computing. Indeed such IPs are notoriously ill-posed and as stressed in the March 2021 issue of SIAM News [1], “Yet DL has an Achilles’ heel. Current implementations can be highly unstable, meaning that a certain small perturbation to the input of a trained neural network can cause substantial change in its output. This phenomenon is both a nuisance and a major concern for the safety and robustness of DL-based systems in critical applications—like healthcare—where reliable computations are essential". Indeed, the Instability Theorem [1] predicts unavoidable lower bound on Lipschitz contants, which may explain the lack of stability of some DNNs, over-performing on training sets. This underlines the need to evaluate precisely a posteriori critical indicators, such as Lipschitz constants of DNNs. However, obtaining an accurate upper bounds on the Lipschitz constant of a DNN is a hard problem, it reduces to prove a global inequality “Ψ(x) ≥ 0 for all x in a domain”,
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
i.e., to provide a proof of positivity for a function Ψ, which has no simple explicit expression. Even for modest size DNNs this task is practically challenging, previous successful attempts [8, 28] were restricted in practice to no more than two hidden layers with less than a hundred nodes. More broadly existing attempts to DNN certification rely either on zonotope calculus, linear programming (LP) or hierarchies of SDP based on positivity certificates from algebraic geometry [34], which may suffer from the curse of dimensionality.
Recently, implicit deep learning [12] arises as a generalization of the recursive rules of traditional feedforward neural networks. Specifically, Deep Equilibrium Models (DEQ) [2] have emerged as a potential alternative to classical DNNs. With their much simpler layer structure, they provide competitive results on machine learning benchmarks [2, 3]. Training DEQs involves solving fix-point equations, which requires further conditions to make the iteration converge. Fortunately, Monotone operator equilibrium network (monDEQ) introduced in [45] satisfies such conditions. Moreover, the authors in [33] provide explicit bounds on global Lipschitz constant of monDEQs (w.r.t. the L2-norm) which can be used for robustness certification.
From a certification point of view, DEQs have the definite advantage of having only one implicit layer compared to DNNs and therefore is potentially more amenable to sophisticated techniques (e.g. algebraic certificates of positivity) which rapidly face their limit even for classical DNNs with modest size (but long depth). Therefore monDEQs constitute a class of DNNs for which robustness certification, uncertainty propagation or Lipschicity could potentially be investigated in a more satisfactory way than classical networks. Contrary to DNNs, for which a variety of tools have been developped, certification of DEQ modeld is relatively open, the only available tool is the Lipschitz bound in [33].
Contribution
We present three general semialgebraic models of ReLU monDEQ for certification (p ∈ Z+∪{+∞}): • Robustness Model for network Lp robustness certification. • Lipschitz Model for network Lipschitz constant estimation with respect to any Lp norm. • Ellipsoid Model for ellipsoidal outer-approximation of the image by the network of a polyhedra or an ellipsoid.
All these models can be used for robustness certification, a common task which we consider experimentally. And the main originality of this work is to successfully apply the proposed techniques to monDEQ networks (especially for the Ellipsoid model which is a novel form of reachability analysis), which was not proposed before. Both Lipschitz and Ellipsoid models can in addition be used for further a posteriori analyses. Interestingly, all three models are given by solutions of semidefinite programs (SDP), obtained by Shor relaxation of a common semialgebraic representation of ReLU monDEQs. Our models are all evaluated to simple ReLU monDEQs on MNIST dataset similar as [45, 33] on the task of robustness certification. We demonstrate that all three models ourperform the approach of [33] and the Robustness Model being the most efficient. Our experiments also suggest that DEQs are much less robust to L∞ perturbations than L2 perturbations, in contrast with classical DNNs [36] which have been shown to be robust on MNIST dataset to the level we are considering (ε = 0.1).
Related works
Neural network certification is a challenging topic in machine learning, contributions include:
Robustness certification of DNNs Even with high test accuracy, DNNs are very sensitive to tiny input perturbations, see e.g. [42, 16]. Robustness to input perturbation has been investigated in many different works with various techniques, including SDP relaxation in SDP-cert [36], first-order dual SDP method [10], abstract interpretation with ERAN [15, 41], multi-neuron convex relaxations with PRIMA [31], GPU-based method GPUPoly [30] which can scale to large networks with one million neurons and 34 layers, scalable quantitative verification framework with PROVERO [4], LP relaxation in Reluplex [25], analytical certification with Fast-lin [44] and CROWN [47], and their extension to convolutional neural networks with CNN-Cert [7]. All these methods are restricted to DNNs or CNNs and do not directly apply to DEQs.
Lipschitz constant estimation of DNNs Lipschitz constant of DNNs can be used for robustness certification [42, 18]. Existing contributions include naive layer-wise product bound [22, 43], controlled ordinary differential equations (ODE) method [20], zonotope-based method [24] which can be further extended to generative models, Mixed-Integer Programming (MIP) [23] for exactly computing the local Lipschitz constant of ReLU networks, the LP relaxations for CNNs [48] and DNNs [28], as well as SDP relaxations [35, 14, 8].
Lipschitz constant estimation of DEQs The authors in [33] provide an optimization-free upper bound of the Lipschitz constant of monDEQs depending on network weights. This bound is valid for L2 norm and we present a general model for arbitrary Lp norm. Based on the works of [33], the authors in [38] introduced new parameterizations of DEQs that admit a Lipschitz bound during training via unconstrained optimization.

heading: 2 Preliminary Background and Notations
text: We consider the monotone operator equilibrium network (monDEQ) [45], with a single implicit hidden layer. The main difference between monDEQ and deep equilibrium network (DEQ) [2] is that strong monotonicity is enforced on weight matrix and activation functions to guarantee the convergence of fixed point iterations. The authors in [45] proposed various structures of implicit layer, we only consider fully-connected layers, investigation of more advanced convolutional layers is in our list of future works.
Network description: Denote by F : Rp0 → RK a fully-connected monDEQ for classification, where p0 is the input dimension and K is the number of labels. Let x0 ∈ Rp0 be the input variable and z ∈ Rp be the variable in the implicit layer. We consider the ReLU activation function, which is simply defined as ReLU(x) = max{0, x}, and the output of the monDEQs can be written as
F (x0) = Cz + c, z = ReLU(Wz + Ux0 + u), (monDEQ)
where W ∈ Rp×p,U ∈ Rp×p0 ,u ∈ Rp,C ∈ RK×p, c ∈ RK are parameters of the network. The vector-valued function F (x0) provides a score for each label i ∈ {1, . . . ,K} associated to the input x0, the prediction corresponds to the highest score, i.e., yx0 = arg maxi=1,...,K F (x0)i. As in [45], the matrix Ip −W is strongly monotone: there is a known m > 0 such that Ip −W mIp, this constraint can be enforced by specific parametrization of the matrix W. With the monotonicity assumption, the solution to equation z = ReLU(Wz + Ux0 + u) is unique and can be evaluated using convergent algorithms, see [45] for more details.
Robustness of monDEQs: Given an input x0 ∈ Rp0 , a norm ‖ · ‖, and a network F : Rp0 → RK , let y0 be the label of input x0, i.e., y0 = arg maxi=1,...,K F (x0)i. For ε > 0, denote by E = B(x0, ε, ‖ · ‖) the ball centered at x0 with radius ε for norm ‖ · ‖. If for all inputs x ∈ E , the label of x equals y0, i.e., y = arg maxi=1,...,K F (x)i = y0, then we say that the network F is ε-robust at input x0 for norm ‖ · ‖. An equivalent way to verify whether the network F is ε-robust is to check that for all labels i 6= y0, F (x)i − F (x)y0 < 0. Semialgebraicity of ReLU function: The key reason why neural networks with ReLU activation function can be tackled using polynomial optimization techniques is semialgebraicity of the ReLU function, i.e., it can be expressed with a system of polynomial (in)equalities. For x, y ∈ R, we have y = ReLU(x) = max{0, x} if and only if y(y−x) = 0, y ≥ x, y ≥ 0. For x,y ∈ Rn, we denote by ReLU(x) the coordinate-wise evaluation of ReLU function, and by xy the coordinate-wise product of x and y. A subset of Rn defined by a finite conjunction of polynomial (in)equalities is called a basic closed semialgebraic set. The graph of the ReLU function is a basic closed semialgebraic set.
Going back to equation (monDEQ), we have the following equivalence:
z = ReLU(Wz + Ux0 + u)⇔ z(z−Wz−Ux0 − u) = 0, z ≥Wz + Ux0 + u, z ≥ 0, (1) where the right hand side is a system of polynomial (in)equalities. For the rest of the paper, mention of the ReLU function will refer to the equivalent polynomial system in (1).
POP and Putinar’s positivity certificate: In general, a polynomial optimization problem (POP) has the form
ρ = max x∈Rn
{f(x) : gi(x) ≥ 0, i = 1, . . . , p} , (POP0)
where f and gi are polynomials whose degree is denoted by deg. The robustness certification model (3.1), Lipischitz constant model (3.2) and ellipsoid model (3.3) are all POPs.
In most cases, the POPs are non-linear and non-convex problems, which makes them NP-hard. A typical approach to reduce the complexity of these problems is replacing the positivity constraints by Putinar’s positivity certificate [34]. The problem (POP0) is equivalent to
ρ = min λ∈R {λ : λ− f(x) ≥ 0, gi(x) ≥ 0, i = 1, . . . , p, ∀x ∈ Rn} . (POP)
In order to reduce the size of the feasible set of problem (POP), we replace the positivity constraint λ− f(x) ≥ 0 by a weighted sum-of-square (SOS) polynomial decomposition, involving the polynomials gi. Let d be a non-negative integer. Denote by σd0(x), σ d i (x) some SOS polynomials of degree at most 2d, for each i = 1, . . . , p. Note that if d = 0, such polynomials are non-negative real numbers. Then the positivity of λ− f(x) is implied by the following decomposition
λ− f(x) = σd0(x) + p∑ i=1 σd−ωii (x)gi(x) , ωi = ddeg gi/2e , ∀x ∈ R n , (2)
for any d ≥ maxi ωi. Equation (2) is called the order-d Putinar’s certificate. By replacing the positivity constraint f(x) − λ ≥ 0 in problem (POP) by Putinar’s certificate (2), we have for d ≥ maxi ωi,
ρd = min λ∈R {λ : λ− f(x) = σd0(x) + p∑ i=1 σd−ωii (x)gi(x) , ωi = ddeg gi/2e,∀x ∈ R n} . (POP-d)
It is obvious that ρd ≥ ρ for all d ≥ maxi ωi. Under certain conditions (slightly stronger than compactness of the set of constraints), it is shown that limd→∞ ρd = ρ [27]. The main advantage of relaxing problem (POP) to (POP-d) is that problem (POP-d) can be efficiently solved by semidefinite programming (SDP). Indeed a polynomial f of degree at most 2d is SOS if and only if there exists a positive semidefinte (PSD) matrix M (called a Gram matrix) such that f(x) = v(x)TMv(x), for all x ∈ Rp, where v(x) is the vector of monomials of degree at most d. Problem (POP-d) is also called the order-d Lasserre’s relaxation. When the input polynomials are quadratic, the order-1 Lasserre’s relaxation is also known as Shor’s relaxation [40]. All our models are obtained using variations of Shor’s relaxation applied to different POPs, see Section 3.3 for more details.

heading: 3 Semialgebraic Models for Certifying Robustness of Neural Networks
text: In this section, we introduce several models for certification of monDEQs. All the models are based on semialgebraicity of ReLU and ∂ReLU (the subgradient of ReLU, see Section 3.2) to translate our targeted problems to POPs. Then we use Putinar’s certificates, defined in Section 2 to relax the non-convex problems to convex SDPs which can be solved efficiently using modern solvers.
Notation: Throughout this section, we consider a monDEQ for classification, denoted by F , with fixed, given parameters, W ∈ Rp×p,U ∈ Rp×p0 ,u ∈ Rp,C ∈ RK×p, c ∈ RK , where p0 is the number of input neurons, p is the number of hidden neurons, and K is the number of labels. For q ∈ Z+ ∪ {+∞}, ‖ · ‖q is the Lq norm defined by ‖x‖q := ( ∑p0 i=1 |xi|q)1/q for all x ∈ Rp0 . Throughout this section > 0 and x0 ∈ Rp0 are fixed, we denote by E := B(x0, ε, ‖ · ‖q) = {x ∈ Rp0 : ‖x − x0‖q ≤ ε} the ball centered at x0 with radius ε for Lq norm, a perturbation region. If q < +∞, i.e., q is a positive integer, ‖x − x0‖q ≤ ε is equivalent to the polynomial inequality ‖x − x0‖qq ≤ εq; if q = ∞, ‖x − x0‖q ≤ ε is equivalent to |x − x0|2 ≤ 2 (where |x| denotes the vector of absolute values of coordinates of x) which is a system of p0 polynomial inequalities. Hence the input set E is a semialgebraic set for all considered Lq norms. For a matrix A ∈ Rm×n, its operator norm induced by the norm ‖ · ‖ is given by |||A||| := inf{λ : ‖Ax‖ ≤ λ‖x‖,∀x ∈ Rn}.

heading: 3.1 Robustness Model
text: Let y0 be the label of x0 and let z ∈ Rp be the variables in the monDEQ implicit layer. The proposed model directly estimates upper bounds on the gap between the score of label y0 and the score of labels
different from y0. Precisely, for i ∈ {1, . . . ,K} such that i 6= y0, denote by ξi = (Ci,:−Cy0,:)T . For x ∈ E , the gap between its score of label i and label y0 is F (x)i − F (x)y0 = ξTi z. The Robustness Model for monDEQ reads:
δi := max x∈Rp0 ,z∈Rp
{ξTi z : z = ReLU(Wz + Ux + u), x ∈ E} . (CertMON-i)
Using the semialgebraicity of both ReLU in (1), and set E , problem (CertMON-i) is a POP for all i. As discussed in Section 2, one is able to derive a sequence of SDPs (Lasserre’s relaxation) to obtain a converging serie of upper bounds of the optimal solution of (CertMON-i). For Robustness Model, we consider only the order-1 Lasserre’s relaxation (Shor’s relaxation), and denote by δ̃i the upper bound of δi by Shor’s relaxation, i.e., δi ≤ δ̃i. Recall that if for all label i different from y0, we have F (x)i < F (x)y0 , then the label of x is still y0. This justifies the following claim:
Certification criterion: If δ̃i < 0 for all i 6= y0, then the network F is ε-robust at x0. Robustness Model for DNNs has already been investigated in [36], where the authors also use Shor’s relaxation as we do. Different from DNNs, we only have one implicit layer in monDEQ. Therefore, the number of variables in problem (CertMON-i) only depends on the number of input neurons p0 and hidden neurons p.

heading: 3.2 Lipschitz Model
text: We bound the Lipschitz constant of monDEQ with respect to input perturbation. Recall that the Lipschitz constant of the vector-valued function F (resp. z) w.r.t. the Lq norm and input ball S ⊃ E , denoted by LqF,S (resp. L q z,S), is the smallest value of L such that ‖F (x)− F (y)‖q ≤ L‖x− y‖q (resp. ‖z(x)− z(y)‖q ≤ L‖x− y‖q) for all x,y ∈ S . For x,x0 ∈ S , with ‖x− x0‖q ≤ , we can estimate the perturbation of the output as follows:
‖F (x)− F (x0)‖q ≤ LqF,S · ‖x− x0‖q ≤ εL q F,S , (3) ‖F (x)− F (x0)‖q ≤ |||C|||q · L q z,S · ‖x− x0‖q ≤ ε|||C|||q · L q z,S . (4)
The authors in [33] use inequality (4) with q = 2, as they provide an upper bound of L2z,S . In contrast, our model provides upper bounds on Lipschitz constants of F or z for arbitrary Lq norms. We directly focus on estimating the value of LqF,S instead of L q z,S .
Since the ReLU function is non-smooth, we define its subgradient, denoted by ∂ReLU, as the set-valued map ∂ReLU(x) = 0 for x < 0, ∂ReLU(x) = 1 for x > 0, and ∂ReLU(x) = [0, 1] for x = 0. Similar to ReLU function, ∂ReLU is also semialgebraic.
Semialgebraicity of ∂ReLU: If x, y ∈ R, we have y ∈ ∂ReLU(x), if and only if y(y−1) ≤ 0, xy ≥ 0, x(y − 1) ≥ 0. If x,y ∈ Rn, then ∂ReLU(x) denotes the coordinate-wise evaluation of ∂ReLU. Going back to monDEQ, let s be any subgradient of the implicit variable: s ∈ ∂ReLU(Wz+Ux+u). We can write equivalently a system of polynomial inequalities:
s(s− 1) ≤ 0, s(Wz + Ux + u) ≥ 0, (s− 1)(Wz + Ux + u) ≥ 0 . (5)
For the following discussion, ∂ReLU will refer to the equivalent polynomial systems (5).
With the semialgebraicity of ReLU in (1) and ∂ReLU in (5), one is able to compute upper bounds of the Lipschitz constant of F via POPs. The proof of the next Lemma is postponed to Appendix A.1.
Lemma 1 Define
L̃qF,S = max{t TUTy : t,x ∈ Rp0 , s, z,y, r ∈ Rp, v,w ∈ RK , x ∈ S,
‖t‖q ≤ 1, wTv ≤ 1, ‖w‖q ≤ 1, r−WTy = CTv, y = diag(s) · r; s ∈ ∂ReLU(Wz + Ux + u), z = ReLU(Wz + Ux + u)} . (LipMON)
Then L̃qF,S is an upper bound of the Lipschitz constant of F w.r.t. the Lq norm, i.e., L q F,S ≤ L̃ q F,S .
Since problem (LipMON) is a POP, we also consider Shor’s relaxation and denote by L̂qF,S the upper bound of L̃qF,S by Shor’s relaxation, i.e., L̃ q F,S ≤ L̂ q F,S . Define δ := εL̂ q F,S . By equations (3), (4), if
E ⊂ S , using Lemma 1 and the fact that ‖ · ‖∞ ≤ ‖ · ‖q , we have ‖F (x)− F (x0)‖∞ ≤ δ, yielding the following criterion:
Certification criterion: Let y0 be the label of x0. Define τ := F (x0)y0 − maxk 6=y0 F (x0)k. If 2δ < τ , then the network F is ε-robust at x0.
Remark: In order to avoid some possible numerical issues, we add some bound constraints and redundant constraints to problem (LipMON), see Appendix A.2 for details. Shor’s relaxation of Lipschitz Model for DNNs has already been extensively investigated in [14, 8]. If one want to certify robustness for several input test examples, then one may choose S to be a big ball containing all such examples with an additional margin of . Choosing a big ball for all input points requires to solve only one optimization problem, while choosing one ball for each input point (say among N ) requires to solve N optimization problems, which is much more costly. In this case it is more favorable to use the Certification Model directly.

heading: 3.3 Ellipsoid Model
text: As above, the input region E is a neighborhood of input x0 ∈ Rp0 with radius ε for the Lq norm, i.e., E = {x ∈ Rp : ‖x− x0‖q ≤ ε}. More generally, E could be other general semialgebraic sets, such as polytopes or zonotopes. Denote by F (E) the image of E by F . In this section, we aim at finding a semialgebraic set C such that F (E) ⊆ C. We choose C to be an ellipsoid which can in turn be used for robustness certification. Our goal is to find such outer-approximation ellipsoid with minimum volume. Let C := {ξ ∈ RK : ‖Qξ + b‖2 ≤ 1} be an ellipsoid in the output space RK parametrized by Q ∈ SK and b ∈ RK , where SK is the set of PSD matrices of size K × K. The problem of finding the minimum-volume ellipsoid containing the image F (E) can be formulated as
max Q∈SK ,b∈RK
{det(Q) : x ∈ E , ‖Q(Cz + c) + b‖2 ≤ 1, z = ReLU(Wz + Ux + u)} .
(EllipMON-POP)
By semialgebraicity of both ReLU in (1) and C, problem (EllipMON-POP) can be cast as a POP (the determinant being a polynomial). We no longer apply Shor’s relaxation, we rather replace the non-negativity output constraint 1−‖Q(Cz+c)+b‖2 ≥ 0 by a stronger Putinar’s certificate related to both ReLU and input constraints. We can then relax the non-convex problem (EllipMON-POP) to a problem with SOS constraints, which can be reformulated by (convex) SDP constraints. This is due to the fact that a polynomial f of degree at most 2d is SOS if and only if there exists a PSD matrix M (called a Gram matrix) such that f(x) = v(x)TMv(x), for all x ∈ Rp, with v(x) being the vector containing all monomials of degree at most d. In summary, we relax problem (EllipMON-POP) to an SOS constrained problem keeping the determinant unchanged:
max Q∈SK ,b∈RK
{det(Q) : 1− ‖Q(Cz + c) + b‖22 = σ0(x, z) + σ1(x, z)T gq(x− x0)
+ τ(x, z)T (z(z−Wz−Ux− u)) + σ2(x, z)T (z−Wz−Ux− u) + σ3(x, z)T z} . (EllipMON-SOS-d)
where gq(x) = εq − ‖x‖qq for q < +∞ and gq(x) = ε2 − |x|2 for q = +∞, σ0 is a vector of SOS polynomials of degree at most 2d, σ1 is a vector of SOS polynomials of degree at most 2(d− dq/2e) for q < +∞ and 2d − 2 for q = +∞, σ2, σ3 are vectors of SOS polynomials of degrees at most 2d− 2, τ is a vector of polynomials of degree at most 2d− 2. Problem (EllipMON-SOS-d) provides an ellipsoid feasible for (EllipMON-POP), that is an ellipsoid which contains F (E). In practice, the determinant is replaced by a log-det objective because there exist efficient solver dedicated to optimize such objectives on SDP constraints. By increasing the relaxation order d in problem (EllipMON-SOS-d), one is able to obtain a hierarchy of log-det objected SDP problems for which the outer-approximation ellipsoids have decreasing volumes.
In this paper, we only consider the case p = 2,∞, and order-1 relaxation (d = 1). Therefore, σi and τ are all (vectors of) real (non-negative) numbers for i = 1, 2, 3, except that σ0 is an SOS polynomial of degree at most 2. In this case, problem (EllipMON-SOS-d) is equivalent to a problem with log-det objective and SDP constraints, as the following lemma states (proof postponed to Appendix A.3):
Lemma 2 For p = 2 or p =∞, problem (EllipMON-SOS-d) with d = 1 is equivalent to max
Q∈SK ,b∈RK ,σ1,σ2,σ3≥0,τ∈Rp {log det(Q) : −M 0} . (EllipMON-SDP)
where M ∈ S(p0+p+1)×(p0+p+1) is a symmetric matrix parametrized by the decision variables (Q,b), the coefficients (σ1, σ2, σ3, τ), and the parameters of the network (W,U,u,C, c).
Since the outer-approximation ellipsoid C = {ξ ∈ RK : ‖Qξ + b‖2 ≤ 1} contains the image F (E), i.e., all possible outputs of the input region E , one is able to certify robustness by solving the following optimization problems.
Certification criterion: Let y0 be the label of x0. For i 6= y0, define δi := maxξ∈RK{ξi − ξy0 : ‖Qξ + b‖2 ≤ 1}. If δi < 0 for all i 6= y0, then the network F is ε-robust at x0. The certification criterion for Ellipsoid Model has a geometric explanation: for i 6= y0, denote by Pi the projection map from output space RK to its 2-dimensional subspace Ry0 × Ri, i.e., Pi(ξ) = [ξy0 , ξi]T for all ξ ∈ RK . Let Li be the line in subspace Ry0 × Ri defined by {[ξy0 , ξi]T ∈ Ry0 × Ri : ξy0 = ξi}. Then the network F is ε-robust if the projection Pi(C) lies strictly below the line Li for all i 6= y0. We give an explicit example in Section 4.3 to visually illustrate this.

heading: 3.4 Summary of the Models
text: We have already presented three models which can all be dedicated to certify robustness of neural networks. However, the size and complexity of each model are different. We summarize the number of variables in each model and the maximum size of PSD matrices in the resulting Shor’s relaxation, see Table 1. The complexity of our models only depends on the number of neurons in the input layer and implicit layer. The size of PSD matrices is a limiting factor for SDP solvers, our models are practically restricted to network for which such size can be handled by SDP solvers. For the popular dataset MNIST [46], whose input dimension is 28× 28 = 784, we are able to apply our model on monDEQs with moderate size implicit layers (87) and report the corresponding computation time.

heading: 4 Experiments
text: In this section, we present the experimental results of Robustness Model, Lipschitz Model and Ellipsoid Model described in Section 3 for a pretrained monDEQ on MNIST dataset. The network we use consists of a fully-connected implicit layer with 87 neurons and we set its monotonicity parameter m to be 20. The training hyperparameters are set to be the same as in Table D1 of [45], where the training code (in Python) is available at https://github.com/locuslab/monotone_op_net. Training is based on the normalized MNIST database in [45], we use the same normalization setting on each test example with mean µ = 0.1307 and standard deviation σ = 0.3081, which means that each input is an image of size 28× 28 with entries varying from −0.42 to 2.82. And for every perturbation ε, we also take the normalization into account, i.e., we use the normalized perturbation ε/σ for each input.
Since all our three models can be applied to certify robustness of neural networks, we first compare the performance of each model in certification of the first 100 test MNIST examples. Then we compare the upper bounds of Lipschitz Model with the upper bounds proposed in [33]. Finally we show that Ellipsoid Model can also be applied for reachability analysis. For Certification model and Lipschitz model, we implement them in Julia [5] with JuMP [11] package; for Ellipsoid model, we implement it in Matlab [39] with CVX [17] package. For all the three models, we use Mosek [29] as a backend to solve the targeted POPs. All experiments are performed on a personal laptop with an Intel 8-Core i7-8665U CPU @ 1.90GHz Ubuntu 18.04.5 LTS, 32GB RAM. The code of all our models is available at https://github.com/NeurIPS2021Paper4075/SemiMonDEQ.

heading: 4.1 Robustness certification
text: We consider ε = 0.1 for the L2 norm and ε = 0.1, 0.05, 0.01 for the L∞ norm. For each model, we compute the ratio of certified test examples among the first 100 test inputs. Following [33], we also compute the projected gradient descent (PGD) attack accuracy using Foolbox library [37], which indicates the ratio of non-successful attacks among our 100 inputs. Note that the ratio of certified examples should always be less or equal than the ratio of non-successful attacks. The gaps between them shows how many test examples there are for which we are neither able to certify robustness nor find adversarial attacks.
Remark: For Lipschitz Model, we use inequality (3) to test robustness, i.e., we compute directly the upper bound of the Lipschitz constant of F rather than z (seen as a function of x) where S is a big ball containing all test examples.
From Table 2, we see that the monDEQ is robust to all the 100 test examples for the L2 norm and ε = 0.1 (the only example that we can not certify is because the label itself is wrong). However, it is not robust for the L∞ norm at the same level of perturbation (all our three models cannot certify any examples), and the PGD algorithm finds adversarial examples for 85% of the inputs. The network becomes robust again for the L∞ norm when we reduce the perturbation ε to 0.01. Overall, we see that Robustness Model is the best model as it provides the highest ratio, Ellipsoid Model is the second best model compared to Robustness Model, and Lipschitz Model provides the lowest ratio. As a trade-off, for each test example, Robustness Model requires to consider at most 9 optimization problems, each one being solved in around 150 seconds, while Ellipsoid Model requires to consider only one problem, which is solved in around 500 seconds. We only need to calculate one (global) Lipschitz constant, which takes around 1500 seconds, so that we are able to certify any number of inputs. Each model we propose provide better or equal certification accuracy compared to [33], and significant improvements for L∞ perturbations.

heading: 4.2 Comparison with Lipschitz constants
text: In this section, we compare the upper bounds of Lipschitz constants computed by Lipschitz Model with the upper bounds proposed in [33]. Notice that the upper bounds in [33] only involve the function z(x), hence we are only able to use inequality (4) to test robustness. In fact the quantity |||C|||q ·L q z,S can be regarded as an upper bound of LqF,S , the Lipschitz constant of F . We denote by UB 2 z the upper bound of the Lipschitz constant of z w.r.t. the L2 norm, given by UB2z = |||U|||2/m according to [33], where U is the parameter of the network and m is the monotonicity factor. We can then compute the upper bound w.r.t. the L∞ norm by UB∞z = √ p0 · UB2z where p0 is the input dimension. The upper bound of Lipschitz constant of F is computed via the upper bound of z: UBqF = |||C|||q · UB q z. Denote similarly by SemiUBqF the upper bounds of Lipschitz constants of F provided by Lipschitz Model, w.r.t. the Lq norm.
From Table 3, we see that Lipschitz Model provides consistently tighter upper bounds than the ones in [33]. Especially for L∞ norm, the upper bound computed by |||C|||∞ ·UB ∞ z is rather crude compared to the bound obtained directly by Lipschitz Model. Therefore, we are able to certify more examples using SemiUBqF than UB q F , see Table 2.

heading: 4.3 Outer ellipsoid approximation
text: In this section, we provide a visible illustration of how Ellipsoid Model can be applied to certify robustness of neural networks.
Take the first MNIST test example (which is classified as 7) for illustration. For ε = 0.1, this example is certified to be robust for the L2 norm but not for the L∞ norm. We show the landscape of the projections onto R7 × R3, i.e., the x-axis indicates label 7 and the y-axis indicates label 3. In Figure 1, the red points are projections of points in the image F (E), for E an L2 or L∞ norm perturbation zone, the black circles are projections of some (successful and unsuccessful) adversarial examples found by the PGD algorithm. Notice that the adversarial examples also lie in the image F (E). The blue curve is the boundary of the projection of the outer-approximation ellipsoid (which is an ellipse), and the blue dashed line plays the role of a certification threshold. Figure 1a shows the landscape for the L2 norm, we see that the ellipse lies strictly below the threshold line, which means that for all points ξ ∈ C, we have ξ3 < ξ7. Hence for all ξ ∈ F (E), we also have ξ3 < ξ7. On the other hand, for the L∞ norm, we see from Figure 1b that the threshold line crosses the ellipse, which means that we are not able to certify robustness of this example by Ellipsoid Model. Indeed, we can find adversarial examples with the PGD algorithm, as shown in Figure 1b by the black circles that lie above the threshold line. The visualization of one of the attack examples is shown in Figure 2 in Appendix A.4.

heading: 5 Conclusion and Future Works
text: In this paper, we introduce semialgebraic representations of monDEQ and propose several POP models that are useful for certifying robustness, estimating Lipschitz constants and computing outerapproximation ellipsoids. For each model, there are several hierarchies of relaxations that allow us to improve the results by increasing the relaxation order. Even though we simply consider the order-1 relaxation, we obtain tighter upper bounds of Lipschitz constants compared to the results in [33]. Consequently, we are able to certify robustness of more examples.
Our models are based on SDP relaxation, hence requires an efficient SDP solver. However, the statof-the-art SDP solver Mosek (by interior-point method) can only handle PSD matrices of moderate size (smaller than 5000). This is the main limitation of our method if the dimension of the input gets larger. Moreover, we only consider the fully-connected monDEQ based on MNIST datasets for illustration. One important and interesting future work is to generalize our model to single and multi convolutional monDEQ, and to other datasets such as CIFAR [26] and SVHN [32]. Directly using off-the-shelves interior point SDP solvers to solve the problems for convolutional networks is not possible because of their current size limitation. Fortunately, a convolutional layer can be regarded as a fully-connected layer with a larger (but sparse) weight matrix. Hence one is able to build similar models via sparse polynomial optimization tools.
The authors in [10] provide an interesting and promising first-order method as an alternative to the costly interior point methods for solving semidefinite relaxations associated with robustness certification of DNNs. What is crucial in [10] is to exploit the network layer structure for efficient back-propagation in gradient computation. In particular, the algorithm requires memory only linear in the total number of network activations and only requires a fixed number of forward/backward passes through the network per iteration. This enables the algorithm to certify robustness of largescaled networks more efficiently and accurately. However, the auto-differentiation technique for computing subgradients in DNNs does not apply directly to monDEQs, since computing subgradients in monDEQs involves solving fixed-point equations. Therefore, adapting the approach of [10] for solving more efficiently our semidefinite relaxations associated with monDEQs is certainly worth considering but not straightforward. It is a topic of further investigation.
The authors in [33] showed that we can train DEQs with small Lipschitz constants for the L2 norm, by controlling the monotonicity of the weight matrix. This guarantees the robustness of monDEQ w.r.t. the L2 norm but not for the L∞ norm. A natural investigation track is to adapt this training technique to the L∞ norm with a better control of the associated Lipschitz constant.
Acknowledgments and Disclosure of Funding
The authors acknowledge the support of AI Interdisciplinary Institute ANITI funding, through the French “Investing for the Future – PIA3” program under the Grant agreement ANR-19-PI3A-0004. Edouard Pauwels acknowledges the support of Air Force Office of Scientific Research, Air Force Material Command, USAF, under grant numbersFA9550-19-1-702 6and ANR MaSDOL 19-CE230017-01. This work was supported by the National Research Foundation, Prime Minister’s Office, Singapore under its Campus for Research Excellence and Technological Enterprise (CREATE) programme. Victor Magron was supported by the FMJH Program PGMO (EPICS project) and EDF, Thales, Orange et Criteo, as well as from the Tremplin ERC Stg Grant ANR-18-ERC2-0004-01 (T-COPS project). This work has benefited from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie Actions.

Review summaries:
This paper introduces a semialgebraic representation of the monotone deep equilibrium model (monDEQ), a type of implicit-depth equilibrium network which guarantees the convergence of fixed-point iterations. This representation is used to derive SDP models for 1) certifying robustness, 2) estimating the Lipschitz constant, and 3) ellipsoid uncertainty propagation. These techniques show improved performance over existing methods for certifying robustness and bounding the Lipschitz constant of monDEQs. In addition, they apply to general $\ell_p$ norms, unlike existing work. 

The authors propose a framework for the robustness verification of monDEQs (monotone operator deep equilibrium networks). In this framework the paper introduces three different methods: one for direct robustness verification, one for Lipschitz constant estimation and one to compute the ellipsoid over-approximating all possible outputs for a set of inputs (in a norm ball). The authors also show how the Lipschitz constant and ellipsoid in turn can be used for robustness verification.  Mathematically this is enabled encoding the monDEQ as the constraints of a (polynomial) optimization problem. In particular, the fix-point $z = ReLU(Wz + Ux_{0} + u)$ is encoded using the semialgebracity of ReLU, similar to [27]*. This encoding as well as encoding $\ell_{2}$ and $\ell_{\infty}$ inputs yields quadratic constraints. (Other $\ell_{p}$-norms $p>2$ yield higher-order constraints.) The resulting constraints are relaxed via Puntinar's positivity certificate and Lasserre's/Shor's relaxations into a SDP (semidefinite program), e.g. approximated lower-degree polynomials.  For the robustness verification approach the authors follow [27] and maximize the difference between any logit and the target logit. If all of these differences are $< 0$ a robustness certificate is obtained. To make the resulting problem tractable Shor's relaxation is applied. The Lipschitz approach bounds the Lipschitz constant by bounding the supremum of the operator norm of the generalized Jacobian. This in turn, via Shor's relaxation gives an SDP. The Ellipsoid model minimizes the volume of an ellipsoid that covers the output. Via a custom relaxation this also yields an SDP.  All methods are evaluated on an MNIST classifier with a single feed-forward fix-point layer where $z \in \mathbb{R}^{87}$. All three approaches allow for robustness certification (although the robustness formulation performs the best). The computed Lipschitz bound are close to those obtained by [24] for the $\ell_{2}$-norm but much tighter for the $\ell_{\infty}$-norm. Lastly, in an example the authors showcase how the output-ellipsoid method allows insights into the model robustness. 

Based on semialgebraic of $\text{ReLU}$ and $\partial\text{ReLU}$, this paper proposes three models that are useful for certifying robustness, estimating Lipschitz constants and computing outer approximation ellipsoids.  The authors also use Putinar’s certificates to relax the non-convex problems to convex SDPs which can be solved efficiently. Experiments on MNIST demonstrate the effectiveness of all three models.