Implicit Deep Learning

TL;DR

Implicit deep learning introduces a framework where neural network predictions are defined by fixed-point equations, simplifying notation and enabling new research directions in architecture, robustness, interpretability, and optimization.

Contribution

The paper formalizes implicit deep learning prediction rules, expanding the theoretical foundation and opening avenues for novel architectures and analysis methods.

Findings

Simplifies deep learning notation

Enables new architectures and algorithms

Facilitates robustness and interpretability analysis

Abstract

Implicit deep learning prediction rules generalize the recursive rules of feedforward neural networks. Such rules are based on the solution of a fixed-point equation involving a single vector of hidden features, which is thus only implicitly defined. The implicit framework greatly simplifies the notation of deep learning, and opens up many new possibilities, in terms of novel architectures and algorithms, robustness analysis and design, interpretability, sparsity, and network architecture optimization.

Tables

Table 1. Table 1: Experimental results of attack success rate against percentage of perturbed inputs on MNIST and CIFAR-10 ( 10 000 samples from test set).

	Sensitivity matrix attack		Gradient-based attack
% of perturbed inputs	MNIST	CIFAR-10	MNIST	CIFAR-10
0.1%	1.01%	3.04%	2.42%	1.75%
1%	13.41%	10.16%	26.92%	6.66%
10%	70.67%	36.21%	74.90%	33.18%
20%	89.82%	57.01%	87.10%	52.57%
30%	90.22%	67.45%	89.82%	66.59%

Equations

$$\hat{y}(u)$$

$$x$$

$$\Vert A{\Vert }_{p\to q}=\underset{\xi }{\max }\Vert A\xi {\Vert }_q:\Vert \xi {\Vert }_p\le 1.$$

$$\Vert A{\Vert }_{\infty }:=\underset{i}{\max }\sum _j|A_{ij}|.$$

$$\eta (z):=(\Vert z_1{\Vert }_{p_1},\dots ,\Vert z_L{\Vert }_{p_L}{)}^{\top }.$$

$$\forall u,v\in {\mathbb{R}}^{n_l}:\Vert {\varphi }_l(u)-{\varphi }_l(v){\Vert }_{p_l}\le {\gamma }_l\Vert u-v{\Vert }_{p_l}.$$

$$\forall u,v\in {\mathbb{R}}^n:|\varphi (u)-\varphi (v)|\le |u-v|,$$

$$z\to \text{mbox}SoftMax(z):={\left(\frac{e^{z_i}}{{\sum }_i{e}^{z_j}}\right)}_{i\in [n]},$$

$$x=\varphi (Ax+b)$$

$$x(0)=0,x(t+1)=\varphi (Ax(t)+b),t=0,1,2,\dots .$$

$$x_n=\frac{(b_n{)}_{+}}{1-A_{nn}},x_i=\frac{1}{1-A_{ii}}(b_i+\sum _{j>i}{A}_{ij}{x}_j{)}_{+},i=n-1,\dots ,1.$$

$$(N(A){)}_{ij}:=\Vert A{\Vert }_{p_j\to p_i}=\underset{\xi }{\max }\Vert A\xi {\Vert }_{p_i}:\Vert \xi {\Vert }_{p_j}\le 1.$$

$${\lambda }_{\mathrm{pf}}(\Gamma N(A))<1,$$

$$A:=\left(\begin{array}{cc}{A}_{11}& A_{12}\\ 0& A_{22}\end{array}\right)$$

$$A:=\left(\begin{array}{cc}{A}_{11}& 0\\ A_{21}& A_{22}\end{array}\right),$$

$$A:=\left(\begin{array}{cc}{A}_{11}& A_{12}\\ 0& A_{22}\end{array}\right)$$

$${\hat{y}}_2=C_2{x}_2+D_2{u}_2\text{mbox}where{u}_2={\hat{y}}_1=C_1{x}_1+D_1{u}_1,\text{mbox}where{x}_i={\varphi }_i(A_i{x}_i+B_i{u}_i),i=1,2.$$

$$\left(\begin{array}[]{c|c}A&B\\ \hline\cr C&D\end{array}\right)=\left(\begin{array}[]{cc|c}A_{2}&B_{2}C_{1}&B_{2}D_{1}\\ 0&A_{1}&B_{1}\\ \hline\cr C_{2}&D_{2}C_{1}&D_{2}D_{1}\end{array}\right).$$

$$(u_1,u_2)\in {\mathbb{R}}^2\to \varphi (u)={\varphi }_1(u_1){\varphi }_2(u_2)$$

$$y_i=C_i{x}_i+D_i{u}_i,x_i={\varphi }_i(A_i{x}_i+B{u}_i),i=1,2.$$

$$\left(\begin{array}[]{c|c}A&B\\ \hline\cr C&D\end{array}\right)=\left(\begin{array}[]{cc|c}A_{1}&B_{1}C_{2}&B_{1}\\ B_{2}C_{1}&A_{2}&0\\ \hline\cr C_{1}&0&0\end{array}\right).$$

$${\lambda }_{\mathrm{pf}}(|A|)=\underset{S}{\inf }\Vert SA{S}^{-1}{\Vert }_{\infty }:S=\mathrm{diag}(s),s>0.$$

$$\left(\begin{array}[]{c|c}A^{\prime}&B^{\prime}\\ \hline\cr C^{\prime}&D^{\prime}\end{array}\right)=\left(\begin{array}[]{c|c}SAS^{-1}&SB\\ \hline\cr CS^{-1}&D^{\prime}\end{array}\right),$$

$$\hat{y}(u)={\varphi }_o(Cx+Du),x=\varphi (Ax+Bu),$$

$$\left(\begin{array}{c}z\\ x\end{array}\right)=\tilde{\varphi }\left(\left(\begin{array}{cc}0& C\\ 0& A\end{array}\right)\left(\begin{array}{c}z\\ x\end{array}\right)+\left(\begin{array}{c}D\\ B\end{array}\right)u\right),y=z.$$

$$y=Cx+Du+d,x=\varphi (Ax+Bu+b),$$

$$y=Cx+Du,x=\varphi (Ax+Bu).$$

$$y=C\tilde{x}+\left(\begin{array}{cc}D& -C\varphi (0)\end{array}\right)\left(\begin{array}{c}u\\ 1\end{array}\right),\tilde{x}=\tilde{\varphi }(A\tilde{x}+\left(\begin{array}{cc}B& -A\varphi (0)\end{array}\right)\left(\begin{array}{c}u\\ 1\end{array}\right)).$$

$$y=Du+d,\text{mbox}where[D,d]:=\mathrm{diag}(\sigma {)}^{-1}[I,-\hat{u}].$$

$$\hat{y}(u)=W_L{x}_L,x_{l+1}={\varphi }_l(W_l{x}_l),x_0=u.$$

$$\left(\begin{array}[]{c|c}A&B\\ \hline\cr C&D\end{array}\right)=\left(\begin{array}[]{cccc|c}0&W_{L-1}&\ldots&0&0\\ &0&\ddots&\vdots&\vdots\\ &&\ddots&W_{1}&0\\ &&&0&W_{0}\\ \hline\cr W_{L}&0&\ldots&0&0\end{array}\right),$$

Full text

{AMS}

690C26, 49M99, 65K10, 62M45, 26B10

1 Introduction

1.1 Implicit prediction rules

In this paper, we consider a new class of deep learning models that are based on implicit prediction rules. Such rules are not obtained via a recursive procedure through several layers, as in current neural networks. Instead, they are based on solving a fixed-point equation in some single “state” vector $x\in{\mathbb{R}}^{n}$. Precisely, for a given input vector $u$, the predicted vector is

[TABLE]

where $\phi\>:\>{\mathbb{R}}^{n}\rightarrow{\mathbb{R}}^{n}$ is a nonlinear vector map (the “activation” map), and matrices $A,B,C,D$ contain model parameters. Fig. 1 provides a block-diagram view of an implicit model, to be read from right to left, so as to be consistent with matrix-vector multiplication rules.

We can think of the vector $x\in{\mathbb{R}}^{n}$ as a “state” corresponding to $n$ “hidden” features that are extracted from the inputs, based on the so-called equilibrium equation Eq. 1b. In general, that equation cannot be solved in closed-form, and the model above provides $x$ only implicitly. This equation is not necessarily well-posed, in the sense that it may not admit a solution, let alone a unique one; we discuss this important issue of well-posedness in Section 2.

For notational simplicity only, our rule does not contain any bias terms; we can easily account for those by considering the vector $(u,1)$ instead of $u$, thereby increasing the column dimension of $B$ by one.

Perhaps surprisingly, as seen in Section 3, the implicit framework includes most current neural network architectures as special cases. Implicit models are a much wider class: they present much more capacity, as measured by the number of parameters for a given dimension of the hidden features; also, they allow for cycles in the network, which is not permitted under the current paradigm of deep networks.

Implicit rules open up the possibility of using novel architectures and prediction rules for deep learning, which are not based on any notion of “network” or “layers”, as is classically understood. In addition, they allow one to consider rigorous approaches to challenging problems in deep learning, ranging from robustness analysis, sparsity and interpretability, and feature selection.

1.2 Contributions and paper outline

Our contributions in this paper, and its outline, are as follows.

• •

Well-posedness and composition (Section 2): In contrast with standard deep networks, implicit models may not be well-posed, in the sense that the equilibrium equation may have no or multiple solutions. We establish rigorous and numerically tractable conditions for implicit rules to be well-posed. These conditions are then used in the training problem, guaranteeing the well-posedness of the learned prediction rule. We also discuss the composition of implicit models, via cascade connections for example.

• •

Implicit models of neural networks (Section 3): We provide details on how to represent a wide variety of neural networks as implicit models, building on the composition rules of Section 2.

• •

Robustness analysis (Section 4): We describe how to analyze the robustness properties of a given implicit model, deriving bounds on the state under input perturbations, and generating adversarial attacks. We also discuss which penalties to include into the training problem so as to encourage robustness of the learned rule.

• •

Interpretability, sparsity, compression and deep feature selection (Section 5): Here we focus on finding appropriate penalties to use in order to improve properties such as model sparsity, or obtain feature selection. We also discuss the impact of model errors.

• •

Training problem: formulations and algorithms (Section 6): Informed by our previous findings, we finally discuss the corresponding training problem. Following the work of [24] and [34], we represent activation functions using so-called Fenchel divergences, in order to relax the training problem into a more tractable form. We discuss several algorithms, including stochastic projected gradients, Frank-Wolfe, and block-coordinate descent.

Finally, Section 7 provides a few experiments supporting the theory put forth in this paper. Our final Section 8 is devoted to prior work and references.

1.3 Notation

For a matrix $U$, $|U|$ (resp. $U_{+}$) denotes the matrix with the absolute values (resp. positive part) of the entries of $U$. For a vector $v$, we denote by $\mathop{\bf diag}(v)$ the diagonal matrix formed with the entries of $v$; for a square matrix $V$, $\mathop{\bf diag}(V)$ is the vector formed with the diagonal elements of $V$. The notation $\mathbf{1}$ refers to the vector of ones, with size inferred from context. The Hadamard (componentwise) product between two $n$-vectors $x,y$ is denoted $x\odot y$. We use $s_{k}(z)$ to denote the sum of the largest $k$ entries of a vector $z$. For a matrix $A$, and integers $p,q\geq 1$, we define the induced norm

[TABLE]

The case when $p=q=\infty$ corresponds to the $l_{\infty}$-induced norm of $A$, also known as its max-row-sum norm:

[TABLE]

We denote the set $\{1,\cdots,L\}$ compactly as $[L]$. For a $n$-vector partitioned into $L$ blocks, $z=(z_{1},\ldots,z_{L})$, with $z_{l}\in{\mathbb{R}}^{n_{l}}$, $l\in[L]$, with $n_{1}+\ldots+n_{L}=n$, we denote by $\eta(z)$ the $L$-vector of norms:

[TABLE]

Finally, any square, non-negative matrix $M$ admits a real eigenvalue that is larger than the modulus of any other eigenvalue; this non-negative eigenvalue is the so-called Perron-Frobenius eigenvalue [39], and is denoted $\lambda_{\rm pf}(M)$.

2 Well-Posedness and Composition

2.1 Assumptions on the activation map

We restrict our attention to activation maps $\phi$ that obey a “Blockwise LIPschitz” (BLIP) continuity condition. This condition is satisfied for most popular activation maps, and arises naturally when “composing” implicit models (see Section 2.4). Precisely, we assume that:

Blockwise: the map $\phi$ acts in a block-wise fashion, that is, there exist a partition of $n$: $n=n_{1}+\ldots+n_{L}$ such that for every vector partitioned into the corresponding blocks: $z=(z_{1},\ldots,z_{L})$ with $z_{l}\in{\mathbb{R}}^{n_{l}}$, $l\in[L]$, we have $\phi(z)=(\phi_{1}(z_{1}),\ldots,\phi_{L}(z_{L}))$ for appropriate maps $\phi_{l}\>:\>{\mathbb{R}}^{n_{l}}\rightarrow{\mathbb{R}}^{n_{l}}$, $l\in[L]$. 2. 2.

Lipschitz: For every $l\in[L]$, the maps $\phi_{l}$ are Lipschitz-continuous with constant $\gamma_{l}>0$ with respect to the $l_{p_{l}}$-norm for some integer $p_{l}\geq 1$:

[TABLE]

In the remainder of the paper, we refer to such maps with the acronym BLIP, omitting the dependence on the underlying structure information (integers $n_{l}$, $p_{l}$, $\gamma_{l}$, $l\in[L]$). We shall consider a special case, referred to a COmponentwise Non-Expansive (CONE) maps, when $n_{l}=1$, $\gamma_{l}=1$, $l\in[L]$. Such CONE maps satisfy

[TABLE]

with inequality and absolute value taken componentwise. Examples of CONE maps include the ReLU (defined as $\phi(\cdot)=\max(0,\cdot)$) and its “leaky” variants, tanh, sigmoid, each applied componentwise to a vector input. Our model also allows for maps that do not operate componentwise, such as the softmax function, which operates on a $n$-vector $z$ as:

[TABLE]

The softmax map is 1-Lipschitz-continuous with respect to the $l_{1}$-norm [21].

2.2 Well-posed matrices

We consider the prediction rule Eq. 1a with input point $u\in{\mathbb{R}}^{p}$ and predicted output vector $\hat{y}(u)\in{\mathbb{R}}^{q}$. The equilibrium equation Eq. 1b does not necessarily have a well-defined, unique solution $x$, as Fig. 2 illustrates in a scalar case. In order to prevent this, we assume that the $n\times n$ matrix $A$ satisfies the following well-posedness property.

Definition 2.1** (Well-posedness property).**

The $n\times n$ matrix $A$ is said to be well-posed for $\phi$ (in short, $A\in\mbox{\rm WP}(\phi)$) if, for any $n$-vector $b$, the equation in $x\in{\mathbb{R}}^{n}$:

[TABLE]

*has a unique solution. *

There are many classes of matrices that satisfy the well-posedness property. As seen next, strictly upper-triangular matrices are well-posed with respect to any activation map that acts componentwise; such a class arises when modeling feedforward neural networks as implicit models, as seen in Section 3.2.

2.3 Tractable sufficient conditions for well-posedness

Our goal now is to understand how we can constrain $A$ to have the well-posedness property, in a numerically tractable way.

A sufficient condition is based on the contraction mapping theorem. The following result focuses on the case when $\phi$ is componentwise non-expansive (CONE), as defined in Section 2.1.

Theorem 2.2** (PF sufficient condition for well-posedness with CONE activation).**

Assume that $\phi$ is a componentwise non-expansive (CONE) map, as defined in Section 2.1. Then, $A$ is well-posed with respect to $\phi$ if $\lambda_{\rm pf}(|A|)<1$, in which case, for any $n$-vector $b$, the solution to the equation Eq. 5 can be computed via the fixed-point iteration

[TABLE]

The theorem is a direct consequence of the global contraction mapping theorem [44, p.83]. We provide a proof in Appendix A. A few remarks are in order.

Remark 2.3**.**

*The fixed-point iteration Eq. 6 has linear convergence; each iteration is a matrix-vector product, hence the complexity is comparable to that of a forward pass through a network of similar size. *

Remark 2.4**.**

*The PF condition $\lambda_{\rm pf}(|A|)<1$ is not convex in $A$, but the convex condition $\|A\|_{\infty}<1$, is sufficient, in light of the bound $\|A\|_{\infty}\geq\lambda_{\rm pf}(|A|)$. *

Remark 2.5**.**

The PF condition of Theorem 2.2 is conservative. For example, a triangular matrix $A$ is well-posed with respect to the ReLU and if only if $\mathop{\bf diag}(A)<\mathbf{1}$, a consequence of the upcoming Theorem 2.9. The corresponding equilibrium equation can then be solved via the backward recursion

[TABLE]

*Such a matrix does not necessarily satisfy the PF condition; we can have in particular $A_{11}<-1$, which implies $\lambda_{\rm pf}(|A|)>1$. *

Remark 2.6**.**

*The well-posedness property is invariant under row and column permutation, provided $\phi$ acts componentwise. Precisely, if $A$ is well-posed with respect to a componentwise CONE map $\phi$, then for any $n\times n$ permutation matrix $P$, $PAP^{\top}$ is well-posed with respect to $\phi$. The PF sufficient condition is also invariant under row and column permutations. *

In some contexts, the map $\phi$ does not satisfy the componentwise non-expansiveness condition, but the weaker blockwise Lipschitz continuity (BLIP) defined in Section 2.1. The previous theorem can be extended to this case, as follows. We partition the $A$ matrix according to the tuple $(n_{1},\ldots,n_{L})$, into blocks $A_{ij}\in{\mathbb{R}}^{n_{i}\times n_{j}}$, $1\leq i,j\leq L$, and define a $L\times L$ matrix of induced norms, with elements for $l,h\in[L]$ given by

[TABLE]

Theorem 2.7** (PF sufficient condition for well-posedness for BLIP activation).**

Assume that $\phi$ satisfies the BLIP condition, as defined in Section 2.1. Then, $A$ is well-posed with respect to $\phi$ if

[TABLE]

*where $\Gamma:=\mathop{\bf diag}(\gamma)$, with $\gamma$ the vector of Lipschitz constants, and $N(\cdot)$ is the matrix of induced norms defined in Eq. 7. In this case, for any $n$-vector $b$, the solution to the equation Eq. 5 can be computed via the fixed-point iteration Eq. 6. *

Proof 2.8**.**

*See Appendix B. *

2.4 Composition of implicit models

Implicit models can be easily composed via matrix algebra. Sometimes, the connection preserves well-posedness, thanks to the following result.

Theorem 2.9** (Well-posedness of block-triangular matrices, componentwise activation).**

Assume that the activation map $\phi$ acts componentwise. The upper block-triangular matrix

[TABLE]

*with $A_{ii}\in{\mathbb{R}}^{n_{i}\times n_{i}}$, $i=1,2$, is well-posed with respect to $\phi$ if and only if its the diagonal blocks $A_{11},A_{22}$ are. *

Proof 2.10**.**

*See Appendix C. *

This result establishes the fact stated previously, that when $\phi$ is the ReLU, an upper-triangular matrix $A\in\mbox{\rm WP}(\phi)$ if and only if $\mathop{\bf diag}(A)<\mathbf{1}$. A similar result holds with the lower block-triangular matrix

[TABLE]

where $A_{12}\in{\mathbb{R}}^{n_{1}\times n_{2}}$ is arbitrary. It is possible to extend this result to activation maps $\phi$ that satisfy the block Lipschitz continuity (BLIP) condition, in which case we need to assume that the partition of $A$ into blocks is consistent with that of $\phi$. As seen later, this feature arises naturally when composing implicit models from well-posed blocks.

Theorem 2.11** (Well-posedness of block-triangular matrices, blockwise activation).**

Assume that the matrix $A$ can be written as

[TABLE]

*with $A_{ii}\in{\mathbb{R}}^{n_{i}\times n_{i}}$, $i=1,2$, and $\phi$ acts blockwise accordingly, in the sense that there exist two maps $\phi_{1},\phi_{2}$ such that $\phi((z_{1},z_{2}))=(\phi_{1}(z_{1}),\phi_{2}(z_{2}))$ for every $z_{i}\in{\mathbb{R}}^{n_{i}}$, $i=1,2$. Then $A$ is well-posed with respect to $\phi$ if and only if for $i=1,2$, $A_{ii}$ is well-posed with respect to $\phi_{i}$. *

Using the above results, we can preserve well-posedness of implicit models via composition. For example, given two models with matrix parameters $(A_{i},B_{i},C_{i},D_{i})$ and activation functions $\phi_{i}$, $i=1,2$, we can consider a “cascaded” prediction rule:

[TABLE]

The above rule can be represented as Eq. 1a, with $x=(x_{2},x_{1})$, $\phi((z_{2},z_{1}))=(\phi_{2}(z_{2}),\phi_{1}(z_{1}))$ and

[TABLE]

Due to Theorem 2.11, the cascaded rule is well-posed for the componentwise map with values $\phi(z_{1},z_{2})=(\phi_{1}(z_{1}),\phi_{2}(z_{2}))$ if and only if each rule is.

A similar result holds if we put two or more well-posed models in parallel, and do a (weighted) sum the outputs. With the above notation, setting $\hat{y}(u_{1},u_{2})=\hat{y}_{1}(u_{1})+\hat{y}_{2}(u_{2})$ leads to a new implicit model that is also well-posed. Other possible connections include concatenation: $\hat{y}(u)=(\hat{y}_{1}(u),\hat{y}_{2}(u))$, and affine transformations (a special case of cascade connection where one of the systems has no activation). We leave the details to the reader.

In both cascade and parallel connections, the triangular structure of the matrix $A$ of the composed system ensures that the PF sufficient condition for well-posedness is satisfied for the composed system if and only if it holds for each sub-system.

Multiplicative connections are in general are not Lipschitz-continuous, unless the inputs are bounded. Precisely, consider two activation maps $\phi_{i}$ that are Lipschitz-continuous with constant $\gamma_{i}$ and are bounded, with $|\phi_{i}(v)|\leq c_{i}$ for every $v$, $i=1,2$; then, the multiplicative map

[TABLE]

is Lipschitz-continuous with respect to the $l_{1}$-norm, with constant $\gamma:=\max\{c_{2}\gamma_{1},c_{1}\gamma_{2}\}$. Such connections arise in the context of attention units in neural networks, which use (bounded) activation maps such as $\tanh$.

Finally, feedback connections are also possible. Consider two well-posed implicit systems:

[TABLE]

Now let us connect them in a feedback connection: the combined system is described by the implicit rule Eq. 1, where $u_{1}=u+y_{2}$, $u_{2}=y_{1}=y$. The feedback system is also an implicit model of the form (1), with appropriate matrices $(A,B,C,D)$, and activation map acting blockwise: $\phi(z_{1},z_{2})=(\phi_{1}(z_{1}),\phi_{2}(z_{2}))$ and state $(x_{1},x_{2})$. In the simplified case when $D_{1}=D_{2}=0$, the feedback connection has the model matrix

[TABLE]

Note that the connection is not necessarily well-posed when, even if $A_{1}\in\text{WP}(\phi_{1}),A_{2}\in\text{WP}(\phi_{2})$.

2.5 Scaling implicit models

Assume that the activation map $\phi$ is componentwise non-expansive (CONE) and positively homogeneous, as is the ReLU, or its leaky version. Consider an implicit model of the form Eq. 1, and assume it satisfies the PF sufficient condition for well-posedness of Theorem 2.2: $\lambda_{\rm pf}(|A|)\leq\kappa$, where $0\leq\kappa<1$ is given. Then, there is another implicit model with the same activation map $\phi$, matrices $(A^{\prime},B^{\prime},C^{\prime},D^{\prime})$, which has the same prediction rule (i.e. $y^{\prime}(u)=y(u),\forall u$), and satisfies $\|A^{\prime}\|_{\infty}<1$.

This result is a direct consequence of the following expression of the PF eigenvalue as an optimally scaled $l_{\infty}$-norm, known as the Collatz-Wielandt formula, see [39, p. 666]:

[TABLE]

When the eigenvalue is simple, the optimal scaling vector is positive: $s>0$, and the new model matrices are obtained by diagonal scaling:

[TABLE]

where $S=\mathop{\bf diag}(s)$, with $s>0$ a Perron-Frobenius eigenvector of $|A|$.

In a training problem, this result allows us to consider the convex constraint $\|A\|_{\infty}<1$ in lieu of its Perron-Frobenius eigenvalue counterpart. This result also allows us to rescale any given implicit model, such as one derived from deep neural networks, so that the norm condition is satisfied; we will exploit this in our robustness analyses Section 4.

3 Implicit models of deep neural networks

A large number of deep neural networks can be modeled as implicit models. Our goal here is to show how to build a well-posed implicit model for a given neural network, assuming the activation map satisfies the componentwise non-expansiveness (CONE), or the more general block Lipschitz-continuity (BLIP) condition, as detailed in 2.1. Thanks to the composition rules of Section 2.4, it suffices to model individual layers, since a neural network is just a cascade connection of such layers. The block Lipschitz structure then emerges naturally as the result of composing the layers.

We will find that the resulting models always have a strictly (block) upper triangular matrix $A$, which automatically implies that these models are well-posed; in fact the equilibrium equation can be simply solved via backwards substitution. In turn, models with strictly upper triangular structure also naturally satisfy the PF sufficient condition for well-posedness: for example, in the case of a componentwise non-expansive map $\phi$, the matrix $|A|$ arising in Theorem 2.2 is also strictly upper triangular, and therefore all of its eigenvalues are zero. A similar result also holds for the case when $\phi$ is block Lipschitz continuous, as defined in Section 2.1.

3.1 Basic operations

Activation at the output

It is common to have an activation at the output level, as in the rule

[TABLE]

with $\phi_{o}$ the output activation function. We can represent this rule as in Eq. 1, by introducing a new state variable: with $\tilde{\phi}:=(\phi_{o},\phi)$,

[TABLE]

The rule is well-posed with respect to $\tilde{\phi}$ if and only if the matrix $A$ is, with respect to $\phi$.

Bias terms and affine rules

The affine rule

[TABLE]

where $d\in{\mathbb{R}}^{q}$, $b\in{\mathbb{R}}^{n}$, is handled by simply appending a $1$ at the end of the input vector $u$.

A related transformation is useful with activation functions $\phi$, such as the sigmoid, that do not satisfy $\phi(0)=0$. Consider the rule

[TABLE]

Defining $\tilde{\phi}(\cdot):=\phi(\cdot)-\phi(0)$ and using the state vector $\tilde{x}:=x+\phi(0)$, we may represent the above rule as

[TABLE]

Again, the rule is well-posed if and only if the matrix $A$ is.

Batch normalization

Batch normalization consists in including in the prediction rule a normalization step using some estimates of input mean $\hat{u}$ and variance ${\sigma}>0$ that are based on a batch of training data. The parameters $\hat{u},{\sigma}$ are given at test time. This step is a simple affine rule:

[TABLE]

3.2 Fully connected feedforward neural networks

Consider the following prediction rule, with $L>1$ fully connected layers:

[TABLE]

Here $W_{l}\in{\mathbb{R}}^{n_{l+1}\times n_{l}}$ and $\phi_{l}\>:\>{\mathbb{R}}^{n_{l+1}}\rightarrow{\mathbb{R}}^{n_{l+1}}$, $l=1,\ldots,L$, are given weight matrices and activation maps, respectively. We can express the above rule as Eq. 1a, with $x=(x_{L},\ldots,x_{1})$, and

[TABLE]

and with an appropriately defined blockwise activation function $\phi$, defined as operating on an $n$-vector $z=(z_{L},\ldots,z_{1})$ as $\phi(z)=(\phi_{L}(z_{L}),\ldots,\phi_{1}(z_{1}))$. Due to the strictly upper triangular structure of $A$, the system is well-posed, irrespective of $A$; in fact the equilibrium equation $x=\phi(Ax+Bu)$ is easily solved via backward block substitution, which corresponds to a simple forward pass through the network.

3.3 Convolutional layers and max-pooling

A single convolutional layer can be represented as a linear map: $y=Du$, where $u$ is the input, $D$ is a matrix that represents the (linear) convolution operator, with a “constant-along-diagonals”, Toeplitz-like structure. For example a 2D convolution with a 2D kernel $K$ takes a $3\times 3$ matrix $U$ and produces a $2\times 2$ matrix $Y$. With

[TABLE]

the convolution can be represented as $y=Du$, with $y,u$ vectors formed by stacking the rows of $U,Y$ together, and

[TABLE]

Often, a convolutional layer is combined with a max-pooling operation. The latter forms a down-sample of an image, which is a smaller image that contains the largest pixel values of specific sub-areas of the original image. Such an operation can be represented as

[TABLE]

In the above, $p=qh$, and the matrices $B_{j}\in{\mathbb{R}}^{h\times p}$, $j\in[q]$, are used to select specific sub-areas of the original image. In the example of Fig. 5, the number of pixels selected in each area is $h=4$, the output dimension is $q=4$, that of the input is $p=qh=16$; vectorizing images row by row:

[TABLE]

where $I_{2}$ is the $2\times 2$ identity matrix.

Define the mapping $\phi\>:\>{\mathbb{R}}^{n}\rightarrow{\mathbb{R}}^{n}$, where $n=p$, as follows. For a $p$-vector $z$ decomposed in $q$ blocks $(z_{1},\ldots,z_{q})$, we set $\phi(z_{1},\ldots,z_{q})=(\displaystyle\max(z_{1}),\ldots,\max(z_{q}),0,\ldots,0)$. (The padded zeroes are necessary in order to make sure the input and output dimensions of $\phi$ are the same). We obtain the implicit model

[TABLE]

Here $C$ is used to select the top $q$ elements,

[TABLE]

The Lipschitz constant of the max-pooling activation map $\phi$, with respect to the $l_{\infty}$-norm, is $1$, hence it verifies the BLIP condition of Section 2.1.

3.4 Residual nets

A building block in residual networks involves the relationship illustrated in Fig. 7.

Mathematically, a residual block combines two linear operations, with non-linearities in the middle and the end, and adds the input to the resulting output:

[TABLE]

where $W_{1},W_{2}$ are weight matrices of appropriate size. The above is a special case of the implicit model Eq. 1: defining the blockwise map $\phi(z_{2},z_{1})=(\phi_{2}(z_{2}),\phi_{1}(z_{1}))$,

[TABLE]

Fig. 7 displays the model matrix $A$ for a $20$-layer residual network. Convolutional layers appear as matrix blocks with Toeplitz (constant along diagonal) structure; residual unit correspond to the straight lines on top of the blocks. The network only uses the ReLU map, except for the last layer, which uses a softmax map.

3.5 Recurrent units

Recurrent neural nets (RNNs) can be represented in an “unrolled” form as shown in [28], which is the perspective we will consider here. The input to an RNN block is a sequence of vectors $\{u_{1},\cdots,u_{T}\}$, where for every $t\in[T]$, $u_{t}\in{\mathbb{R}}^{p}$. At each time step, the network takes in a input $u_{t}$ and the previous hidden state $h_{t-1}$ to produce the next hidden state $h_{t}$; the hidden state $h_{t}$ defines the state space or “memory” of the network. The rule can be written as,

[TABLE]

Then, equations Eq. 13 can be expressed as an implicit model  Eq. 10, with $x=(h_{T},\ldots,h_{0})$, $u=(u_{T},\dots,u_{1})$, and weight matrices

[TABLE]

where $A,B$ are strictly upper block triangular and share the same block diagonal sub-matrices, $W_{H}$ and $W_{I}$ respectively, shared across all hidden states and inputs.

Note that the above approach leads to matrices $A,B,C,D$ with a special structure; during training and test time, it is possible to exploit that structure in order to avoid “unrolling” the recurrent layers.

3.6 Multiplicative units: LSTM, attention mechanisms and variants

Some deep learning network architectures use multiplicative units. As mentioned in Section 2.4, multiplication between variables in not Lipschitz-continuous in general, unless the inputs are bounded. Luckily, most of the multiplicative units used in practice involve bounded inputs.

In the case of Long Short-Term Memory (LSTM) [55] or gated recurrent units (GRUs) [15], a basic building block involves the product of two input variables after each one passes through a bounded non-linearity. Precisely, the output takes the form

[TABLE]

where $u_{1},u_{2}\in{\mathbb{R}}$, and $\phi_{1},\phi_{2}$ are both bounded (scalar) non-linearities.

Attention models [7] use componentwise vector multiplication, usually involving a softmax operation Eq. 4:

[TABLE]

where $\phi_{1}$ is a bounded, componentwise Lipschitz-continuous activation map, such as the sigmoid; and $W$ is a matrix of weights. The map $\phi$ is then Lipschitz-continuous with respect to the $l_{1}$-norm. As shown in Section 2.2, we can formulate these multiplicative units as well-posed implicit models.

4 Robustness

In this section, our goal is to analyze the robustness properties of a given implicit model Eq. 1a. We seek to bound the state, output and loss function, under unknown-but-bounded inputs. This robustness analysis task is of interest in itself for diagnosis or for generating adversarial attacks. It will also inform our choice of appropriate penalties or constraints in the training problem. We assume that the activation map is Blockwise Lipschitz, and that the matrix $A$ of the implicit model satisfies the sufficient conditions for well-posedness outlined in Theorem 2.7.

4.1 Input uncertainty models

We assume that the input vector is uncertain, and only known to belong to a given set ${\cal U}\subseteq{\mathbb{R}}^{p}$. Our results apply to a wide variety of sets ${\cal U}$; we will focus on the following two cases.

A first case corresponds to a box:

[TABLE]

Here, the $p$-vector $\sigma_{u}>0$ is a measure of componentwise uncertainty affecting each data point, while $u^{0}$ corresponds to a vector of “nominal” inputs. The following variant limits the number of changes in the vector $u$:

[TABLE]

where $\mathop{\bf Card}$ denotes the cardinality (number of non-zero components) in its vector argument, and $k<p$ is a given integer.

4.2 Box bounds on the state vector

Assume first that $\phi$ is a CONE map, and that the input is only known to belong to the box set ${\cal U}^{\rm box}$ Eq. 22. We seek to find componentwise bounds on the state vector $x$, of the form $|x-x^{0}|\leq\sigma_{x}$, with $x$ and $x^{0}$ the unique solution to $\xi=\phi(A\xi+Bu)$ and $\xi=\phi(A\xi+Bu^{0})$ respectively, and $\sigma_{x}>0$. We have

[TABLE]

which shows in particular that

[TABLE]

hence, provided $\|A\|_{\infty}<1$, we have:

[TABLE]

With the cardinality constrained uncertainty set ${\cal U}^{\rm card}$ Eq. 23, we obtain

[TABLE]

with $e_{i}$ the $i$-th unit vector in ${\mathbb{R}}^{n}$, and $s_{k}$ the sum of the top $k$ entries in a vector.

We may refine the analysis above, with a “box” (componentwise) bound.

Theorem 4.1** (Box bound on state, CONE map).**

Assuming that $\phi$ is componentwise non-expansive, and $\lambda_{\rm pf}(|A|)<1$. then, $I-|A|$ is invertible, and

[TABLE]

Proof 4.2**.**

*See Appendix D. *

Note that the box bound can be computed via $(I-|A|)\sigma_{x}=|B|\sigma_{u}$, as the limit point of the fixed-point iteration

[TABLE]

which converges since $\lambda_{\rm pf}(|A|)<1$.

The case when $\phi$ is block Lipschitz (BLIP) involves the matrix of norms $N(A)$ defined in Eq. 7, as well as a related matrix defined for the input matrix $B$. Decomposing $B$ into blocks $B=(B_{li})_{l\in[L],i\in[p]}$, with $B_{li}\in{\mathbb{R}}^{n_{l}}$, $l\in[L]$, we define the $L\times p$ matrix of norms

[TABLE]

Theorem 4.3** (Box bound on the vector norms of the state, BLIP map).**

Assuming that $\phi$ is BLIP, and the corresponding sufficient well-posedness condition of Theorem 2.7 is satisfied. Define $\Gamma:=\mathop{\bf diag}(\gamma)$, where $\gamma$ is the vector of Lipschitz constants. Then, $I-\Gamma N(A)$ is invertible, and

[TABLE]

*where the vector of norms function $\eta(\cdot)$ is defined in Eq. 2. *

Proof 4.4**.**

*See Appendix D. *

4.3 Bounds on the output and the sensitivity matrix

The above allows us to analyze the effect of effect of input noise on the output vector $y$. Let us assume that the function $\phi$ satisfies the CONE (componentwise non-expansiveness) condition Eq. 3. In addition, we assume that the stronger condition for well-posedness, $\|A\|_{\infty}<1$, is satisfied. (we can always rescale the model so as to ensure that property, provided $\lambda_{\rm pf}(|A|)<1$.) For the implicit prediction rule Eq. 1, we then have

[TABLE]

where

[TABLE]

The above shows that the prediction rule is Lipschitz-continuous, with a constant bounded above by $\kappa$. This result motivates the use of the $\|\cdot\|_{\infty}$ norm as a penalty on model matrices $A,B,C,D$, for example via a convex penalty that bounds the Lipschitz constant,

[TABLE]

We can refine this analysis with the following theorem.

Theorem 4.5** (Box bound on the output, CONE map).**

Assuming that $\phi$ is componentwise non-expansive, and $\lambda_{\rm pf}(|A|)<1$. Then, $I-|A|$ is invertible, and

[TABLE]

where the (non-negative) matrix

[TABLE]

*is a “sensitivity matrix” of the implicit model with a CONE map. *

The sensitivity matrix is a way to visualize the input-output properties of a given model. Fig. 8 shows the sensitivity matrix for a given neural network for classification that has $q=10$ classes, with each row represented as an image. We observe that some classes have higher sensitivity values than others, which can be used to focus an attack on the model.

Our refined analysis suggests a “natural” penalty to use during the training phase in order to improve robustness, namely $\|S\|_{\infty}$.

As before, the approach can be generalized to the case when $\phi$ is block Lipschitz (BLIP). Decomposing $C$ into column blocks $C=(C_{1},\ldots,C_{L})$, with $C_{l}\in{\mathbb{R}}^{q\times n_{l}}$, $l\in[L]$, we define the matrix of (dual) norms

[TABLE]

where $p^{*}_{l}:=1/(1-1/p_{l})$, $l\in[L]$. Also recall the corresponding matrix of norms related to $B$ Eq. 27.

Theorem 4.6** (Box bound on the output, BLIP map).**

Assuming that $\phi$ is a BLIP map, and that the sufficient condition for well-posedness $\lambda_{\rm pf}(\Gamma N(A))<1$ is satisfied, with $\Gamma:=\mathop{\bf diag}(\gamma)$ containing the Lipschitz constants associated with $\phi$. Then, $I-\Gamma N(A)$ is invertible, and

[TABLE]

where the (non-negative) $q\times p$ matrix

[TABLE]

*is a “block sensitivity matrix” of the implicit model with a BLIP map. *

Proof 4.7**.**

*See Appendix D. *

4.4 Worst-case loss function

In this section, we assume that the map $\phi$ satisfies the CONE property Eq. 3. We seek to find bounds on the loss function evaluated between a given “target” $y\in{\mathbb{R}}^{p}$ and a prediction $\hat{y}$, using the “box” bounds on the state and output of Section 4.2 and Section 4.3.

Consider the case of squared Euclidean loss function, that is:

[TABLE]

The worst-case loss can be computed based on the box bound Eq. 31, as

[TABLE]

We can also consider the cross entropy loss function, that is:

[TABLE]

where $y\geq 0$, $\mathbf{1}^{\top}y=1$ is given.

Assume that $D=0$ for simplicity, so that $\hat{y}^{0}:=Cx^{0}$ is the nominal output. Using the norm bound Eq. 24, we have $\|x-x^{0}\|_{\infty}\leq\rho$, for an appropriate $\rho>0$. The corresponding worst-case loss is:

[TABLE]

with $c_{i}^{\top}$ the $i$-th row of $C$. The above may be hard to compute, but we can work with a bound, based on evaluating the maximum above for each of the two terms independently. We obtain

[TABLE]

Since $y\geq 0$ and $\mathbf{1}^{\top}y=1$, we have $\|C^{\top}y\|_{1}\leq\|C\|_{\infty}$, and we obtain the bound

[TABLE]

which again involves the norm $\|C\|_{\infty}$ as a penalty. We can refine the analysis, based on the box bound Eq. 31:

[TABLE]

4.5 LP relaxation for CONE maps

The previous bounds does not provide a direct way to generate an adversarial attack, that is, a feasible point $u\in{\cal U}$ that has maximum impact on the state. In some cases it may be possible to refine our previous box bounds via an LP relaxation, which has the advantage of suggesting a specific adversarial attack. Here we restrict our attention to the ReLU activation: $\phi(z)=\max(z,0)=z_{+}$, which is a CONE map.

We consider the problem

[TABLE]

where $f_{i}$’s are arbitrary functions. Setting $f_{i}(\xi)=(\xi-x_{i}^{0})^{2}$, $i\in[n]$, leads to the problem of finding the largest discrepancy (measured in $l_{2}$-norm) between $x$ and $x^{0}$; setting $f_{i}(\xi)=-\xi$, $i\in[n]$, finds the minimum $l_{1}$-norm state. By construction, our bound can only improve on the previous state bound, in the sense that

[TABLE]

Our result is expressed for general sets ${\cal U}$, based on the support function $\sigma_{\cal U}$ with values for $b\in{\mathbb{R}}^{p}$ given by

[TABLE]

Note that this function depends only on the convex hull of the set ${\cal U}$. In the case of the box set ${\cal U}^{\rm box}$ given in Eq. 22, there is a convenient closed-form expression:

[TABLE]

Likewise for the cardinality-constrained set ${\cal U}^{\rm card}$ defined in LABEL:{eq:U-def-card}, we have

[TABLE]

where $s_{k}$ is the sum of the top $k$ entries of its vector argument—a convex function.

The only coupling constraint in Eq. 33 is the affine equation, which suggests the following relaxation.

Theorem 4.8** (LP bound on the state).**

An upper bound on the objective of problem Eq. 33 is given by

[TABLE]

where $s_{\cal U}$ is the support function defined in Eq. 34, and $g_{i}$, $i\in[n]$, are the convex functions

[TABLE]

If the functions $g_{i}$, $i\in[n]$ are closed, we have the bidual expression

[TABLE]

*where $g_{i}^{*}$ is the conjuguate of $g_{i}$, $i\in[n]$. *

Proof 4.9**.**

*See Appendix E. *

Consider for example $f_{i}(\xi)=c_{i}\xi$, $i\in[n]$, where $c\in{\mathbb{R}}^{n}$ is given: the problem is

[TABLE]

Let $\underline{x}:=(x^{0}-\sigma_{x})_{+}$, $\overline{x}:=x^{0}+\sigma_{x}$, so that the constraint $x\geq 0$, $|x-x^{0}|\leq\sigma_{x}$ writes $\underline{x}\leq x\leq\overline{x}$. We have

[TABLE]

This leads to the dual

[TABLE]

For every $i\in[n]$, the function $g_{i}$ is closed, with conjugate given by

[TABLE]

We obtain a “natural” relaxation:

[TABLE]

We consider a variant of the previous problem, where we seek a sparse adversarial attack: the cardinality of changes in the input is constrained in the set ${\cal U}^{\rm card}$.

[TABLE]

where ${\cal U}^{\rm card}$ is the set Eq. 23, with $k\in[p]$ given. The relaxation expresses as

[TABLE]

The bidual writes

[TABLE]

where $\mathop{\bf diag}(\cdot)$ is a diagonal matrix formed with the entries of its vector argument. The above provides a heuristic for a sparse adversarial attack, based on any vector $u$ that is optimal for the above.

4.6 SDP relaxations for maximal state discrepancy with ReLU activation

Assume again that $\phi$ is the ReLU, and consider the box uncertainty set ${\cal U}^{\rm box}$ Eq. 22. Our goal is to find the largest state discrepancy, measured in the Euclidean norm:

[TABLE]

Based on a quadratic representation of $\phi$, we can arrive at a SDP relaxation of (36). We first use the fact that for any two $n$-dimensional vectors $x,z$:

[TABLE]

where $\odot$ is the Hadamard (component-wise) product. Next, we use the fact that any bound the form $l\leq z-z_{0}\leq u$ where $z,l,u$ are vectors can be rewritten as

[TABLE]

These two reformulations allow us to represent (36) as a non-convex QCQP:

[TABLE]

We can derive an upper bound on the value of this QP by solving the Lagrange dual of the problem, leading to a semidefinite program. However, it is a well known result of quadratic optimization (for example see Appendix B of [10]) that the bidual of a non-convex QCQP simply becomes the canonical rank relaxed version of (36) when we redefine the objective in terms of a new variable $Z=zz^{\top}$, where $z=[x,\;u,\;1]^{\top}$. Making this substitution and relaxing the equality constraint into a semidefinite constraint, the bidual of (36) becomes

[TABLE]

where used a rank relaxation involving

[TABLE]

By construction, (41) is convex. One immediate reason why we wish to solve the bidual as opposed to the Lagrange dual of (36) is that we are able to get a candidate state $x$ and attack $u$, which can be viewed, respectively, as the perturbed state resulting from an adversarial attack on our network. This attack can be retrieved from the bidual in different ways. The first one would be to look at the optimal $u$ and the second would be to do a rank one decomposition of $X$ and extract $u$ from it.

5 Sparsity and Model Compression

In this section, we examine the role of sparsity and low-rank structure in implicit deep learning, specifically in the model parameter matrix

[TABLE]

We assume that the activation map $\phi$ is a CONE map, as defined in Section 2.1. Most of our results can be generalized to BLIP maps.

Since a CONE map acts componentwise, the prediction rule Eq. 1a is invariant under permutations of the state vector, in the sense that, for any $n\times n$ permutation matrix, the matrix $\mathop{\bf diag}(P,I)M\mathop{\bf diag}(P^{\top},I)$ represents the same prediction rule as $M$ given above.

Various kinds of sparsity of $M$ can be encouraged in the training problem, with appropriate penalties. For example, we can use penalties that encourage many elements in $M$ to be zero; the advantage of such “element-wise” sparsity is, of course, computational, since sparsity in matrices $A,B,C,D$ will allow for computational speedups at test time. Another interesting kind of sparsity is rank sparsity, which refers to the case when model matrices are low-rank.

Next, we examine the benefits of row- (or, column-) sparsity, which refers to the fact that entire rows (or, columns) of a matrix are zero. Note that column sparsity in a matrix $N$ can be encouraged with a penalty in the training problem, of the form

[TABLE]

where $\alpha>1$. Row sparsity can be handled via ${\cal P}(N^{\top})$.

5.1 Deep feature selection

We may use the implicit model to select features. Any zero column in the matrix $(B^{\top},D^{\top})^{\top}$ means that the corresponding element in an input vector does not play any role in the prediction rule. We may thus use a column-norm penalty in the training problem, in order to encourage such a sparsity pattern:

[TABLE]

with $\alpha>1$.

5.2 Dimension reduction via row- and column-sparsity

Assume that the matrix $A$ is row-sparse. Without loss of generality, using permutation invariance, we can assume that $M$ writes

[TABLE]

where $A_{11}$ is square of order $n_{1}<n$. We can then decompose $x$ accordingly, as $x=(x_{1},x_{2})$ with $x_{1}\in{\mathbb{R}}^{n_{1}}$, and the above implies $x_{2}=\phi(B_{2}u)$. The prediction rule for an input $u\in{\mathbb{R}}^{p}$ then writes

[TABLE]

The rule only involves $x_{1}$ as a true hidden feature vector. In fact, the row sparsity of $A$ allows for a computational speedup, as we simply need to solve a fixed-point equation for the vector with reduced dimensions, $x_{1}$.

Further assume that $(A,B)$ is row-sparse. Again without loss of generality we may put $M$ in the above form, with $B_{2}=0$. Then the prediction rule can be written

[TABLE]

This means that the dimension of the state variable can be fully reduced, to $n_{1}<n$. Thus, row sparsity of $(A,B)$ allows for a reduction in the dimension of the prediction rule.

Assume now that the matrix $A$ is column-sparse. Without loss of generality, using permutation invariance, we can assume that $M$ writes

[TABLE]

where $A_{11}$ is square of order $n_{1}<n$. We can then decompose $x$ accordingly, as $x=(x_{1},x_{2})$ with $x_{1}\in{\mathbb{R}}^{n_{1}}$. The above implies that the prediction rule for an input $u\in{\mathbb{R}}^{p}$ writes

[TABLE]

Thus, column-sparsity allows for a computational speedup, since $x_{2}$ can be directly expressed as closed-form function of $x_{1}$.

Now assume that $(A^{\top},C^{\top})^{\top}$ is column-sparse. Again without loss of generality we may put $M$ in the above form, with $C_{2}=0$. We obtain that the prediction rule does not need $x_{2}$ at all, so that the computation of the latter vector can be entirely avoided. This means that the dimension of the state variable can be fully reduced, to $n_{1}<n$. Thus, column sparsity of $(A^{\top},C^{\top})^{\top}$ allows for a reduction in the dimension of the prediction rule.

To summarize, row or column sparsity of $A$ allows for a computational speedup; if the corresponding rows of $B$ (resp. columns of $C$) are zero, then the prediction rule involves only a vector of reduced dimensions.

5.3 Rank sparsity

Assume that the matrix $A$ is rank $k\ll n$, and that a corresponding factorization is known: $A=LR^{\top}$, with $L,R\in{\mathbb{R}}^{n\times k}$. In this case, for any $p$-vector $u$, the equilibrium equation $x=\phi(Ax+Bu)$ can be written as $x=\phi(Lz+Bu)$, where $z:=R^{\top}x$. Hence, we can obtain a prediction for a given input $u$ via the solution of a low-dimensional fixed-point equation in $z\in{\mathbb{R}}^{k}$:

[TABLE]

It can be shown that, when $\phi$ is a CONE map, the above rule is well-posed if $\lambda_{\rm pf}(|R|^{T}|L|)<1$. Once a solution $z$ is found, we simply set the prediction to be $\hat{y}(u)=C\phi(Lz+Bu)+Du$.

At test time, if we use fixed-point iterations to obtain our predictions, then the computational savings brought about by the low-rank representation of $A$ can be substantial, with a per-iteration cost going from $O(n^{2})$, to $O(kn)$ if we use the above.

5.4 Model error analysis

The above suggests to replace a given model matrix $A$ with a low-rank or sparse approximation, denoted $A^{0}$. The resulting state error can be then bounded, as follows. Assume that $|A-A^{0}|\leq E$, where $E\geq 0$ is a known upper bound on the componentwise error in $A$. The following theorem, proved in Appendix F, provides relative error bounds on the state, provided the perturbed system satisfies the well-posedness condition $\lambda_{\rm pf}(|A|)<1$. As before, we denote by $x^{0},x$ the (unique) solutions to the unperturbed and perturbed equilibrium equations $\xi=\phi(A^{0}\xi+Bu)$ and $\xi=\phi(A\xi+Bu)$, respectively.

Theorem 5.1** (Effect of errors in $A$).**

Assuming that $\phi$ is a CONE map, and that $\lambda_{\rm pf}(|A^{0}+E|)<1$. Then:

[TABLE]

Proof 5.2**.**

*See Appendix F. *

6 Training Implicit Models

6.1 Setup

We are now given an input data matrix $U=[u_{1},\ldots,u_{m}]\in{\mathbb{R}}^{p\times m}$ and response matrix $Y=[y_{1},\ldots,y_{m}]\in{\mathbb{R}}^{q\times m}$, and seek to fit a model of the form Eq. 1a, with $A$ well-posed with respect to $\phi$, which we assume to be a BLIP map. We note that the rule Eq. 1a, when applied to a collection of inputs $(u_{i})_{1\leq i\leq m}$, can be written in matrix form, as

[TABLE]

where $U=[u_{1},\ldots,u_{m}]\in{\mathbb{R}}^{p\times m}$, and $\hat{Y}(U)=[\hat{y}(u_{1}),\ldots,\hat{y}(u_{m})]\in{\mathbb{R}}^{q\times m}$.

We consider a training problem of the form

[TABLE]

In the above, ${\cal L}$ is a loss function, assumed to be convex in its second argument, and ${\cal P}$ is a convex penalty function, which can be used to enforce a given (linear) structure (such as, $A$ strictly upper block triangular) on the parameters, and/or encourage their sparsity.

In practice, we replace the well-posedness condition by the sufficient condition of Theorem 2.7. As argued in Section 2.5, the latter can be further replaced without loss of generality with the easier-to-handle (convex) constraint; in the case of CONE maps, this condition is $\|A\|_{\infty}\leq\kappa$, where $\kappa\in(0,1)$ is a hyper-parameter.

Examples of loss functions

For regression tasks, we may use the squared Euclidean loss: for $Y,\hat{Y}\in{\mathbb{R}}^{q\times m}$,

[TABLE]

For multi-class classification, a popular loss is a combination of negative cross-entropy with the soft-max: for two $q$-vectors $y,\hat{y}$, with $y\geq 0$, $y^{\top}\mathbf{1}=1$, we define

[TABLE]

We can extend the definition to matrices, by summing the contribution to all columns, each corresponding to a data point: for $Y,Z\in{\mathbb{R}}^{q\times m}$,

[TABLE]

where both the $\log$ and the exponential functions apply componentwise.

Examples of penalty functions

Via an appropriate definition of ${\cal P}$, we can make sure that the model matrix $A$ is well-posed with respect to $\phi$, either imposing an upper triangular structure for $A$, or via an $l_{\infty}$-norm constraint $\|A\|_{\infty}<1$, or a Perron-Frobenius eigenvalue constraint. Note that in the case of the CONE maps such as the ReLU, due to scale invariance seen in Section 2.5, we can always replace a Perron-Frobenius eigenvalue constraint with a $l_{\infty}$-norm constraint.

Beyond well-posedness, the penalty can be used to encourage desired properties of the model. For robustness, the convex penalty Eq. 30 can be used, provided we also enforce $\|A\|_{\infty}<1$. We may also use an $l_{\infty}$-norm bound on the sensitivity matrices $S$ appearing in Theorem 4.5 or Theorem 4.6. For feature selection, an appropriate penalty may involve the block norms Eq. 48; sparsity of the model matrices can be similarly handled with ordinary $l_{1}$-norm penalties on the elements of the model matrix $M=(A,B,C,D)$.

Fenchel divergence formulation

Problem (50) can be equivalently written

[TABLE]

where $F_{\phi}$ is the so-called Fenchel divergence adapted to $\phi$ [24], applied column-wise to matrix inputs. In the case of the ReLU activation, for two given vectors $x,z$ of the same size, we have

[TABLE]

with $\odot$ the componentwise multiplication. We can then define $F_{\phi}(X,Z)$ with $X,Z$ two matrix inputs having the same number of columns, by summing over these columns. As seen in [24], a large number of popular activation maps can be expressed similarly.

We may also consider a relaxed form:

[TABLE]

where $\lambda>0$ is a $n$-vector parameter.

6.2 Gradient Methods

Assuming that $\phi$ is a CONE map for simplicity, we consider the problem

[TABLE]

where $\kappa<1$ is given.

We can solve Eq. 53 using (stochastic) projected gradient descent. The basic idea is to update the model parameters $M=(A,B,C,D)$ using stochastic gradients, by differentiating through the equilibrium equation. It turns out that this differentiation requires solving an equilibrium equations involving a matrix variable, which can be very efficiently solved via fixed-point iterations. More precisely, the main effort in computing the gradient for a mini-batch of size say $b$ consists in solving matrix equations in a $n\times b$ matrix $V$:

[TABLE]

where the columns of $L$ contains gradients of the loss with respect to $\hat{y}$, and $\Psi$ is a matrix whose columns contain the derivatives of the activation map, both evaluated at a specific training point. Due to the fact taht $A$ satisfies the PF sufficient condition for well-posedness with respect to $\phi$, the equation above has a unique solution; the matrix $V$ can be computed as the limit point of the convergent fixed-point iterations

[TABLE]

The proof of the following result is given in Appendix G.

Theorem 6.1** (Computing gradients).**

*Consider an implicit model ($\phi$, $A,B,C,D$), with $\phi$ a differentiable BLIP map. If $A$ is well-posed with respect to $\phi$, the gradients with respect to the model parameter matrices $A,B,C$ and $D$ exist, and are uniquely given via the solution of a fixed-point equation, as detailed in Appendix G. *

In order to handle the well-posedness constraint $\|A\|_{\infty}\leq\kappa$, the projected gradient method requires a projection at each step. This step corresponds to a sub-problem of the form:

[TABLE]

with matrix $A^{0}$ given. The above problem is decomposable across rows, and can be efficiently solved using (vectorized) bisection, as detailed in Appendix H.

In some cases, it is desirable to obtain a (block) sparse model matrix $M$. One way to efficiently achieve sparsity is the use of a conditional gradient (Frank-Wolfe) algorithm that typically only adds a few non-zeros elements to the model matrix $M$ at each step. When updating $A$ for example, considering $B,C,D$ fixed for simplicity, the method requires solving a problem of the form

[TABLE]

where $G$ is a gradient matrix. The problem can be solved in closed-form.

6.3 Block-coordinate descent methods

Block-coordinate descent methods use cyclic updates, optimizing one matrix variable at a time, fixing all the other variables. These methods are particularly interesting when the updates require solving convex problems. For instance, considering the training problem Eq. 53, given $X$, the problem is convex in the model matrix $M$. Then, given the model matrix $M$, finding $X$ consists in a feasibility problem that can be solved with fixed-point iterations.

In the case of the Fenchel divergence problem relaxation, all the updates involve solving convex problems, as shown in [24], as the Fenchel divergence $F_{\phi}$, for most of activation maps, is bi-convex in its two arguments—meaning that given the first argument fixed, it is convex in the second and vice-versa. We refer to [24, 51] for more on Fenchel divergences in the context of implicit deep learning and neural networks. We illustrate and give more detailed examples of such methods in the next section.

7 Numerical experiments

7.1 Learning real nonlinear functions via regression

We start by illustrating the ability of the gradient method, as presented in Section 6.2, to learn the parameters of the implicit model, with well-posedness enforced via a max-row-sum norm constraint, $\|A\|_{\infty}\leq 0.5$. We aim at learning a real function, as an example we focus on

[TABLE]

We select the input $u_{i}$, $i\in\{1,\cdots,m\}$ uniformly at random between $-5$ and $5$ with $m=200$; we add a random noise to the output, $y(u)=f(u)+w$ with $w$ taken uniformly at random between $-1$ and $1$, hence the standard deviation for $y(u)$ is $1/\sqrt{(}3)\simeq 0.57$. We consider an implicit model of order $n=75$. We learn the parameters of the model by doing only two successive block updates: first, we update $(A,B)$ using stochastic projected gradient descent, the gradient being obtained with the implicit chain rule described in Appendix G. The RMSE across iterations for this block-update is shown in Fig. 10. After this first update we achieve a RMSE of $1.77$. Second, we update $(C,D)$ using linear regression. After this update we achieve a RMSE of $0.56$. For comparison purposes, we also train a neural network with $3$ hidden layers of width $n/3=25$ using ADAM, mini-batches, and a tuned learning rate. We run Adam until convergence. We get a RMSE $=0.65$, which is slightly above that of our implicit model. This first simple experiment shows the ability to fit nonlinear functions with implicit models as illustrated in Fig. 9.

7.2 Comparison with neural networks

In that section we compare the performance of implicit models with that of neural networks. Experiments on both synthetic datasets and real datasets are conducted. The experiment results show that implicit model has the potential of matching or exceeding the performance of neural networks in various settings. To simplify the experiments, we do not apply any specific network structure or regularization during training. When it comes to training neural networks we will always use ADAM with a grid-search tuned step-size. Similarly, we will always use stochastic projected gradient descent as detailed in Appendix G. The choice of the number of hidden features $n$ for implicit models is always aligned with the neural network architecture for fair comparison as explained in Section 3.2.

7.2.1 Synthetic datasets

We first consider two synthetic classification datasets: one generated from a neural network and an other from an implicit model. For each dataset, we then aim at fitting both an implicit model and a neural network. Each data point in the datasets contains an input $u\in\mathbb{R}^{5}$ ($p=5$) and output $y\in\mathbb{R}^{2}$ ($q=2$). The dataset is generated as follows.

We chose the input datapoints $u_{i}\in{\mathbb{R}}^{5}$ i.i.d. by sampling each entry independently uniformly between $-1$ and $1$. The output data $y_{i}\in{\mathbb{R}}^{2}$ is the one-hot representation of the argmax of the output from the generating model, whose parameters are chosen uniformly at random between $-1$ and $1$.

• •

For a neural network generating model, we us a fully connected 3-layer (5-5-5-2) feedforward neural network with ReLU activation.

• •

For an implicit model, we consider an implicit model with $n=10$ to match the number of hidden nodes in the neural network model. To ensure well-posedness, after sampling matrix the $A$, we proceed with a projection on the unit-sized infinity norm ball. (See Appendix H).

For each dataset, we consider 20 training and test datapoints. The size of datasets are kept small to maintain the over-parameterized settings where the number of parameters is larger than the number of data points. Deep learning falls under the such regime [9]. We run 5 separate experiments for both the neural network and implicit model generated dataset with the training model having the same hyperparameters and initialization as the generating models.

We find that the implicit model outperforms neural networks in both synthetic experiments. This may be explained by the increased modeling capacity of the implicit model, given similar parameter size, with respect to its neural network counterpart as mentioned in Section 1.1.

7.2.2 MNIST dataset

In the digit classification dataset MNIST, the input data points are $28\times 28$ pixels images of hand written digits, the output is the corresponding digit label. For training purposes, each image is reshaped into 784 dimensional vectors and normalized before training. There are $m=5\times 10^{4}$ training data points and $10^{4}$ testing data points.

The architecture we use for the neural network as a reference is a three-layer feedforward neural network (784-60-40-10) with ReLU activation. For the implicit model, we set $n=100$. We choose a batchsize of $100$ for both algorithms. The accuracy with respect to iterations is shown in Fig. 14. We observe that the accuracy of the implicit model matches that of its neural network counterpart.

7.2.3 German Traffic Sign Recognition Benchmark

In the German Traffic Sign Recognition Benchmark (GTSRB) [49], the input data points are $32\times 32$ images with rgb channels of traffic signs consisting of 43 classes. Each image is turned into gray-scale before being reshaped into a 1024 dimensional vector and re-scaled to be between 0 and 1. There are 34,799 training data points and 12,630 testing data points.

For reference, we use a (1024-300-100-43) feedforward neural network with ReLU activations. We choose $n=400$ for the implicit model. We choose a batchsize of 100. The accuracy with respect to iterations is shown in Fig. 14.

Similar to what we observed with synthetic datasets, the implicit model is capable of matching and even outperform classical neural networks.

7.3 Adversarial attack

7.3.1 Attack via the sensitivity matrix

Our analysis above highlights the use of sensitivity matrix as a measure for robustness. In this section, we show how sensitivity matrix can be used to generate effective attacks on two public datasets, MNIST and CIFAR-10. We compare our method against commonly used gradient-based attacks [22, 41]. In this experiment, we consider two models: 1) feed-forward network trained on the MNIST dataset (98% clean accuracy) and 2) ResNet-20 [27] trained on the CIFAR-10 dataset (92% clean accuracy).

Figure 15 shows the sensitivity values for a particular class on MNIST and CIFAR-10 dataset. The sensitivity values are obtain from the implicit representation of the feed-forward network and ResNet-20. The 784 and 3072 input dimensions are arranged to correspond to the $28\times 28$ (grey scale) and $32\times 32$ (color) image pixel alignment. Brighter colors correspond to features with a higher impact on the output when perturbed. As a result, the sensitivity matrix can be used to generate adversarial attacks.

We compare our method with commonly used gradient-based attacks. Precisely, for a given function $F$ (prediction rule) learned by a deep neural network, a benign sample $u\in{\mathbb{R}}^{p}$ and the target $y$ associated with $u$, we compute the gradient of the function $F$ with respect to the given sample $u$, $\nabla_{u}F(u,y)$. We then take the absolute value of the gradient as an indication of which input features an adversary should perturb, similar to the saliency map technique [47, 41]. The absolute value of the gradient can be seen as a “local” version of the sensitivity matrix; however, unlike the gradient, the sensitivity matrix does not depend on the input data, making it a more general measurement of robustness for any given model.

Table 1 presents the experimental results of an adversarial attack using the sensitivity matrix and the absolute value of gradient on MNIST and CIFAR-10. For sensitivity matrix attack, we start from perturbing the input features that have the highest values according to the sensitivity matrix. For gradient-based attack, we do the same according to the absolute value of the gradient. We perturb the input features into small random values. Our experiments show that the sensitivity matrix attack is as effective as the gradient-based attack, while being very simple to implement.

Interestingly, our attack does not rely on any input samples that the gradient-based attack needs. An adversary with the model parameters could easily craft adversarial samples using the sensitivity matrix. In the absence of access to the model parameters, an adversary can rely on the principle of tranferrability [35] and train a surrogate model to obtain the sensitivity matrix. Fig. 16 displays adversarial images generated using the sensitivity matrix. An interesting case is to use the sensitivity matrix to generate a sparse attack as seen in Fig. 16.

7.3.2 Attack with LP relaxation for CONE maps

Although one can use sensitivity matrix to generate an effective adversarial example, one may wish to perform a more sophisticated attack by exploiting the weakness of an individual data point. This can be done by considering the LP relaxation (Theorem 4.8), which has the advantage of generating a specific adversarial example for a given input data. The experiment in this section is again on MNIST and CIFAR-10 images. Here, the problem outlined in Eq. 33 is solved by LP relaxation, with the function $f_{i}(\xi)=(\xi-x_{i}^{0})^{2}$. The optimization problem then finds a perturbed image that leads to the largest discrepancy between the perturbed state $x$ and the nominal state $x^{0}$. Fig. 17 shows five example images. The clean images are shown in the top row. The perturbed images generated by the LP relaxation, as shown in the bottom row, appear quite similar to the naked eye; however, the model fails to predict these images correctly. The examples shown here are all correctly predicted by the model when no perturbation is added.

Our framework also allows for sparse adversarial attack by adding a cardinality constraint. Fig. 18 shows three examples of perturbed images under non-sparse and sparse attack. Images on the left are the results of non-sparse attack and those on the right are the results of sparse attack. The model fails to predict the label correctly under both conditions. These results illustrate how the implicit prediction rules can be used to generate powerful adversarial attacks. It is also useful for adversarial training as a large amount of adversarial examples can be generated using the technique and be added back to the training data.

8 Prior Work

8.1 In implicit deep learning

Recent works have considered versions of implicit models, and demonstrated their potential in deep learning. In particular, recent ground-breaking work by Kolter and collaborators [8, 29] demonstrated success of an entirely implicit framework, which they call Deep Equilibrium Models, for the task of sequence modeling. Paper [14] uses implicit methods to solve and construct a general class of models known as neural ordinary differential equations, while [17] uses implicit models to construct a differentiable physics engine that enables gradient-based learning and high sample efficiency. Furthermore, many papers explore the concept of integrating implicit models with modern deep learning methods in a variety of ways. For example, [52] show promise in integrating logical structures into deep learning by incorporating a semidefinite programming (SDP) layer into a network in order to solve a (relaxed) MAXSAT problem; see also [53]. In [1] the authors propose to include a model predictive control as a differentiable policy class for deep reinforcement learning, which can be seen as a kind of implicit architecture. In [2] the authors introduced implicit layers where the activation is the solution of some quadratic programming problem; in [18], the authors incorporate stochastic optimization formulation for end-to-end learning task, in which the model is trained by differentiating the solution of a stochastic programming problem.

8.2 In lifted models

In implicit learning, there is usually no way to express the state variable in closed-form, which makes the task of computing gradients with respect to model parameters challenging. Thus, a natural idea in implicit learning is to keep the state vector as a variable in the training problem, resulting in a higher-dimensional (or, “lifted”) expression of the training problem. The idea of lifting the dimension of the training problem in (non-implicit) deep learning by introducing “state” variables has been studied in a variety of works; a non-extensive list includes [50], [5], [24], [56], [57], [12] and [34]. Lifted models are trained using block coordinate descent methods, Alternating Direction Method of Multipliers (ADMM) or iterative, non-gradient based methods. In this work, we introduce a novel aspect of lifted models, namely the possibility of defining a prediction rule implicitly.

8.3 In robustness analysis

The issue of robustness in deep learning is generating quite a bit of attention, due to the fact that many deep learning models suffer from the lack of robustness. Prior relevant work have demonstrated that deep learning models are vulnerable to adversarial attacks [22, 30, 41, 31]. The vulnerability issue of deep learning models have motivated the study of corresponding defense strategies [37, 42, 43, 23, 46, 54, 16]. However, many of the defense strategies are later shown to be ineffective [6, 11], suggesting the needs for the theoretical understanding of robustness evaluations for deep learning model. In this work, we formalize the robustness analysis of deep learning via the lens of the implicit model. A large number of deep learning architectures can be modeled using implicit prediction rules, making our robustness evaluation a versatile analysis tool.

In [43], the authors arrive at a similar formulation (which we discuss below) but do not explore the quality of the attacks that can be extracted by solving (41).

The SDP-based formulation given in Section 4.6 parallels the formulation in [43] where the authors represent interval sets ($l\leq x\leq u$) as the quadratic constraint $x\circ x\leq(l+u)\circ x-l\circ u$ and work on the rank-relaxed version of their problem to arrive at a SDP. Upon further inspection, our formulation is similar to, and therefore extends, the formulation proposed in [43] for feed-forward neural networks. Indeed, for a neural network, $A$ is block strictly upper triangular and $(I-|A|)^{-1}$ has a closed form expression; in the case of three layers we have:

[TABLE]

leading to

[TABLE]

The difference between our proposed SDP and the SDP in [43] is how the input-uncertainty constraints are dealt with at every layer. In [43], the authors strengthen their formulation by reducing the size of the uncertainty set at each layer. That is, they propagate the input uncertainty and bound the size of the uncertainty at each layer of the neural network (see Sections 5.1 and 5.2 of [43]). As written, (41) does not include these constraints and hence can be tightened by adding extra box-constraints as is done in the aforementioned reference. This allows to extend, to implicit models, the results that were derived for traditional feed-forward networks.

8.4 In sparsity, compression and deep feature selection

Sparsity and compression, which are well understood in classical settings, have found their place in deep learning and are an active branch of research. Early work in pruning dates back to as early as the 90s [33, 26] and has since gained interest. In [48], the authors showed that by randomly dropping units (i.e. increasing the sparsity level of the network or compressing the network) reduces overfitting and improved the generalization performance of networks. Recently, more sophisticated ways of pruning networks have been proposed, in an effort to reduce the overall size of the model, while retaining or accepting a modest decrease in accuracy: a non-extensive list of works include [58, 40, 25, 45, 3, 32, 13, 36, 20]. Sparsity and model compression is very closely tied with other techniques to reduce model size. This includes, but is not limited to, low-rank matrix factorization, group sparsity regularization, quantization techniques and low precision networks. While this sub-branch of research is far from being complete, sparsity and compression offer promise in the way of low memory footprints and generalization performance for deep networks.

9 Concluding Remarks

We conclude with a few perspectives for future research on implicit models.

9.1 Cousins of implicit models

Implicit models rely on a representation of the prediction rule where the linear operations are clearly separated from the (parameter-free) nonlinear ones, leading to the block-diagram in Fig. 1.

This idea is strongly reminiscent of earlier representations arising in systems and control theory. Perhaps the most famous example lies with linear time-invariant (LTI) systems, where the idea is essentially equivalent to the state-space representation of transfer functions [4]. In that context, the block-diagram representation involves an integrator (in continuous or discrete time) in place of the static nonlinear map $\phi$. The idea is also connected to the linear-fractional representation (or, transformation, LFT) arising in uncertain systems [19], where now the linear operations actually represent a dynamic system (that is, the matrix $M$ is an LTI system itself), and the map $\phi$ is replaced by an uncertain matrix of parameters or LTI system.

These connections raise the prospect of a unified theory tackling deep networks inside the loop of a dynamical system, where such networks can be composed (via feedback connections) with LTI or more generally uncertain systems. With the machinery of Linear Matrix Inequalities or the more general Integral Quadratic Constraints framework [38], one can develop rigorous analyses performed on systems with deep networks in the loop, so that bounds on, say, stability margins can be computed.

9.2 Minimal representations

The prediction rule $u\rightarrow\hat{y}(u)$ provided by a given implicit model can be in general obtained by many alternative models—scaling the model matrices via the diagonal scaling given in Eq. 8 illustrates this non-unicity. An interesting theoretical study would focus on obtaining representations that are minimal from the point of view of size, precisely the order $n$, which would extend the well-established theory of minimal representations of LTI systems.

Acknowledgments

We would like to acknowledge useful comments and suggestions by Mert Pilanci and Zico Kolter.

Appendix A Proof of Theorem 2.2

Let $b\in{\mathbb{R}}^{n}$. We first prove the existence of a solution $\xi\in{\mathbb{R}}^{n}$ to the equation $\xi=\phi(A\xi+b)$. Consider the Picard iteration Eq. 6. We have for every $t\geq 1$:

[TABLE]

which implies that for every $t,\tau\geq 0$:

[TABLE]

where

[TABLE]

Here we exploited the fact that, due to $\lambda_{\rm pf}(|A|)<1$, $I-|A|$ is invertible, and the series above converges. Since $\lim_{t\rightarrow 0}|A|^{t}=0$, we obtain that $x(t)$ is a Cauchy sequence, hence it has a limit point, $x_{\infty}$. By continuity of $\phi$ we further obtain that $x_{\infty}=\phi(Ax_{\infty}+b)$, which establishes the existence of a solution.

To prove unicity, consider $x^{1},x^{2}\in{\mathbb{R}}^{n}_{+}$ two solutions to the equation. Using the hypotheses in the theorem, we have, for any $k\geq 1$:

[TABLE]

The fact that $|A|^{k}\rightarrow 0$ as $k\rightarrow+\infty$ then establishes unicity.

Appendix B Proof of Theorem 2.7

We follow the steps taken in the proof of Theorem 2.2. Let $b\in{\mathbb{R}}^{n}$. Our first step is to establish that for the Picard iteration Eq. 6, we have

[TABLE]

for an appropriate non-negative matrix $M\geq 0$. Here, $\eta$ is a vector of norms, as defined in Eq. 2. For every $l\in[L]$, $t\geq 1$:

[TABLE]

which establishes the desired bound, with $M:=\mathop{\bf diag}(\gamma)N(A)$, where $\gamma$ is the vector of Lipschitz constants, and $N(\cdot)$ is the $L\times L$ matrix of induced norms Eq. 7.

Assume now that $\lambda_{\rm pf}(M)<1$, as posited in the Theorem. Then $I-M$ is invertible. We proceed to prove existence by showing that the sequence of Picard iterates is Cauchy: for every $t,\tau\geq 0$,

[TABLE]

where

[TABLE]

The above proves the existence.

To prove unicity, consider $x^{1},x^{2}\in{\mathbb{R}}^{n}_{+}$ two solutions to the equation. Using the hypotheses in the theorem, we have, for any $k\geq 1$:

[TABLE]

The fact that $M^{k}\rightarrow 0$ as $k\rightarrow+\infty$ then establishes unicity.

Appendix C Proof of Theorem 2.9

Express the equation $x=\phi(Ax+b)$ as

[TABLE]

where $b=(b_{1},b_{2})$, $x=(x_{1},x_{2})$, with $b_{i}\in{\mathbb{R}}^{n_{i}}$, $x_{i}\in{\mathbb{R}}^{n_{i}}$, $i=1,2$. Here, since $\phi$ acts componentwise, we use the same notation $\phi$ in the two equations.

Now assume that $A_{11}$ and $A_{22}$ are well-posed with respect to $\phi$. Since $A_{22}$ is well-posed for $\phi$, the second equation has a unique solution $x_{2}^{*}$; plugging $x_{2}=x_{2}^{*}$ into the second equation, and using the well-posedness of $A_{11}$, we see that the first equation has a unique solution in $x_{1}$, hence $A$ is well-posed.

To prove the converse direction, assume that $A$ is well-posed. The second equation above must have a unique solution $x_{2}^{*}$, irrespective to the choice of $b_{2}$, hence $A_{22}$ must be well-posed. To prove that $A_{11}$ must be well-posed too, set $b_{2}=0$, $b_{1}$ arbitrary, leading to the system

[TABLE]

Since $A_{22}$ is well-posed for $\phi$, there is a unique solution $x_{2}^{*}$ to the second equation; the first equation then reads $x_{1}=\phi(A_{11}x_{1}+b_{1}+A_{12}x_{2}^{*})$. It must have a unique solution for any $b_{1}$, hence $A_{11}$ is well-posed.

Appendix D Proofs of Theorem 4.1, Theorem 4.3 and Theorem 4.6

First let us prove Theorem 4.1. Due to the non-expansivity of the CONE map $\phi$:

[TABLE]

where $b:=|B|\sigma_{u}\geq 0$. Since $\lambda_{\rm pf}(|A|)<1$, the matrix $I-|A|$ is invertible and its inverse is componentwise non-negative; for any $n$-vector $b\geq 0$, we have

[TABLE]

Further, any $n$-vector $\sigma$ such that $\sigma\leq|A|\sigma+b$, we have, for every integer $N$:

[TABLE]

Since $\lambda_{\rm pf}(|A|)<1$, we have $|A|^{N}\rightarrow 0$ as $N\rightarrow+\infty$, which yields $\sigma\leq\bar{\sigma}$. Applying this to $\sigma=|x-x^{0}|$ yields the desired result.

Now turn to the proof of Theorem 4.3. For every $l\in[L]$:

[TABLE]

which establishes the desired bound.

The result of Theorem 4.6 is obtained similarly. For given $i\in[q]$, we have

[TABLE]

Appendix E Proof of Theorem 4.8

We have

[TABLE]

which establishes the first part of the theorem. If we further assume that the functions $g_{i}$, $i\in[n]$ are closed, strong duality holds, so that

[TABLE]

where $g_{i}^{*}$ is the conjugate of $g_{i}$, $i\in[n]$.

Appendix F Proof of Theorem 5.1

We have

[TABLE]

Applying a technique similar to that employed in the proof of Theorem 4.1, we obtain the desired relative error bounds.

Appendix G Gradient Equations

In this section, we detail gradient computations with respect to the model parameter matrix $M=(A,B,C,D)$, in the context of the training problem 53. We assume that the map $\phi$ is BLIP with parameters $n_{l}$, $p_{l}$, $\gamma_{l}$, $l\in[L]$, and that $A$ satisfies the corresponding PF well-posedness condition with respect to $\phi$; we also assume that the latter is differentiable.

For simplicity, we first consider mini-batches of size $1$ (that is, $X=x\in\mathbb{R}^{n}$ and $U=u\in\mathbb{R}^{p}$), and we define $\hat{y}=Cx+Du$, $z=Ax+Bu$. We wish to calculate

[TABLE]

The difficult part are the terms $\nabla_{A}\mathcal{L}$ and $\nabla_{B}\mathcal{L}$ due to the presence of the equilibrium equality constraint. We deal with this via implicit differentiation. In order to establish the existence of gradients we will need the following lemma.

Lemma G.1**.**

Assume that the map $\phi$ is a BLIP map with parameters $n_{l}$, $p_{l}$, $\gamma_{l}$, $l\in[L]$, and that $A$ is well-posed with respect to $\phi$. Define the block-diagonal matrix $\Phi:=\frac{\partial\phi(z)}{\partial z}$, then the equation in the $n\times n$ matrix $G$: $G=\Phi(AG+I)$ has a unique solution, which can be computed as the limit point of the recursion

[TABLE]

Proof G.2**.**

*The result follows from the fact that if $\phi$ is BLIP with parameters $n_{l}$, $p_{l}$, $\gamma_{l}$, $l\in[L]$, then $\Phi=\mathop{\bf diag}(\Phi_{1},\ldots,\Phi_{L})$, where $\|\Phi_{l}\|_{n_{l}}<1$, $l\in[L]$ and Theorem 2.7. *

Calculating $\nabla_{A}\mathcal{L}$

For a given index pair $(j,k)$, we have

[TABLE]

We calculate $\nabla_{z}\mathcal{L}$ via implicit differentiation:

[TABLE]

where $\Phi:=\frac{\partial\phi(z)}{\partial z}$ is a block diagonal matrix. Thanks to Lemma G.1, the matrix $G:=(I-\Phi A)^{-1}\Phi$ exists and can be obtained via fixed-point iterations 57. Note that the gradient of the loss function $\nabla_{\hat{y}}\mathcal{L}$ can be easily computed, and we have

[TABLE]

Calculating the other gradients

The other gradients are easily computed, as seen next, where $(j,k)$ denotes a generic index pair. From the above it follows that

[TABLE]

where $\nabla_{z}\mathcal{L}$ is given in 58. Likewise,

[TABLE]

and

[TABLE]

Finally, the gradient with respect to input $u$ is obtained as

[TABLE]

where $\nabla_{z}\mathcal{L}$ is given in 58.

To summarize, we have obtained the following closed form evaluation of gradients for ($A,B,C,D$):

[TABLE]

where $\nabla_{\hat{y}}\mathcal{L}$ can be easily computed and $\nabla_{z}\mathcal{L}$ is given by the following equation:

[TABLE]

Let $v:=\nabla_{z}\mathcal{L}$, we have

[TABLE]

which is equivalent to

[TABLE]

Again Eq. 59 can be computed via the recursion (54). In practice we use Eq. 59 to compute $\nabla_{z}\mathcal{L}$ without explicitly forming $G$ via Eq. 57 where each iteration would involve a matrix-matrix product instead of a matrix-vector product in Eq. 59.

Mini-batch calculations

We may efficiently compute the gradient corresponding to a whole mini-batch of size greater than $1$. In the case of componentwise maps $\phi$, this is based on expressing the equation (59) as

[TABLE]

where columns of $L$ contains gradients of the loss with respect to $\hat{y}$, and $\Psi$ is a matrix that contains the derivatives of the activation, one column corresponding to the (diagonal) elements of $\Phi$, where $\Phi$ is defined in the previous section.

Appendix H Projection on $l_{\infty}$ Matrix Norm Ball

We address problem Eq. 55, which we write as

[TABLE]

where $A^{0}\in{\mathbb{R}}^{n\times n}$ is given. The problem is decomposable across the rows of the matrices involved, leading to $n$ sub-problems of the form

[TABLE]

which $a^{0}_{i}\in{\mathbb{R}}^{n}$ the $i$-th row of $A^{0}$.

The problem cannot be solved in closed form, but a bisection method can be applied to the dual:

[TABLE]

where, for $\lambda\geq 0$ given:

[TABLE]

A subgradient of the objective is

[TABLE]

Observe that $p^{*}\geq 0$, hence at optimum:

[TABLE]

The bisection can be initialized with the interval $\lambda\in[0,\lambda^{\rm max}]$.

Returning to the original problem Eq. 55, we see that all the iterations can be expressed in a “vectorized” form, where updates for the different rows of $A$ are done in parallel. The dual variables corresponding to each row are collected in a vector $\lambda\in{\mathbb{R}}^{n}$. We initialize the bisection with a vector interval $[\lambda_{l},\lambda_{u}]$, with $\lambda^{l}=0$, $\lambda^{u}_{i}=\frac{1}{2}\|a^{0}_{i}\|_{2}^{2}/\kappa$, $i\in[n]$. We update the current vector interval as follows:

Set $\lambda=(\lambda_{l}+\lambda_{u})/2$. 2. 2.

Form a vector $g(\lambda)$ containing the sub-gradients corresponding to each row, evaluated at $\lambda_{i}$, $i\in[n]$:

[TABLE] 3. 3.

For every $i\in[n]$, reset $\lambda^{u}_{i}=\lambda_{i}$ if $g_{i}(\lambda)>0$, $\lambda^{l}_{i}=\lambda_{i}$ if $g_{i}(\lambda)\leq 0$.