# Monotone, Bi-Lipschitz, and Polyak-Łojasiewicz Networks

Ruigang Wang
The University of Sydney
<EMAIL_ADDRESS>
Krishnamurthy (Dj) Dvijotham
Google DeepMind
<EMAIL_ADDRESS>
Ian R. Manchester
The University of Sydney
<EMAIL_ADDRESS>
###### Abstract

This paper presents a new _bi-Lipschitz_ invertible neural network, the
BiLipNet, which has the ability to control both its _Lipschitzness_ (output
sensitivity to input perturbations) and _inverse Lipschitzness_ (input
distinguishability from different outputs). The main contribution is a novel
invertible residual layer with certified strong monotonicity and
Lipschitzness, which we compose with orthogonal layers to build bi-Lipschitz
networks. The certification is based on incremental quadratic constraints,
which achieves much tighter bounds compared to spectral normalization.
Moreover, we formulate the model inverse calculation as a three-operator
splitting problem, for which fast algorithms are known. Based on the proposed
bi-Lipschitz network, we introduce a new scalar-output network, the PLNet,
which satisfies the Polyak-Łojasiewicz condition. It can be applied to learn
non-convex surrogate losses with favourable properties, e.g., a unique and
efficiently-computable global minimum.

##  1 Introduction

In many applications, it is desirable to learn neural networks with
_guaranteed input-output behaviors_ , i.e., certain properties that are
ensured by design. For example, Lipschitz-bounded networks have proven to be
beneficial to stabilization of generative adversarial network (GAN) training
Arjovsky et al., (2017); Gulrajani et al., (2017), certifiable robustness
against adversarial attacks Tsuzuku et al., (2018); Revay et al., (2020);
Singla and Feizi, (2021); Pauli et al., (2021); Prach and Lampert, (2022);
Araujo et al., (2023); Wang and Manchester, (2023) and robust reinforcement
learning Russo and Proutiere, (2021).

Model | inv. Lip. | Lip. | loss
---|---|---|---
i-ResNet | 0.80 | 4.69 | 0.2090
i-DenseNet | 0.82 | 4.66 | 0.2091
BiLipNet | 0.11 | 9.97 | 0.0685
Optimal | 0.10 | 10.0 | 0.0677

Figure 1: Fitting a step function, which is not Lipschitz, with certified
$(0.1,10)$-Lipschitz models. Compared to the analytically-computed optimum,
the proposed BiLipNet achieves much tighter bounds than models based on
spectral normalization.

Another input-output property – _invertibility_ has received much attention in
the deep learning literature since the introduction of _normalizing flows_
Dinh et al., (2015) for probability-density learning. Invertible neural
networks (INNs) have been applied in applications such as generative modeling
Dinh et al., (2017); Kingma and Dhariwal, (2018), probabilistic inference
Bauer and Mnih, (2019); Ward et al., (2019); Louizos and Welling, (2017),
solving inverse problems Ardizzone et al., (2018) and uncertainty estimation
Liu et al., (2020). A common way to construct INNs is to compose simple
invertible transformations and sophisticated invertible layers, including
coupling flows Dinh et al., (2017); Kingma and Dhariwal, (2018), auto-
regressive models Huang et al., (2018); De Cao et al., (2020); Ho et al.,
(2019), invertible residual layers Chen et al., (2019); Behrmann et al.,
(2019), monotone networks Ahn et al., (2022), and neural ordinary differential
equations Grathwohl et al., (2019), see also in the surveys Papamakarios et
al., (2021); Kobyzev et al., (2020).

However, Behrmann et al., (2021) observed that commonly-used INNs suffer from
exploding inverses and are thus prone to becoming numerically non-invertible.
This observation motivates the input-output property of _bi-Lipschitzness_. A
layer ${\mathcal{F}}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is said to be
bi-Lipschitz with bound of $(\mu,\nu)$, or simply $(\mu,\nu)$-Lipschitz, if
the following inequalities hold for all $x,x^{\prime}\in\mathbb{R}^{n}$:

$\mu\|x-x^{\prime}\|\leq\|{\mathcal{F}}(x)-{\mathcal{F}}(x^{\prime})\|\leq\nu\|x-x^{\prime}\|,$
(1)

where $\|\cdot\|$ is the Euclidean norm. The bound $\nu$ controls the output
sensitivity to input perturbations while $\mu$ controls the input
distinguishability from different outputs. Here $\mu$ is also called the
_inverse Lipschitz_ bound as ${\mathcal{F}}^{-1}$ exists and is
$1/\mu$-Lipschitz. For applications involving uncertainty estimation, this
prevents the output from being unnecessarily invariant to the semantically
meaningful changes in the input Liu et al., (2020). The ratio $\tau:=\nu/\mu$
is called _distortion_ Liang et al., (2023), which is the upper bound of the
condition number of the Jacobian matrix of ${\mathcal{F}}$. A larger
distortion implies more expressive flexibility in the model.

In this paper we argue that the bi-Lipschitz property is also useful for
learning of surrogate loss (or reward) functions. Given some input/output
pairs of a loss function, the objective is to learn a function which matches
the observed data and is “easy to optimize” in some sense. This problem
appears in many areas, including Q-learning with continuous action spaces, see
e.g. Gu et al., (2016); Amos et al., (2017); Ryu et al., (2019), learning
reward models in inverse reinforcement learning Arora and Doshi, (2021), and
data-driven surrogate losses for engineering process optimization Cozad et
al., (2014); Misener and Biegler, (2023). An important contribution was the
input convex neural network (ICNN) Amos et al., (2017). However, the
requirement of input convexity could be too strong in many applications.

###  1.1 Contribution

  * •

We propose a novel strongly monotone and Lipschitz residual layer of the form
${\mathcal{F}}(x)=\mu x+{\mathcal{H}}(x)$. For the nonlinear block
${\mathcal{H}}$, we introduce a new architecture – _feed-through network_
(FTN), which takes a multi-layer perceptron (MLP) as its backbone and adds
connections from each hidden layer to the input and output variables. For deep
networks, this architecture can improve the model expressivity without
suffering from vanishing gradients. We establish the end-to-end _strong
monotonicity_ (a special case of inverse Lipschitzness) and Lipschitzness for
${\mathcal{F}}$ via the integral quadratic constraint (IQC) framework
Megretski and Rantzer, (1997), which has recently been shown to enable very
tight Lipschitz bounds for network analysis Fazlyab et al., (2019) and
synthesis Araujo et al., (2023); Wang and Manchester, (2023), see also the toy
example in Figure 1.

  * •

We formulate the model inversion ${\mathcal{F}}^{-1}$ as a three-operator
splitting problem, which admits a numerically efficient solver Davis and Yin,
(2017).

  * •

We introduce a new scalar-output network
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}$, which we call a Polyak-Łojasiewicz
network (or PLNet) since it satisfies the condition of the same name Polyak,
(1963); Lojasiewicz, (1963). It consists of a bi-Lipschitz network composed
with a quadratic potential, and automatically satisfies favourable properties
for surrogate loss learning, in particular a unique global optimum which is
efficiently computable.

###  1.2 Related work

#### Bi-Lipschitz invertible layer.

In the literature, there are two types of invertible layers closely related to
our models. The first is the invertible residual layer
${\mathcal{F}}(x)=x+{\mathcal{H}}(x)$ Chen et al., (2019); Behrmann et al.,
(2019), where the nonlinear block ${\mathcal{H}}$ is a shallow network with
Lipschitz bound of $c<1$. In Perugachi-Diaz et al., (2021), ${\mathcal{H}}$ is
further extended to a deep MLP. It is easy to show that ${\mathcal{F}}$ is
$(1-c)$-inverse Lipschitz and $(1+c)$-Lipschitz. In both cases, the Lipschitz
regularization is via spectral normalization Miyato et al., (2018), which we
observe to be very conservative (see Figure 1). Alternatively, bi-Lipschitz
${\mathcal{F}}$ can be defined by an implicit equation Lu et al., (2021); Ahn
et al., (2022). However, these require an iterative solver for both the
forward and inverse model inference. In contrast, our model has an explicit
forward pass and iterative solution is only required for the inverse.

#### IQC-based Lipschitz estimation and training.

In Fazlyab et al., (2019), the IQC framework of Megretski and Rantzer, (1997)
was first applied to obtain accurate Lipschitz bound estimation of deep
networks with slope-restricted activations. Direct (i.e. unconstrained)
parameterizations based on IQC were were proposed in Revay et al., (2020) for
deep equilibrium networks, in Araujo et al., (2023) for residual networks, for
deep MLPs and CNNs in Wang and Manchester, (2023), and recurrent models in
Revay et al., (2023). It was pointed out by Havens et al., (2023) that many
recent Lipschitz model parameterizations Meunier et al., (2022); Prach and
Lampert, (2022); Araujo et al., (2023); Wang and Manchester, (2023) are
special cases of Revay et al., (2020). In a recent work Pauli et al., (2024),
the IQC-based Lipschitz estimation was recently extended to more general
activations such as GroupSort and MaxMin. All of these are for one-sided
(upper) Lipschitzness, whereas our work applies the IQC framework for
monotonicity and bi-Lipschitzness.

#### Surrogate constraint learning.

Liang et al., (2023) considers the problem of learning a surrogate constraint
set. The basic idea is to (i) learn a minimum _distortion_ homeomorphic
mapping between the non-convex constraint set ${\mathbb{X}}$ and a unit ball
${\mathbb{B}}$ using an INN ${\mathcal{F}}$, and then (ii) perform a simple
bisection operation to find $y\in{\mathbb{B}}$ so that its inverse
${\mathcal{F}}^{-1}(y)$ is feasible w.r.t. ${\mathbb{X}}$ with minor
distortion-induced optimality loss. The distortion estimation in their work is
based on data samples while our work offers certified distortion bounds.

##  2 Preliminaries

Figure 2: This figure depicts the possible ranges of $\Delta
y={\mathcal{F}}(x^{\prime})-{\mathcal{F}}(x)$ on $\mathbb{R}^{2}$ for a given
$\Delta x=x^{\prime}-x$. The ring (blue area) is for $(\mu,\nu)$-Lipschitz
${\mathcal{F}}$ while the half moon (red area) is for a $\mu$-strongly
monotone and $\nu$-Lipschitz ${\mathcal{F}}$. The largest angle between
$\Delta x$ and $\Delta y$ satisfies $\cos\alpha=\tau^{-1}$ with $\tau=\nu/\mu$
as the distortion.

We give some definitions for a mapping
${\mathcal{F}}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$.

######  Definition 2.1.

${\mathcal{F}}$ is said to be $\mu$-strongly monotone with $\mu>0$ if for all
$x,x^{\prime}\in\mathbb{R}^{n}$ we have

$\left\langle{{\mathcal{F}}(x)-{\mathcal{F}}(x^{\prime})},{x-x^{\prime}}\right\rangle\geq\mu\|x-x^{\prime}\|^{2},$
(2)

where $\left\langle{\cdot},{\cdot}\right\rangle$ is the Euclidean inner
product: $\left\langle{a},{b}\right\rangle=a^{\top}b$. ${\mathcal{F}}$ is
monotone if the above condition holds for $\mu=0$.

######  Definition 2.2.

${\mathcal{F}}$ is said to be $\nu$-Lipschitz with $\nu>0$ if

$\|{\mathcal{F}}(x)-{\mathcal{F}}(x^{\prime})\|\leq\nu\|x-x^{\prime}\|,\quad\forall
x,x^{\prime}\in\mathbb{R}^{n}.$

${\mathcal{F}}$ is said to be $\mu$-inverse Lipschitz with $\mu>0$ if

$\|{\mathcal{F}}(x)-{\mathcal{F}}(x^{\prime})\|\geq\mu\|x-x^{\prime}\|,\quad\forall
x_{1},x_{2}\in\mathbb{R}^{n}.$

${\mathcal{F}}$ is said to be bi-Lipschitz with $\nu\geq\mu>0$, or simply
$(\mu,\nu)$-Lipschitz, if it is $\mu$-inverse Lipschitz and $\nu$-Lipschitz.

For any $(\mu,\nu)$-Lipschitz mapping ${\mathcal{F}}$, its inverse
${\mathcal{F}}^{-1}$ is well-defined and $(1/\nu,1/\mu)$-Lipschitz Yeh,
(2006). By the Cauchy–Schwarz inequality, strong monotonicity is a special
case of inverse Lipschitzness, see Figure 2. A notable difference between
monotonicity and bi-Lipschitzness is their composition behaviour. Given two
bi-Lipschitz mappings ${\mathcal{F}}_{1},{\mathcal{F}}_{2}$, their composition
${\mathcal{F}}={\mathcal{F}}_{2}\circ{\mathcal{F}}_{1}$ is also bi-Lipschitz
with bound of $(\mu_{1}\mu_{2},\nu_{1}\nu_{2})$ where $(\mu_{1},\nu_{1})$ and
$(\mu_{2},\nu_{2})$ are the bi-Lipschitz bounds of ${\mathcal{F}}_{1}$ and
${\mathcal{F}}_{2}$, respectively. However, given two strongly monotone
${\mathcal{F}}_{1},{\mathcal{F}}_{2}$ with monotonicity bounds
$\mu_{1},\mu_{2}$, the composition
${\mathcal{F}}={\mathcal{F}}_{2}\circ{\mathcal{F}}_{1}$ does not need to be
strongly monotone. However, it is still $\mu_{1}\mu_{2}$-inverse Lipschitz.

######  Definition 2.3.

${\mathcal{F}}$ satisfies a distortion bound $\tau$ with $\tau\geq 1$ if
${\mathcal{F}}$ is $(\mu,\nu)$-Lipschitz with $\tau=\nu/\mu$.

For example, an invertible affine mapping ${\mathcal{F}}(x)=Px+q$ has a
distortion bound of $\kappa$, the condition number of $P$. An orthogonal
mapping (i.e., $P^{\top}P=I$) has the smallest possible model distortion
$\tau=1$. Distortion bounds satisfy a composition property, i.e., if
${\mathcal{F}}_{1},{\mathcal{F}}_{2}$ have distortion bound of
$\tau_{1},\tau_{2}$, then ${\mathcal{F}}_{2}\circ{\mathcal{F}}_{1}$ satisfies
a distortion bound of $\tau_{1}\tau_{2}$. Both ${\mathcal{F}}$ and
${\mathcal{F}}^{-1}$ have the same distortion.

##  3 Monotone and bi-Lipschitz Networks

In this section we first present the main contribution: construction of
$\mu$-strongly monotone and $\nu$-Lipschitz residual layers of the form
${\mathcal{F}}(x)=\mu x+{\mathcal{H}}(x)$. Then, we use it to construct bi-
Lipschitz networks by composition with orthogonal layers.

###  3.1 Feed-through network

Figure 3: The proposed invertible residual network ${\mathcal{F}}(x)=\mu
x+{\mathcal{H}}(x)$ where the nonlinear block ${\mathcal{H}}$ is a feed-
through network, whose hidden layers are directly connected to the input and
output.

For the nonlinear block ${\mathcal{H}}$, we introduce a network architecture,
called _feed-through network_ (FTN), which takes an MLP as its backbone and
then connects each hidden layer to input and output variables, see Figure 3.
To be specific, the residual layer ${\mathcal{F}}(x)=\mu x+{\mathcal{H}}(x)$
can be written as

$\begin{split}z_{k}&=\sigma(W_{k}z_{k-1}+U_{k}x+b_{k}),\;z_{0}=0\\\ y&=\mu
x+\sum_{k=1}^{L}Y_{k}z_{k}+b_{y}\end{split}$ (3)

where $z_{k}\in\mathbb{R}^{m_{k}}$ are the hidden variables,
$U_{k},W_{k},Y_{k}$ and $b_{k},b_{y}$ are the learnable weights and biases,
respectively. Throughout the paper we assume that the activation $\sigma$ is a
scalar nonlinearity with slope restricted in $[0,1]$, which is satisfied
(possibly with rescaling) by common activation functions such as ReLU,
smoothed/leaky ReLU, tanh, and sigmoid.

######  Remark 3.1.

FTN contains both short paths $x\rightarrow y$ and $x\rightarrow
z_{i}\rightarrow y$ preventing vanishing gradients, as well as long paths
$x\rightarrow z_{i}\rightarrow\cdots\rightarrow z_{j}\rightarrow y$ improving
model expressivity (see Figure 3).

###  3.2 SDP conditions for monotonicity and Lipschitzness

We rewrite ${\mathcal{F}}$ in a compact form:

$z=\sigma(Wz+Ux+b),\quad y=\mu x+Yz+b_{y}$ (4)

where
$z=\bigl{[}\,z_{1}^{\top}\;\cdots\;z_{L}^{\top}\,\bigr{]}^{\top},\,b=\bigl{[}\,b_{1}^{\top}\;\cdots\;b_{L}^{\top}\,\bigr{]}^{\top}$,
and

$\begin{split}W&=\begin{bmatrix}0&\\\ W_{2}&0\\\ &\ddots&\ddots\\\
&&W_{L}&0\end{bmatrix},\quad U=\begin{bmatrix}U_{1}\\\ U_{2}\\\ \vdots\\\
U_{L}\end{bmatrix},\\\
Y&=\begin{bmatrix}Y_{1}&Y_{2}&\cdots&Y_{L}\end{bmatrix}.\end{split}$

We first establish strong monotonicity and Lipschitzness for ${\mathcal{F}}$
via semidefinite programming (SDP) conditions.

######  Theorem 3.2.

${\mathcal{F}}$ is $\mu$-strongly monotone and $\nu$-Lipschitz if there exists
a $\Lambda\in{\mathbb{D}}_{+}^{m}$, where ${\mathbb{D}}_{+}^{m}$ is the set of
positive diagonal matrices, such that the following conditions hold:

$Y=U^{\top}\Lambda,\quad 2\Lambda-\Lambda
W-W^{\top}\Lambda\succeq\frac{2}{\gamma}Y^{\top}Y$ (5)

where $\gamma=\nu-\mu>0$.

######  Remark 3.3.

The above conditions are obtained by applying the IQC theory Megretski and
Rantzer, (1997) to (4).

###  3.3 Model parameterization

Let $\Theta$ be the set of all $\theta=\\{\Lambda,U,W,Y\\}$ such that
Condition (5) holds. Since it is generally not scalable to train a model with
SDP constraints, we instead construct a direct parameterization, i.e. both
unconstrained and complete:

######  Definition 3.4.

A direct parameterization of a constraint set $\Theta$ is a differentiable
mapping ${\mathcal{M}}$ such that for any $\phi\in\mathbb{R}^{N}$ we have
${\mathcal{M}}(\phi)\in\Theta$, and the image of $\mathbb{R}^{N}$ maps onto
$\Theta$, i.e. ${\mathcal{M}}(\mathbb{R}^{N})=\Theta$.

A direct parameterization allows us to replace a constrained optimization over
$\theta\in\Theta$ with an unconstrained optimization over
$\phi\in\mathbb{R}^{N}$, using standard first-order optimization algorithms
such as SGD or ADAM Kingma and Ba, (2015).

We now construct the direct parameterization, see Appendix A for details.
First, we introduce the free parameter

$\phi=\\{d,F_{k}^{a},F_{k}^{b},F^{q},F^{\star}\\},\quad k=1,\ldots,L$

where $d\in\mathbb{R}^{m}$, $F_{k}^{a}\in\mathbb{R}^{m_{k}\times m_{k}}$,
$F_{k}^{b}\in\mathbb{R}^{m_{k-1}\times m_{k}}$, $F^{q}\in\mathbb{R}^{m\times
n}$ and $F^{\star}\in\mathbb{R}^{n\times n}$ with $m_{0}=0$. Then, we compute
some intermediate variables $\Psi=\mathrm{diag}\bigl{(}e^{\psi}\bigl{)}$ and

$\begin{bmatrix}A_{k}^{\top}\\\
B_{k}^{\top}\end{bmatrix}=\operatorname{Cayley}\left(\begin{bmatrix}F_{k}^{a}\\\
F_{k}^{b}\end{bmatrix}\right),\;\begin{bmatrix}Q\\\
\star\end{bmatrix}=\operatorname{Cayley}\left(\begin{bmatrix}F^{q}\\\
F^{\star}\end{bmatrix}\right)$

where $\operatorname{Cayley}:\mathbb{R}^{n\times
p}\rightarrow\mathbb{R}^{n\times p}$ with $n\geq p$ is defined by

$J=\operatorname{Cayley}\left(\begin{bmatrix}G\\\
H\end{bmatrix}\right):=\begin{bmatrix}(I+Z)(I-Z)^{-1}\\\
-2V(I-Z)^{-1}\end{bmatrix}$ (6)

where $Z=G^{\top}-G+H^{\top}H$. It is easy to verify that $J^{\top}J=I$ for
any $G\in\mathbb{R}^{p\times p}$ and $H\in\mathbb{R}^{(n-p)\times p}$. The
Cayley transformation has been applied to construct orthogonal layers Helfrich
et al., (2018); Li et al., (2020); Trockman and Kolter, (2021). By introducing
$V_{k}=2B_{k}A_{-1}^{\top}$ and $S_{k}=A_{k}Q_{k}-B_{k}Q_{k-1}$ with
$B_{1}=0$, $Q_{0}=0$ and
$Q=\bigl{[}\,Q_{1}^{\top}\;\cdots\;Q_{L}^{\top}\,\bigr{]}^{\top}$, we finally
construct $\theta={\mathcal{M}}(\phi)$ as follows:

$\begin{split}U_{k}&=\sqrt{2\gamma}\Psi_{k}^{-1}S_{k},\;W_{k}=\Psi_{k}^{-1}V_{k}\Psi_{k-1},\\\
Y_{k}&=\sqrt{\frac{\gamma}{2}}S_{k}^{\top}\Psi_{k},\;\Lambda_{k}=\frac{1}{2}\Psi_{k}^{2}.\end{split}$
(7)

######  Proposition 3.5.

The model parameterization ${\mathcal{M}}$ defined in (7) is a direct
parameterization for the set $\Theta$, i.e. all models (4) satisfying
Condition (5).

This means that we can learn the free parameter $\phi$ using first-order
methods without any loss of model expressivity.

The construction is now done, but we note that $\Psi_{k}$ is shared between
layers $k$ and $k+1$. To have a modular implementation, we introduce the
following scaled variables $\hat{z}=\Psi z$ and bias $\hat{b}=\Psi b$. Then,
(4) can be rewritten as follows

$\hat{z}=\hat{\sigma}\bigl{(}V\hat{z}+\sqrt{2\gamma}Sx+\hat{b}\bigr{)},\;y=\mu
x+\sqrt{\gamma/2}S^{\top}\hat{z}+b_{y}$ (8)

where $\hat{\sigma}(x):=\Psi\sigma\left(\Psi^{-1}x\right)$ is a
$(0,1)$-Lipschitz layer with learnable scaling $\Psi$. The weight matrices are

$S=\begin{bmatrix}S_{1}\\\ S_{2}\\\ \vdots\\\ S_{L}\end{bmatrix},\quad
V=\begin{bmatrix}0&\\\ V_{2}&0\\\ &\ddots&\ddots\\\ &&V_{L}&0\end{bmatrix}.$
(9)

###  3.4 Bi-Lipschitz networks

By composing strongly monotone and Lipschitz layers, we can obtain bi-
Lipschitz networks (referred as BiLipNets):

${\mathcal{G}}={\mathcal{O}}_{K+1}\circ{\mathcal{F}}_{K}\circ{\mathcal{O}}_{K}\circ{\mathcal{F}}_{K-1}\circ\cdots\circ{\mathcal{O}}_{2}\circ{\mathcal{F}}_{1}\circ{\mathcal{O}}_{1}$
(10)

where ${\mathcal{O}}_{k}(x)=P_{k}x+q_{k}$ with $P_{k}^{\top}P_{k}=I$ is an
orthogonal layer and ${\mathcal{F}}_{k}$ is a $\mu_{k}$-strongly monotone and
$\nu_{k}$-Lipschitz layer (8). By the composition rule, the above BiLipNet is
$(\mu,\nu)$-Lipschitz with $\mu=\prod_{k=1}^{L}\mu_{k}$ and
$\nu=\prod_{k=1}^{L}\nu_{k}$. The orthogonal matrix $P$ can be parameterized
via the Cayley transformation (6) or Householder transformation Singla et al.,
(2022). Since the distortion of ${\mathcal{O}}_{k}$ is 1, it can improve
network expressivity without increasing model distortion.

In some applications, e.g., normalising flows Dinh et al., (2015);
Papamakarios et al., (2021), we need to compute the inverse of
${\mathcal{G}}$, which can be done in a backward manner:

${\mathcal{G}}^{-1}(y)={\mathcal{O}}_{1}^{-1}\circ{\mathcal{F}}_{1}^{-1}\circ\cdots\circ{\mathcal{O}}_{K}^{-1}\circ{\mathcal{F}}_{K}^{-1}\circ{\mathcal{O}}_{K+1}^{-1}(y),$
(11)

where ${\mathcal{O}}_{k}$ has an explicit inverse
${\mathcal{O}}_{k}^{-1}(y)=P_{k}^{\top}(y-q_{k})$. Computing the inverse
${\mathcal{F}}_{k}^{-1}(y)$ requires an iterative solver, which will be
addressed in the next section.

#### Partially bi-Lipschitz networks.

A neural network
$\tilde{\mathcal{G}}:\mathbb{R}^{n}\times\mathbb{R}^{l}\rightarrow\mathbb{R}^{n}$
is said to be _partially bi-Lipschitz_ if for any fixed value of
$p\in\mathbb{R}^{l}$, the mapping $y=\tilde{\mathcal{G}}(x;p)$ is a
$(\mu,\nu)$-Lipschitz network from $x$ to $y$. We can construct such mappings
via $\tilde{\mathcal{G}}(x;p)={\mathcal{G}}_{h(p)}(x)$ where
${\mathcal{G}}_{\phi}$ is a$(\mu,\nu)$-Lipschitz network for any free
parameter $\phi\in\mathbb{R}^{N}$ and $h:p\rightarrow\phi$ is a new learnable
function. Since the dimension of $\phi$ is often very high, a practical
approach is to make $\phi$ partially depend on $p$. For instance, we can learn
$p$-dependent bias via an MLP with weight matrices independent of $p$.

##  4 Model inverse via operator splitting

In this section we give an efficient algorithm to compute
${\mathcal{F}}^{-1}(y)$ where ${\mathcal{F}}$ is a $\mu$-strongly monotone and
$\nu$-Lipschitz layer (8). First, we write its model inverse
${\mathcal{F}}^{-1}$ as

$\begin{split}\hat{z}&=\hat{\sigma}\left(\left(V-\frac{\gamma}{\mu}SS^{\top}\right)\hat{z}+b_{z}\right)\\\
x&=\frac{1}{\mu}(y-b_{y}-\sqrt{\gamma/2}S^{\top}\hat{z})\end{split}$ (12)

with $b_{z}=\sqrt{2\gamma}/\mu S(y-b_{y})+\hat{b}$. Both ${\mathcal{F}}$ and
${\mathcal{F}}^{-1}$ can be treated as special cases of _equilibrium networks_
Revay et al., (2020). The difference is that ${\mathcal{F}}$ has an explicit
formula due to the strictly lower-triangular $V$ while ${\mathcal{F}}^{-1}$ is
an implicit equation as $SS^{\top}$ is a full matrix. A natural question for
(12) is its well-posedness, i.e., for any $y\in\mathbb{R}^{n}$, does there
exists a unique $\hat{z}\in\mathbb{R}^{m}$ satisfying (12)?

######  Proposition 4.1.

${\mathcal{F}}^{-1}$ is well-posed if $V,S$ are given by (9).

A simple way to compute ${\mathcal{F}}^{-1}$ is the forward step method (FSM)
see, e.g., Ryu and Boyd, (2016)

$x^{k+1}=x^{k}-\alpha({\mathcal{F}}(x^{k})-y)$ (13)

which has a convergence rate of $1-\mu^{2}/\nu^{2}$ if $\alpha=\mu/\nu^{2}$.
In Winston and Kolter, (2020); Revay et al., (2020), certain classes of
equilibrium networks were solved via two-operator splitting problems. We
follow a similar strategy, but our structure requires a three-operator
splitting, see Proposition 4.2 with background in Appendix B. To state the
result, we first recall that for a monotone and 1-Lipschitz activation
function, there exists a proper convex function
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ satisfying
$\hat{\sigma}(\cdot)=\mathbf{prox}_{f}^{1}(\cdot)$ with

$\mathbf{prox}_{f}^{\alpha}(x)=\arg\min_{z\in\mathbb{R}^{n}}\;1/2\|x-z\|^{2}+\alpha
f(z).$

A list of $f$ for popular activations is given in Section B.1. In particular,
the ReLU activation satisfies
$\hat{\sigma}(\cdot)=\mathbf{prox}_{f}^{\alpha}(\cdot)$ for any $\alpha>0$,
which simplifies implementation of the proposed algorithm.

######  Proposition 4.2.

Finding a solution $\hat{z}\in\mathbb{R}^{m}$ to (12) is equivalent to finding
a zero to the three-operator splitting problem
$0\in{\mathcal{A}}(z)+{\mathcal{B}}(z)+{\mathcal{C}}(z)$ where
${\mathcal{A}},{\mathcal{B}},{\mathcal{C}}$ are monotone operators defined by

$\begin{split}{\mathcal{A}}(z)=(I-V)z-b_{z},\;{\mathcal{B}}(z)=\partial
f(z),\;{\mathcal{C}}(z)=\frac{\gamma}{\mu}SS^{\top}z\end{split}$

where $f$ satisfies $\hat{\sigma}(\cdot)=\mathbf{prox}_{f}^{1}(\cdot)$.

For three-operator problems, the Davis-Yin splitting algorithm (DYS) Davis and
Yin, (2017) can be applied, obtaining the following fixed-point iteration:

$\begin{split}z^{k+1/2}&=\mathbf{prox}_{f}^{\alpha}(u^{k})\\\
u^{k+1/2}&=2z^{k+1/2}-u^{k}\\\
z^{k+1}&=R_{{\mathcal{A}}}(u^{k+1/2}-\alpha{\mathcal{C}}(z^{k+1/2}))\\\
u^{k+1}&=u^{k}+z^{k+1}-z^{k+1/2}\end{split}$ (14)

where $R_{{\mathcal{A}}}(v)=((1+\alpha)I-\alpha V)^{-1}(v+\alpha b_{z})$.
Since $V$ is strictly lower triangular, we can solve $R_{{\mathcal{A}}}(v)$
using forward substitution. Furthermore, we can show that (14) is guaranteed
to converge with $\alpha\in\bigl{(}0,\frac{1}{\tau-1}\bigr{)}$, where $\tau$
is the model distortion.

##  5 Polyak-Łojasiewicz Networks

We call a network $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ a _Polyak-
Łojasiewicz (PL) network_ , or PLNet for short, if it satisfies the following
PL condition (Polyak,, 1963; Lojasiewicz,, 1963):

$\frac{1}{2}\|\nabla_{x}f(x)\|^{2}\geq m(f(x)-\min_{x}f(x)),\,\forall
x\in\mathbb{R}^{n},$ (15)

where $m>0$. The PL condition is significant in optimization since it is
weaker than convexity, but still implies that gradient methods converge to a
global minimum with a linear rate (Karimi et al.,, 2016), making PL networks a
promising candidate for learning a surrogate loss models.

######  Proposition 5.1.

If ${\mathcal{G}}$ is $\mu$-inverse Lipschitz, then

$f(x)=\frac{1}{2}\|{\mathcal{G}}(x)\|^{2}+c,\quad c\in\mathbb{R}$ (16)

is a PL network with $m=\mu^{2}$.

By construction, $f$ has a global minimum value of $c$, obtained uniquely at
the point $x^{\star}={\mathcal{G}}^{-1}(0).$

######  Remark 5.2.

We can relax the quadratic assumption:
$f(x)=h\bigl{(}{\mathcal{G}}(x)\bigr{)}$ is a PL network if
$h:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is strongly convex Karimi et al.,
(2016).

######  Remark 5.3.

For parametric optimization problem, one can learn a surrogate loss via
$f(x;p)=1/2\|\tilde{\mathcal{G}}(x;p)\|^{2}+c$ where $p\in\mathbb{R}^{m}$ is
the problem-specific parameter and $\tilde{\mathcal{G}}$ is a partially bi-
Lipschitz network.

######  Remark 5.4.

Any non-empty sub-level set ${\mathbb{L}}_{\alpha}=\\{x:f(x)<\alpha\\}$ is
homeomorphic to a unit ball, making PL networks suitable for learning Lyapunov
functions of stable dynamical systems Wilson Jr, (1967).

#### Computing global optimum of a PL network

If ${\mathcal{G}}$ is also $\nu$-Lipschitz (i.e., it is $(\mu,\nu)$-bi-
Lipschitz) then $f$ has a unique global optimum
$x^{\star}={\mathcal{G}}^{-1}(0)$ with ${\mathcal{G}}^{-1}$ given by (11).
This can be efficiently computed by analytical inversion orthogonal layers and
applying the DYS algorithm (14) to monotone and Lipschitz layers.

#### Limitations of gradient descent for finding global optimum.

An alternative way to compute the global optimum $x^{\star}$ is the standard
gradient descent (GD) method $x^{k+1}=x^{k}-\alpha\nabla_{x}f(x^{k})$. If
$\nabla_{x}f$ is $L$-Lipschitz, then the above GD solver with $\alpha=1/L$ has
a linear global convergence rate of $1-m/L$ with $m=\mu^{2}$ (Karimi et al.,,
2016). However, this method has two drawbacks. First, the gradient function
$\nabla_{x}f$ may not be globally Lipschitz, see Example 5.5. Secondly, even
if a global Lipschitz bound exists, it is still hard to obtain a good
estimation of $L$ for high-dimensional problems.

######  Example 5.5.

Consider a scalar function $f(x)=0.5g^{2}(x)$ with $g(x)=2x+\sin x$. Then,
$\partial f/\partial x=(2+\cos x)(2x+\sin x)$ is not globally Lipschitz due to
the term $2x\cos x$.

##  6 Experiments

Here we present experiments which explore the expressive quality of the
proposed models, regularisation via model distortion, and performance of the
DYS solution method.

###  6.1 Uncertainty quantification via neural Gaussian process

As shown in Liu et al., (2020), a model’s ability to quantify the distance of
a testing example from training data is a necessary condition for neural
networks to achieve high quality uncertainty estimation. The distance-
awareness can be achieved by imposing bi-Lipschitz smoothness for each
residual layer ${\mathcal{F}}(x)=x+{\mathcal{H}}(x)$ and replace the output
layer with a Gaussian process layer. The bi-Lipschitz regularization is done
by imposing Lipschitz bound of $0<c<1$ for ${\mathcal{H}}$ via spectral
normalization. The resulting model is called Spectral-normalized Neural
Gaussian Process (SNGP).

We replace the hidden space mapping with our bi-Lipschitz models and keep the
same output layer. We compare our $(\mu,\nu)$-Lipschitz network to a 3-layer
i-ResNet under the same bi-Lipschitz bound constraint, i.e., $\mu=(1-c)^{3}$
and $\nu=(1+c)^{3}$. Empirical results on the two-moon dataset are shown in
Figure 4. For the lower-distortion case (i.e., small $c=0.1$), SNGP fails to
completely separate the train and out-of-distribution (OOD) data due to its
loose Lipschitz bound. Our model can distinguish the OOD examples from
training dataset and the uncertainty surface is close to the SNGP with much
higher distortion ($c=0.9$). As the model distortion increases, our model can
have an uncertainty surface very close the dataset. The uncertainty surface of
SNGP does not change much form $c=0.1$ to $c=0.9$, see Figure 4 and additional
results in Section D.2.

Figure 4: The uncertainty plot of neural Gaussian process with different
models for the same bi-Lipschitz bounds. For the case with small distortion,
our model can still distinguish the train and OOD data, achieving similar
results of SNGP with large distortion.

###  6.2 Surrogate loss learning

We explore the PL Network’s performance with the Rosenbrock function
$r(x,y)=1/200(x-1)^{2}+0.5(y-x^{2})^{2}$ and its higher-dimensional
generalisations. The classical Rosenbrock function is a non-convex function
which has a unique global minimum point $(1,1)$. It has been widely used as a
test example for various of gradient-descent based optimization algorithms due
to the fact that its Hessian near the optimum is very poorly conditioned. We
also consider learning a surrogate loss for the sum of the Rosenbrock function
and a 2D sinewave function, chosen such that the overall loss still has a
unique global minimum at $(1,1)$ while having many local minima, see Section
D.1.

We learned models based on the form (16) where ${\mathcal{G}}$ is
parameterized by MLP, i-ResNet Behrmann et al., (2019), i-DenseNet Perugachi-
Diaz et al., (2021) and the proposed BiLipNet (10). We also trained the ICNN,
a scalar-output model which is convex w.r.t. inputs Amos et al., (2017).

Based on the results in Figure 7, we have the following observations. The
unconstrained MLP can achieve very small test errors. However, it has many
local minima near the valley $y=x^{2}$ since it is not inverse Lipschitz. This
phenomena is more easily visible for the Rosenbrock+Sine case but also occurs
in the plain Rosenbrock case. The ICNN model has a unique global minimum but
the fitting error is large as its sub-level sets are convex. For i-DenseNet,
the sub-level sets become mildly non-convex but their bi-Lipschitz bound is
quite conservative, so they do not capture the overall shape. In contrast, our
proposed BiLipNet is more flexible and captures the non-convex shape while
maintaining a unique global minimum. Additional results are in Section D.2.

#### Partially PL network.

We fit a parameterized Rosenbrock function:
$r(x,y;p)=1/200(x-a)^{2}+0.5(y-bx^{2})^{2}$ with $p=(a,b)$ by learning a
partially PL network (see Remark 5.3.) with $p$-dependent bias terms. The
results in Figure 8 indicate that the approach can be effective even if only
bias terms are modified by $p$, and not weights.

---
Figure 5: Surrogate loss learning for 20-dimensional Rosenbrock function.
Comparison of training and test error vs model distortion for PLNet with
different bi-Lipschitz models.

#### High-dimensional case.

We now turn to scalability of the approach to higher-dimensional problems and
analyse convergence of the DYS method for computing the global minimum. We
apply the approach to a $N$=20-dimensional version of the Rosenbrock function:

$r(x)=\frac{1}{N-1}\sum_{i=1}^{N-1}r(x_{i},x_{i+1})$ (17)

which has a global minimum of zero at $x=(1,1,...,1)$ but is non-convex and
has spurious local minima Kok and Sandrock, (2009).

We sample 10,000 training points uniformly over $[-2,2]^{20}$. Note that, in
contrast to the 2D example above, this is very sparse sampling of
20-dimensional space. We fit PL-networks of the form (16) with ${\mathcal{G}}$
constructed as the proposed BiLipNet as well as i-Resnet and i-DenseNet for
comparison.

A comparison of train and test error vs model distortion is shown in Figure 5.
It can be seen that our proposed BiLipNet model achieves far better fits than
iResNet Behrmann et al., (2021) and iDenseNet Perugachi-Diaz et al., (2021),
which can not achieve small training error for any value of the distortion
parameter. Furthermore, for our network, the distortion parameter appears to
act as an effective regularizer. Note that the best test error occurs after
training error drops to near zero ($\sim 10^{-8}$) but distortion is still
relatively small.

#### Solver comparison.

Given the surrogate loss function learned by BiLipNet, we now compare methods
to compute the location of its global minimum. In Figure 6 we compare the
proposed DYS solver to the FSM algorithm (13) and Adam Kingma and Ba, (2015)
applied directly to the surrogate loss. We take two values of the distortion
parameter: $\tau=5$ (optimal) and $\tau=50$. In both cases, the proposed DYS
method converges significantly faster than the alternatives, and the results
illustrate an additional benefit of regularising via distortion, besides
improving the test error: the $\tau=5$ case converges significantly faster
than $\tau=50$.

At the computed point $x^{\star}={\mathcal{G}}^{-1}(0)$ for the $\tau=5$ case,
the true function (17) takes a value of $r(x^{\star})=0.041$ . This is more
than an order of magnitude better than the smallest value of $r(x)$ over the
training data, which ranged over $[0.475,6.532]$, indicating that some useful
structure has been learned.

Figure 6: Solver comparison for finding the global minimum of a PLNet. We try
a range of rates $[0.1,0.5,1.0,2.0,5.0]$ for ADAM and present the best result.
The proposed back solve method with DYS algorithm (14) converges much faster
than ADAM applied to $f$ or back solve method with FSM algorithm (13).  Figure
7: Learning a surrogate loss for the Rosenbrock and Rosenbrock+Sine functions,
which is non-convex and has many local minima. The first row contains the true
functions while the remaining rows show learned functions and errors for
various surrogate loss models.
---
Figure 8: Learning a parameterized Rosenbrock function via a partially-PL
network constructed with a partially-BiLipNet.

##  7 Conclusion

In this paper we have introduced a “feed-through” layer with certified bounds
for strong monotonicity and Lipschitzness. By composing with orthogonal layers
we obtain a bi-Lipschitz network structure (BiLipNet) which has much tighter
bounds than those based on spectral normalization.

Building on these models, we have introduced the PLNet which is a scalar-
output neural network satisfying the Polyak-Łojasiewicz condition. PLNet has
many features making is suitable for learning surrogate loss functions,
including as a unique global minimum and a fast algorithm to compute it.

## References

  * Ahn et al., (2022) Ahn, B., Kim, C., Hong, Y., and Kim, H. J. (2022).  Invertible monotone operators for normalizing flows.  Advances in Neural Information Processing Systems, 35:16836–16848.
  * Amos et al., (2017) Amos, B., Xu, L., and Kolter, J. Z. (2017).  Input convex neural networks.  In International Conference on Machine Learning (ICML), pages 146–155. PMLR.
  * Araujo et al., (2023) Araujo, A., Havens, A. J., Delattre, B., Allauzen, A., and Hu, B. (2023).  A unified algebraic perspective on lipschitz neural networks.  In The Eleventh International Conference on Learning Representations (ICLR).
  * Ardizzone et al., (2018) Ardizzone, L., Kruse, J., Rother, C., and Köthe, U. (2018).  Analyzing inverse problems with invertible neural networks.  In International Conference on Learning Representations (ICLR).
  * Arjovsky et al., (2017) Arjovsky, M., Chintala, S., and Bottou, L. (2017).  Wasserstein generative adversarial networks.  In International conference on machine learning (ICML), pages 214–223. PMLR.
  * Arora and Doshi, (2021) Arora, S. and Doshi, P. (2021).  A survey of inverse reinforcement learning: Challenges, methods and progress.  Artificial Intelligence, 297:103500.
  * Bauer and Mnih, (2019) Bauer, M. and Mnih, A. (2019).  Resampled priors for variational autoencoders.  In The 22nd International Conference on Artificial Intelligence and Statistics, pages 66–75. PMLR.
  * Behrmann et al., (2019) Behrmann, J., Grathwohl, W., Chen, R. T., Duvenaud, D., and Jacobsen, J.-H. (2019).  Invertible residual networks.  In International conference on machine learning (ICML), pages 573–582. PMLR.
  * Behrmann et al., (2021) Behrmann, J., Vicol, P., Wang, K.-C., Grosse, R., and Jacobsen, J.-H. (2021).  Understanding and mitigating exploding inverses in invertible neural networks.  In International Conference on Artificial Intelligence and Statistics, pages 1792–1800. PMLR.
  * Chen et al., (2019) Chen, R. T., Behrmann, J., Duvenaud, D. K., and Jacobsen, J.-H. (2019).  Residual flows for invertible generative modeling.  Advances in Neural Information Processing Systems, 32.
  * Coleman et al., (2017) Coleman, C., Narayanan, D., Kang, D., Zhao, T., Zhang, J., Nardi, L., Bailis, P., Olukotun, K., Ré, C., and Zaharia, M. (2017).  Dawnbench: An end-to-end deep learning benchmark and competition.  Training, 100(101):102.
  * Cozad et al., (2014) Cozad, A., Sahinidis, N. V., and Miller, D. C. (2014).  Learning surrogate models for simulation-based optimization.  AIChE Journal, 60(6):2211–2227.
  * Davis and Yin, (2017) Davis, D. and Yin, W. (2017).  A three-operator splitting scheme and its optimization applications.  Set-valued and variational analysis, 25:829–858.
  * Davis, (2006) Davis, T. A. (2006).  Direct methods for sparse linear systems.  SIAM.
  * De Cao et al., (2020) De Cao, N., Aziz, W., and Titov, I. (2020).  Block neural autoregressive flow.  In Uncertainty in artificial intelligence, pages 1263–1273. PMLR.
  * Dinh et al., (2015) Dinh, L., Krueger, D., and Bengio, Y. (2015).  Nice: Non-linear independent components estimation.  ICLR Workshop Track.
  * Dinh et al., (2017) Dinh, L., Sohl-Dickstein, J., and Bengio, S. (2017).  Density estimation using real nvp.  In International Conference on Learning Representations (ICLR).
  * Fazlyab et al., (2019) Fazlyab, M., Robey, A., Hassani, H., Morari, M., and Pappas, G. (2019).  Efficient and accurate estimation of lipschitz constants for deep neural networks.  In Advances in Neural Information Processing Systems, pages 11427–11438.
  * Grathwohl et al., (2019) Grathwohl, W., Chen, R. T., Bettencourt, J., and Duvenaud, D. (2019).  Scalable reversible generative models with free-form continuous dynamics.  In International Conference on Learning Representations (ICLR).
  * Gu et al., (2016) Gu, S., Lillicrap, T., Sutskever, I., and Levine, S. (2016).  Continuous deep q-learning with model-based acceleration.  In International conference on machine learning (ICML), pages 2829–2838. PMLR.
  * Gulrajani et al., (2017) Gulrajani, I., Ahmed, F., Arjovsky, M., Dumoulin, V., and Courville, A. C. (2017).  Improved training of wasserstein gans.  Advances in neural information processing systems, 30.
  * Havens et al., (2023) Havens, A. J., Araujo, A., Garg, S., Khorrami, F., and Hu, B. (2023).  Exploiting connections between lipschitz structures for certifiably robust deep equilibrium models.  In Thirty-seventh Conference on Neural Information Processing Systems.
  * Helfrich et al., (2018) Helfrich, K., Willmott, D., and Ye, Q. (2018).  Orthogonal recurrent neural networks with scaled cayley transform.  In International Conference on Machine Learning (ICML), pages 1969–1978. PMLR.
  * Ho et al., (2019) Ho, J., Chen, X., Srinivas, A., Duan, Y., and Abbeel, P. (2019).  Flow++: Improving flow-based generative models with variational dequantization and architecture design.  In International Conference on Machine Learning (ICML), pages 2722–2730. PMLR.
  * Huang et al., (2018) Huang, C.-W., Krueger, D., Lacoste, A., and Courville, A. (2018).  Neural autoregressive flows.  In International Conference on Machine Learning (ICML), pages 2078–2087. PMLR.
  * Karimi et al., (2016) Karimi, H., Nutini, J., and Schmidt, M. (2016).  Linear convergence of gradient and proximal-gradient methods under the polyak-łojasiewicz condition.  In Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2016, Riva del Garda, Italy, September 19-23, 2016, Proceedings, Part I 16, pages 795–811. Springer.
  * Kingma and Ba, (2015) Kingma, D. P. and Ba, J. (2015).  Adam: A method for stochastic optimization.  In International Conference on Learning Representations (ICLR).
  * Kingma and Dhariwal, (2018) Kingma, D. P. and Dhariwal, P. (2018).  Glow: Generative flow with invertible 1x1 convolutions.  Advances in neural information processing systems, 31.
  * Kobyzev et al., (2020) Kobyzev, I., Prince, S. J., and Brubaker, M. A. (2020).  Normalizing flows: An introduction and review of current methods.  IEEE transactions on pattern analysis and machine intelligence, 43(11):3964–3979.
  * Kok and Sandrock, (2009) Kok, S. and Sandrock, C. (2009).  Locating and Characterizing the Stationary Points of the Extended Rosenbrock Function.  Evolutionary Computation, 17(3):437–453.
  * Li et al., (2019) Li, J., Fang, C., and Lin, Z. (2019).  Lifted proximal operator machines.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 4181–4188.
  * Li et al., (2020) Li, J., Li, F., and Todorovic, S. (2020).  Efficient riemannian optimization on the stiefel manifold via the cayley transform.  In International Conference on Learning Representations (ICLR).
  * Liang et al., (2023) Liang, E., Chen, M., and Low, S. (2023).  Low complexity homeomorphic projection to ensure neural-network solution feasibility for optimization over (non-) convex set.
  * Liu et al., (2020) Liu, J., Lin, Z., Padhy, S., Tran, D., Bedrax Weiss, T., and Lakshminarayanan, B. (2020).  Simple and principled uncertainty estimation with deterministic deep learning via distance awareness.  Advances in Neural Information Processing Systems, 33:7498–7512.
  * Lojasiewicz, (1963) Lojasiewicz, S. (1963).  A topological property of real analytic subsets.  Coll. du CNRS, Les équations aux dérivées partielles, 117(87-89):2.
  * Louizos and Welling, (2017) Louizos, C. and Welling, M. (2017).  Multiplicative normalizing flows for variational bayesian neural networks.  In International Conference on Machine Learning (ICML), pages 2218–2227. PMLR.
  * Lu et al., (2021) Lu, C., Chen, J., Li, C., Wang, Q., and Zhu, J. (2021).  Implicit normalizing flows.  In International Conference on Learning Representations (ICLR).
  * Megretski and Rantzer, (1997) Megretski, A. and Rantzer, A. (1997).  System analysis via integral quadratic constraints.  IEEE Transactions on Automatic Control, 42(6):819–830.
  * Meunier et al., (2022) Meunier, L., Delattre, B. J., Araujo, A., and Allauzen, A. (2022).  A dynamical system perspective for lipschitz neural networks.  In International Conference on Machine Learning (ICML), pages 15484–15500. PMLR.
  * Misener and Biegler, (2023) Misener, R. and Biegler, L. (2023).  Formulating data-driven surrogate models for process optimization.  Computers & Chemical Engineering, 179:108411.
  * Miyato et al., (2018) Miyato, T., Kataoka, T., Koyama, M., and Yoshida, Y. (2018).  Spectral normalization for generative adversarial networks.  In International Conference on Learning Representations (ICLR).
  * Papamakarios et al., (2021) Papamakarios, G., Nalisnick, E., Rezende, D. J., Mohamed, S., and Lakshminarayanan, B. (2021).  Normalizing flows for probabilistic modeling and inference.  The Journal of Machine Learning Research, 22(1):2617–2680.
  * Pauli et al., (2024) Pauli, P., Havens, A., Araujo, A., Garg, S., Khorrami, F., Allgöwer, F., and Hu, B. (2024).  Novel quadratic constraints for extending lipsdp beyond slope-restricted activations.  In International Conference on Learning Representations (ICLR).
  * Pauli et al., (2021) Pauli, P., Koch, A., Berberich, J., Kohler, P., and Allgöwer, F. (2021).  Training robust neural networks using lipschitz bounds.  IEEE Control Systems Letters, 6:121–126.
  * Perugachi-Diaz et al., (2021) Perugachi-Diaz, Y., Tomczak, J., and Bhulai, S. (2021).  Invertible densenets with concatenated lipswish.  Advances in Neural Information Processing Systems, 34:17246–17257.
  * Polyak, (1963) Polyak, B. (1963).  Gradient methods for minimizing functionals (in russian).  USSR Computational Mathematics and Mathematical Physics, 3(4):643–653.
  * Prach and Lampert, (2022) Prach, B. and Lampert, C. H. (2022).  Almost-orthogonal layers for efficient general-purpose lipschitz networks.  In European Conference on Computer Vision, pages 350–365. Springer.
  * Rantzer, (1996) Rantzer, A. (1996).  On the kalman—yakubovich—popov lemma.  Systems & control letters, 28(1):7–10.
  * Revay et al., (2020) Revay, M., Wang, R., and Manchester, I. R. (2020).  Lipschitz bounded equilibrium networks.  arXiv preprint arXiv:2010.01732.
  * Revay et al., (2023) Revay, M., Wang, R., and Manchester, I. R. (2023).  Recurrent equilibrium networks: Flexible dynamic models with guaranteed stability and robustness.  IEEE Transactions on Automatic Control.
  * Russo and Proutiere, (2021) Russo, A. and Proutiere, A. (2021).  Towards optimal attacks on reinforcement learning policies.  In 2021 American Control Conference (ACC), pages 4561–4567. IEEE.
  * Ryu and Boyd, (2016) Ryu, E. K. and Boyd, S. (2016).  Primer on monotone operator methods.  Appl. comput. math, 15(1):3–43.
  * Ryu et al., (2019) Ryu, M., Chow, Y., Anderson, R., Tjandraatmadja, C., and Boutilier, C. (2019).  Caql: Continuous action q-learning.  In International Conference on Learning Representations (ICLR).
  * Singla and Feizi, (2021) Singla, S. and Feizi, S. (2021).  Skew orthogonal convolutions.  In International Conference on Machine Learning (ICML), pages 9756–9766. PMLR.
  * Singla et al., (2022) Singla, S., Singla, S., and Feizi, S. (2022).  Improved deterministic l2 robustness on cifar-10 and cifar-100.  In International Conference on Learning Representations (ICLR).
  * Trockman and Kolter, (2021) Trockman, A. and Kolter, J. Z. (2021).  Orthogonalizing convolutional layers with the cayley transform.  In International Conference on Learning Representations (ICLR).
  * Tsuzuku et al., (2018) Tsuzuku, Y., Sato, I., and Sugiyama, M. (2018).  Lipschitz-margin training: Scalable certification of perturbation invariance for deep neural networks.  In Advances in neural information processing systems, pages 6541–6550.
  * Wang and Manchester, (2023) Wang, R. and Manchester, I. (2023).  Direct parameterization of lipschitz-bounded deep networks.  In International Conference on Machine Learning (ICML), pages 36093–36110. PMLR.
  * Ward et al., (2019) Ward, P. N., Smofsky, A., and Bose, A. J. (2019).  Improving exploration in soft-actor-critic with normalizing flows policies.  ICML Workshop on Invertible Neural Networks and Normalizing Flows,.
  * Wilson Jr, (1967) Wilson Jr, F. (1967).  The structure of the level surfaces of a lyapunov function.
  * Winston and Kolter, (2020) Winston, E. and Kolter, J. Z. (2020).  Monotone operator equilibrium networks.  Advances in neural information processing systems, 33:10718–10728.
  * Yeh, (2006) Yeh, J. (2006).  Real analysis: theory of measure and integration second edition.  World Scientific Publishing Company.

##  Appendix A Model Parameterization

Model parameterization is a mapping ${\mathcal{M}}:\phi\rightarrow\theta$
where $\phi\in\mathbb{R}^{N}$ is a free learnable parameter while $\theta$
includes the model weights $U\in\mathbb{R}^{m\times n},W\in\mathbb{R}^{m\times
m},Y\in\mathbb{R}^{n\times m}$ and IQC multiplier
$\Lambda\in{\mathbb{D}}_{+}^{m}$ with $m,n$ as the dimensions of hidden unit
and input, respectively. The aim of this section is to construct a
parameterization such that the large-scale SDP constraint (5) holds, i.e.,
$Y=U^{\top}\Lambda$ and

$\begin{split}H=2\Lambda-W^{\top}\Lambda-\Lambda
W=\begin{bmatrix}2\Lambda_{1}&-{W}_{2}^{\top}\Lambda_{2}\\\
-\Lambda_{2}W_{2}&2\Lambda_{2}&-{W}_{3}^{\top}\Lambda_{3}\\\
&\ddots&\ddots&\ddots\\\
&&-\Lambda_{L-1}W_{L-1}&2\Lambda_{L-1}&-W_{L}^{\top}\Lambda_{L}\\\
&&&-\Lambda_{L}{W}_{L}&2\Lambda_{L}\end{bmatrix}\geq\frac{2}{\gamma}Y^{\top}Y\end{split}$
(18)

Since $H\succeq 0$ has band structure, it can be represented by $H=XX^{\top}$
Davis, (2006). Moreover, from Lemma 3 of Rantzer, (1996) we have that any
$U,Y$ satisfying $Y=U^{\top}\Lambda$ and $XX^{\top}\succeq 2/\gamma Y^{\top}Y$
can be represent by

$U=\sqrt{\gamma/2}\Lambda^{-1}XQ,\quad Y=\sqrt{\gamma/2}Q^{\top}X^{\top}$ (19)

where $Q\in\mathbb{R}^{m\times n}$ with $QQ^{\top}\preceq I$. The remaining
task is to find $X$ such that $H=XX^{\top}$ has the same sparse structure as
(18), which was solved by Wang and Manchester, (2023). For self-contained
purpose, we provide detail construction as follows. First, we further
parameterize $X=\Psi P$, where $\Psi=\mathrm{diag}(\Psi_{1},\ldots,\Psi_{L})$
with $\Psi_{k}\in{\mathbb{D}}_{+}^{m_{k}}$ and

$\begin{split}P=\begin{bmatrix}A_{1}\\\ -B_{2}&A_{2}\\\ &\ddots&\ddots\\\
&&-B_{L}&A_{L}\end{bmatrix}.\end{split}$

By comparing $H=\Psi PP^{\top}\Psi$ with (18) we have

$\displaystyle
H_{kk}=\Psi_{k}(B_{k}B_{k}^{\top}+A_{k}A_{k}^{\top})\Psi_{k}=2\Lambda_{k},\quad
H_{k-1,k}=-\Psi_{k}B_{k}A_{k-1}^{\top}=-\Lambda_{k}W_{k},$

which further leads to

$\displaystyle\Psi_{k}^{2}=2\Lambda_{k},\quad
B_{k}B_{k}^{\top}+A_{k}A_{k}^{\top}=I,\quad
W_{k}=2\Psi_{k}^{-1}B_{k}A_{k-1}^{\top}\Psi_{k-1}\quad k=1,\ldots,L,$ (20)

with $B_{1}=0$. We have converted the large-scale SDP constraint (18) into
many simple and small-scale constraints such as

$\Psi_{k}^{2}=2\Lambda_{k},\quad R_{k}R_{k}^{\top}=I,\quad QQ^{\top}\preceq I$
(21)

with $R_{k}=\begin{bmatrix}B_{k}&A_{k}\end{bmatrix}$, which further can be
easily parameterized via the Cayley transformation (6), see Section 3.3. From
Proposition 3.5 we have that the transformation between (18) and (21) is
equivalent.

#### An equivalent model representation.

The model weight parameter $U,Y,W$ defined in (7) can be rewritten as
$U=\sqrt{2\gamma}\Psi^{-1}S$, $Y=\sqrt{\gamma/2}S^{\top}\Psi^{-1}$ and
$W=\Psi^{-1}W\Psi$ with

$S=\begin{bmatrix}S_{1}\\\ S_{2}\\\ \vdots\\\
S_{L}\end{bmatrix}=\begin{bmatrix}A_{1}Q_{1}\\\ A_{2}Q_{2}-B_{2}Q_{1}\\\
\vdots\\\ A_{L}Q_{L}-B_{L}Q_{L-1}\end{bmatrix},\quad V=\begin{bmatrix}0&\\\
V_{2}&0\\\ &\ddots&\ddots\\\ &&V_{L}&0\end{bmatrix}=\begin{bmatrix}0&\\\
2B_{2}A_{1}^{\top}&0\\\ &\ddots&\ddots\\\
&&2B_{L}A_{L-1}^{\top}&0\end{bmatrix}$ (22)

where $Q=\begin{bmatrix}Q_{1}^{\top}&\cdots&Q_{L}^{\top}\end{bmatrix}^{\top}$.
Then, the network (4) can be written as

$z=\sigma(\Psi^{-1}V\Psi z+\sqrt{2\gamma}\Psi^{-1}Sx+b),\quad y=\mu
x+\sqrt{\gamma/2}S^{\top}\Psi z.$ (23)

By introducing the new hidden state $\hat{z}=\Psi z$ and bias $\hat{b}=\Psi
b$, we obtain an equivalent form:

$\hat{z}=\hat{\sigma}\bigl{(}V\hat{z}+\sqrt{2\gamma}Sx+\hat{b}\bigr{)},\quad
y=\mu x+\sqrt{\gamma/2}S^{\top}\hat{z}+b_{y}.$ (24)

This representation is useful for model inverse via monotone operator
splitting, see the next section. We now give a lemma which will used later for
proving some propositions.

######  Lemma A.1.

For the matrices $V,S$ defined in (22) we have

$2I-V-V^{\top}\succeq 0,\quad 2I-SS^{\top}\succeq 0.$ (25)

###### Proof.

First, we have

$2I-(V+V^{\top})=2\begin{bmatrix}I&-A_{1}B_{2}^{\top}\\\
-B_{2}A_{1}^{\top}&I&-A_{2}B_{3}^{\top}\\\
&-B_{3}^{\top}A_{2}&\ddots&\ddots\\\ &&\ddots&\ddots\end{bmatrix}\succeq 0$
(26)

where the inequality is obtained by sequentially applying the fact
$A_{k}A_{k}^{\top}+B_{k}B_{k}^{\top}=I$ and Schur complement to the top
diagonal block. For the inequality on $S$, we have

$\begin{split}2I-SS^{\top}&=2I-PQQ^{\top}P^{\top}\succeq 2I-PP^{\top}\\\
&=2I-\begin{bmatrix}A_{1}\\\ -B_{2}&A_{2}\\\ &\ddots&\ddots\\\
&&-B_{L}&A_{L}\end{bmatrix}\begin{bmatrix}A_{1}\\\ -B_{2}&A_{2}\\\
&\ddots&\ddots\\\
&&-B_{L}&A_{L}\end{bmatrix}^{\top}=\begin{bmatrix}I&A_{1}B_{2}^{\top}\\\
B_{2}A_{1}^{\top}&I&A_{2}B_{3}^{\top}\\\ &B_{3}^{\top}A_{2}&\ddots&\ddots\\\
&&\ddots&\ddots\end{bmatrix}\succeq 0.\end{split}$

Similarly, the last inequality can be established by sequentially applying the
Schur complement to the top diagonal block. ∎

##  Appendix B Monotone Operator Splitting for Computing Model Inverse

Inspired by Winston and Kolter, (2020); Revay et al., (2020), we try to
compute $x={\mathcal{F}}^{-1}(y)$ via operator splitting method. We first
present some background of monotone operator theory based on the survey Ryu
and Boyd, (2016), and then reformulate the model inverse as a three-operator
splitting problem.

###  B.1 Monotone operator

An operator is a set-valued or single-valued map defined by a subset of the
space ${\mathcal{A}}\subseteq\mathbb{R}^{n}\times\mathbb{R}^{n}$; we use the
notation ${\mathcal{A}}(x)=\\{y\mid(x,y)\in{\mathcal{A}}\\}$. For example, the
affine operator is defined by ${\mathcal{L}}(x)=\\{(x,Wx+b)\mid
x\in\mathbb{R}^{n}\\}$. Another important example is the subdifferential
operator $\partial f=\\{(x,\partial f(x))\\}$ for a proper function
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}\cup\\{\infty\\}$ with $f(z)=\infty$ for
$z\notin\mathop{\bf dom}f$, where $\partial f(x)=\\{g\in\mathbb{R}^{n}\mid
f(y)\geq f(x)+\left\langle{y-x},{g}\right\rangle,\,\forall
y\in\mathbb{R}^{n}\\}$. An operator ${\mathcal{A}}$ has a Lipschitz bound of
$L$ if $\|u-v\|\leq L\|x-y\|$ for all $(x,u),(y,v)\in{\mathcal{A}}$. It is
non-expansive if $L=1$ and contractive if $L<1$. ${\mathcal{A}}$ is strongly
monotone with $m>0$ if

$\left\langle{u-v},{x-y}\right\rangle\geq
m\|x-y\|,\quad\forall(x,u),(y,v)\in{\mathcal{A}}.$ (27)

If the above inequality holds for $m=0$, we call ${\mathcal{A}}$ a monotone
operator. Similarly, ${\mathcal{A}}$ is said to be _inverse monotone_ with
$\rho$ if $\left\langle{u-v},{x-y}\right\rangle\geq\rho\|u-v\|,\
\forall(x,u),(y,v)\in{\mathcal{A}}$. An operator is called maximal monotone if
no other monotone operator strictly contains it. The linear operator
${\mathcal{L}}$ is $m$-strongly monotone if $W+W^{\top}\succeq 2mI$, and
$\rho$-inverse monotone if $W+W^{\top}\succeq 2\rho W^{\top}W$. A
subdifferential $\partial f$ is maximal monotone if and only if $f$ is a
convex closed proper (CCP) function. Here are some basic operations for
operators:

  * •

the operator sum
${\mathcal{A}}+{\mathcal{B}}=\\{(x,y+z)\mid(x,y)\in{\mathcal{A}},\,(x,z)\in{\mathcal{B}}\\}$;

  * •

the composition ${\mathcal{A}}{\mathcal{B}}=\\{(x,z)\mid\exists
y\;\mathrm{s.t.}\;(x,y)\in{\mathcal{A}},(y,z)\in{\mathcal{B}}\\}$ ;

  * •

the inverse operator
${\mathcal{A}}^{-1}=\\{(y,x)\mid(x,y)\in{\mathcal{A}}\\}$;

  * •

the resolvent operator $R_{{\mathcal{A}}}=(I+\alpha{\mathcal{A}})^{-1}$ with
$\alpha>0$;

  * •

the Cayley operator $C_{{\mathcal{A}}}=2R_{{\mathcal{A}}}-I$.

Note that the resolvent and Cayley operators are non-expansive for any maximal
monotone ${\mathcal{A}}$, and are contractive if ${\mathcal{A}}$ is strongly
monotone. For a linear operator ${\mathcal{L}}$ we have
$R_{{\mathcal{L}}}(x)=(I+\alpha W)^{-1}(x-\alpha b)$. For a subdifferential
operator $\partial f$, its resolvent is $R_{\partial
f}(x)=\mathbf{prox}_{f}^{\alpha}(x):=\operatorname*{arg\,min}_{z}1/2\|x-z\|+\alpha
f(z)$, which is also called the _proximal operator_.

#### Activation as proximal operator.

As shown in Li et al., (2019); Revay et al., (2020), many popular slope-
restricted scalar activation functions can also be treated as proximal
operators. To be specific, if $\sigma:\mathbb{R}\rightarrow\mathbb{R}$ is
slope-restricted in $[0,1]$, then there exists a convex proper function $f$
such that $\sigma(\cdot)=\mathbf{prox}_{f}^{1}(\cdot)$. For self-contained
purpose, we provide a list of common activations and their associated convex
proper functions in Table 1, which can also be found in Revay et al., (2020);
Li et al., (2019).

Table 1: A list of common activation functions and their associated convex proper $f(z)$ whose proximal operator is $\sigma(x)$ Revay et al., (2020). For $z\notin\mathop{\bf dom}f$, we have $f(z)=\infty$. In the case of Softplus activation, $\mathrm{Li}_{s}(z)$ is the polylogarithm function. Activation | $\sigma(x)$ | Convex $f(z)$ | $\mathop{\bf dom}f$
---|---|---|---
ReLu | $\max(x,0)$ | $0$ | $[0,\infty)$
LeakyReLu | $\max(x,0.01x)$ | $\frac{99}{2}\min(z,0)^{2}$ | $\mathbb{R}$
Tanh | $\tanh(x)$ | $\frac{1}{2}\left[\ln(1-z^{2})+z\ln\left(\frac{1+z}{1-z}\right)-z^{2}\right]$ | $(-1,1)$
Sigmoid | $1/(1+e^{-x})$ | $z\ln z+(1-z)\ln(1-z)-\frac{z^{2}}{2}$ | $(0,1)$
Arctan | $\arctan(x)$ | $-\ln(|\cos z|)-\frac{z^{2}}{2}$ | $(-1,1)$
Softplus | $\ln(1+e^{x})$ | $-\mathrm{Li}_{2}(e^{z})-i\pi z-z^{2}/2$ | $(0,\infty)$

###  B.2 Operator splitting

Many optimization problems (e.g. convex optimization) can be formulated as one
of finding a zero of an appropriate monotone operator ${\mathcal{F}}$, i.e.,
find $x\in\mathbb{R}^{n}$ such that $0\in{\mathcal{F}}(x)$. Note that $x$ is a
solution if and only if it is a fixed point $x={\mathcal{T}}(x)$ with
${\mathcal{T}}=I-\alpha{\mathcal{F}}$ for any nonzero $\alpha\in\mathbb{R}$.
The corresponding fixed point iteration is
$x^{k+1}={\mathcal{T}}(x^{k})=x^{k}-\alpha{\mathcal{F}}(x^{k})$. If
${\mathcal{F}}$ is $m$-strongly monotone and $L$-Lipschitz, then this
iteration converges by choosing $\alpha\in(0,2m/L^{2})$. The optimal
convergence rate is $1-(m/L)^{2}$, given by $\alpha=m/L^{2}$.

If ${\mathcal{F}}$ contains some non-smooth components, we then split
${\mathcal{F}}$ into two or three maximal operators:

two-operator splitting problem: $\displaystyle
0\in{\mathcal{A}}(x)+{\mathcal{B}}(x)$ (28) three-operator splitting problem:
$\displaystyle 0\in{\mathcal{A}}(x)+{\mathcal{B}}(x)+{\mathcal{C}}(x)$ (29)

where ${\mathcal{A}},{\mathcal{B}},$ and ${\mathcal{C}}$ are maximal monotone.
The main benefit of such splitting is that the resolvent or Cayley operators
for individual operator are easy to evaluate, which further leads to more
computationally efficient algorithms. For two-operator splitting problem, some
popular algorithms include

  * •

forward-backward splitting (FBS)
$x=R_{{\mathcal{B}}}(I-\alpha{\mathcal{A}})(x)$

  * •

forward-backward-forward splitting (FBFS)
$x=((I-\alpha{\mathcal{A}})R_{{\mathcal{B}}}(I-\alpha{\mathcal{A}})+\alpha{\mathcal{A}})(x)$

  * •

Peaceman-Rachford splitting (PRS)
$z=C_{{\mathcal{A}}}C_{{\mathcal{B}}}(z),\;x=R_{{\mathcal{B}}}(z)$

  * •

Douglas-Rachford splitting (DRS)
$z=(1/2I+1/2C_{{\mathcal{A}}}C_{{\mathcal{B}}})(z),\;x=R_{{\mathcal{B}}}(z)$

where the corresponding fixed-point iterations, the choices of hyper-parameter
$\alpha$ and convergence results can be found in Ryu and Boyd, (2016). For
three-operator splitting problem, the Davis-Yin splitting (DYS) Davis and Yin,
(2017) can be expressed by $z={\mathcal{T}}(z),\,x=R_{{\mathcal{B}}}(z)$ where
${\mathcal{T}}=C_{{\mathcal{A}}}(C_{{\mathcal{B}}}-\alpha{\mathcal{C}}R_{{\mathcal{B}}})-\alpha{\mathcal{C}}R_{{\mathcal{B}}}$.

###  B.3 Operator splitting perspective for ${\mathcal{F}}^{-1}$

As shown in the proof of Proposition 4.2, By applying the forward-backward
splitting with parameter $\alpha=1$, we can compute the solution $z$ via the
following iteration:

$\begin{split}z^{k+1}=R_{{\mathcal{B}}}(z^{k}-\widehat{{\mathcal{A}}}\bigl{(}z^{k}\bigr{)})=\hat{\sigma}\left(\left(V-\gamma/\mu
SS^{\top}\right)z^{k}+b_{z}\right).\end{split}$

It is worth pointing out that the above iteration may not converge for the
choice of $\alpha=1$. In practice we often use more stable and faster two-
operator splitting algorithms (e.g., PRS or DRS), see Winston and Kolter,
(2020); Revay et al., (2020). In this work, the motivation for further
decomposing the monotone operator $\widehat{\mathcal{A}}$ into two monotone
operators ${\mathcal{A}},{\mathcal{C}}$ is that $R_{{\mathcal{A}}}$ is a
large-scale linear equation with nice sparse structure while
$R_{\widehat{\mathcal{A}}}$ is dense due to the full weight matrix in
${\mathcal{C}}$.

#### Fixed-point iteration.

We now apply the DYS algorithm from Davis and Yin, (2017) to
$0\in{\mathcal{A}}(z)+{\mathcal{B}}(z)+{\mathcal{C}}(z)$, resulting in the
following fixed-point iteration:

$\begin{split}z^{k+1/2}&=R_{\mathcal{B}}(u^{k})=\mathbf{prox}_{f}^{\alpha}(u^{k})\\\
u^{k+1/2}&=2z^{k+1/2}-u^{k}\\\
z^{k+1}&=R_{\mathcal{A}}(u^{k+1/2}-\alpha{\mathcal{C}}(z^{k+1/2}))\\\
u^{k+1}&=u^{k}+z^{k+1}-z^{k+1/2}\end{split}$ (30)

where the third line is a large-scale sparse linear equation of the form

$\begin{bmatrix}(1+\alpha)I\\\ -\alpha V_{21}&(1+\alpha)I\\\ &\ddots&\ddots\\\
&&-\alpha V_{L,L-1}&(1+\alpha)I\end{bmatrix}\begin{bmatrix}z_{1}^{k+1}\\\
z_{2}^{k+1}\\\ \vdots\\\
z_{L}^{k+1}\end{bmatrix}=u^{k+1/2}+\alpha\left(b_{z}-\frac{\gamma}{\mu}SS^{\top}z^{k+1/2}\right).$

By introducing $v^{k+1/2}=b_{z}-\gamma/\mu SS^{\top}z^{k+1/2}$, we have

$z_{0}^{k+1}=0,\quad
z_{l}^{k+1}=\frac{\alpha}{1+\alpha}\left(V_{l,l-1}z_{l-1}^{k+1}+v_{l}^{k+1/2}\right)+\frac{1}{1+\alpha}u_{l}^{k+1/2},\quad
l=1,\ldots,L.$ (31)

#### Convergence range for the hyper-parameter $\alpha$.

From the previous paragraph, we know that (14) is equivalent to the FPI (30).
From Theorem 1.1 of Davis and Yin, (2017), we have that (30) converges for any
$\alpha\in(0,2\beta)$ with $\beta$ as the inverse-monotone bound of
${\mathcal{C}}$. From Lemma A.1 we have $2I\succeq S^{\top}S$ and

$\frac{2\gamma}{\mu}SS^{\top}\succeq\frac{\gamma}{\mu}S(S^{\top}S)S^{\top}=\frac{\mu}{\gamma}(\gamma/\mu
SS^{\top})^{2}=2\beta(\gamma/\mu SS^{\top})^{2}$

i.e., ${\mathcal{C}}(z)$ is inverse monotone with $\beta=\mu/(2\gamma)$.
Therefore, (14) converges for any $\alpha\in(0,\mu/\gamma)$. Since
$\gamma=\nu-\mu$ and $\tau=\nu/\mu$, we then obtain the convergence range in
term of model distortion $\tau$, i.e., $\alpha\in(0,1/(\tau-1))$. Larger
$\alpha$ often implies faster convergence rate, see Figure 9.

Figure 9: Forward step method (FSM) vs Davis-Yin splitting (DYS) applied to
randomly initialized $\mu$-monotone and $\nu$-Lipschitz layer with different
input and hidden unit dimensions. We obverse that DYS converges faster for
lager $\alpha$. If $\alpha$ is close to the bound $\mu/\gamma$ with
$\gamma=\nu-\mu$, it has a much faster rate than FSM.

##  Appendix C Proofs

###  C.1 Proof of Theorem 3.2

We consider the neural network ${\mathcal{H}}:x\rightarrow\tilde{y}$ defined
by

$v=Wz+Ux+b,\quad z=\sigma(v),\quad\tilde{y}=Yz+b_{y}.$ (32)

Since ${\mathcal{F}}(x)=\mu x+{\mathcal{H}}(x)$, then ${\mathcal{F}}$ is
$\mu$-strongly monotone and $\nu$-Lipschitz if ${\mathcal{H}}$ is monotone and
$\gamma$-Lipschitz with $\gamma=\nu-\mu$.

For any pair of solutions $s_{1}=(x_{1},v_{1},z_{1},\tilde{y}_{1})$ and
$s_{2}=(x_{2},v_{2},z_{2},\tilde{y}_{2})$, their difference $\Delta
s=s_{1}-s_{2}$ satisfies

$\Delta v=W\Delta z+U\Delta x,\quad\Delta z=J_{\sigma}(v_{1},v_{2})\Delta
v,\quad\Delta\tilde{y}=Y\Delta z$ (33)

where $J_{\sigma}$ is a diagonal matrix with $[J_{\sigma}]_{ii}\in[0,1]$ since
$\sigma$ is an elementwise activation with slope restricted in $[0,1]$. For
any $\Lambda\in\mathbb{D}_{+}^{m}$ we have

$\left\langle{\Delta v-\Delta z},{\Lambda\Delta z}\right\rangle=\Delta
v^{\top}(I-J_{\sigma})\Lambda J_{\sigma}\Delta_{v}\geq 0,\quad\forall\Delta
v\in\mathbb{R}^{m}.$ (34)

Based on (33), (34) and Condition (5) we have

$\begin{split}&\left\langle{\Delta
x},{\Delta\tilde{y}}\right\rangle-\left\langle{\Delta v-\Delta
z},{\Lambda\Delta z}\right\rangle\\\ =&\left\langle{\Delta x},{Y\Delta
z}\right\rangle-\left\langle{(W-I)\Delta z+U\Delta x},{\Lambda\Delta
z}\right\rangle\\\ =&\left\langle{\Delta x},{Y\Delta
z}\right\rangle-\left\langle{\Delta x},{U^{\top}\Lambda\Delta
z}\right\rangle+\left\langle{(I-W)\Delta z},{\Lambda\Delta z}\right\rangle\\\
=&\bm{1}/2\Delta z^{\top}\left(\Lambda(I-W)+(I-W^{\top})\Lambda\right)\Delta
z\geq\|Y\Delta z\|^{2}\geq 0,\end{split}$ (35)

which further implies $\left\langle{\Delta x},{\Delta
y}\right\rangle-\mu\|\Delta x\|^{2}\geq\left\langle{\Delta v-\Delta
z},{\Lambda\Delta z}\right\rangle\geq 0$. Thus, ${\mathcal{H}}$ is monotone.
We can use the similar technique to derive the Lipschitz bound of
${\mathcal{H}}$. Firstly we have

$\begin{split}&\gamma\|\Delta
x\|^{2}-\frac{1}{\gamma}\|\Delta\tilde{y}\|^{2}-2\left\langle{\Delta v-\Delta
z},{\Lambda\Delta z}\right\rangle\\\ =&\gamma\|\Delta
x\|^{2}-\frac{1}{\gamma}\|\Delta\tilde{y}\|^{2}+2\left\langle{(I-W)\Delta
z},{\Lambda\Delta z}\right\rangle-2\left\langle{U\Delta x},{\Lambda\Delta
z}\right\rangle\\\ =&\gamma\|\Delta x\|^{2}-2\left\langle{\Delta
x},{\Delta\tilde{y}}\right\rangle-\frac{1}{\gamma}\|\Delta\tilde{y}\|^{2}+\Delta
z^{\top}(2\Lambda-\Lambda W-W^{\top}\Lambda)\Delta z\\\ \geq&\gamma\|\Delta
x\|^{2}-2\left\langle{\Delta
x},{\Delta\tilde{y}}\right\rangle-\frac{1}{\gamma}\|\Delta\tilde{y}\|^{2}+\frac{2}{\gamma}\|Y\Delta
z\|^{2}=\left\|\sqrt{\gamma}\Delta
x-\frac{1}{\sqrt{\gamma}}\Delta\tilde{y}\right\|^{2}\geq 0.\end{split}$ (36)

Due to (34) we can further obtain $\gamma^{2}\|\Delta
x\|^{2}\geq\|\Delta\tilde{y}\|^{2}$, i.e., ${\mathcal{H}}$ is
$\gamma$-Lipschitz.

###  C.2 Proof of Proposition 3.5

Sufficient part: (7) $\Rightarrow$ (5). From (19) we have that
$Y=U^{\top}\Lambda$. We check the inequality part of (5) as follows:

$\begin{split}2\Lambda-W^{\top}\Lambda-\Lambda W=\Psi
PP^{\top}\Psi=XX^{\top}\succeq
XQQ^{\top}X^{\top}=\frac{2}{\gamma}Y^{\top}Y.\end{split}$

Necessary part: (5) $\Rightarrow$ (7). Since $H\succeq 0$ has band structure,
then it can be decomposed into $H=XX^{\top}$ where $X$ has the following block
lower triangular structure Davis, (2006):

$X=\begin{bmatrix}X_{11}\\\ X_{21}&X_{22}\\\ &\ddots&\ddots\\\
&&X_{L,L-1}&X_{LL}\end{bmatrix}.$

For this special case, a way to construct $X$ from $\Lambda,W$ and further
computation of the free parameters $d,F_{k}^{a},F_{k}^{b},F^{q},F^{\star}$ can
be found in Wang and Manchester, (2023). Finally, we need to show that
$XX^{\top}\succeq 2/\gamma Y^{\top}Y$ is equivalent to
$Y=\sqrt{\gamma/2}Q^{\top}X^{\top}$ for some $QQ^{\top}\preceq I$, which can
be directly followed by Lemma 3 of Rantzer, (1996).

###  C.3 Proof of Proposition 4.1

From Lemma A.1 we have

$2I-(V-\gamma/\mu SS^{\top})-(V-\gamma/\mu
SS^{\top})^{\top}=2I-V-V^{\top}+2\gamma/\mu SS^{\top}\succeq 2\gamma/\mu
SS^{\top}\succeq 0.$ (37)

Then, the equilibrium network (12) is well-posed by Theorem 1 of Revay et al.,
(2020).

###  C.4 Proof of Proposition 4.2

We first show that $0\in{\mathcal{A}}(z)+{\mathcal{B}}(z)+{\mathcal{C}}(z)$ is
a monotone operator splitting problem. It is obvious that
${\mathcal{B}},{\mathcal{C}}$ are maximal monotone operators. From Lemma A.1
we have $(I-V)+(I-V)^{\top}\succeq 0$, i.e. ${\mathcal{A}}$ is also monotone.
Then, we show that the above operator splitting problem shares the same set of
equilibrium points with the model inverse (12). First, we rewrite it into a
two-operator splitting problem
$0\in\widehat{{\mathcal{A}}}(z)+{\mathcal{B}}(z)$ where
$\widehat{{\mathcal{A}}}={\mathcal{A}}+{\mathcal{C}}$. By applying the
forward-backward splitting with parameter $\alpha=1$, we can compute the
solution $z$ via the following iteration:

$\begin{split}z^{k+1}=&R_{{\mathcal{B}}}(z^{k}-\widehat{{\mathcal{A}}}\bigl{(}z^{k}\bigr{)})=\mathbf{prox}_{f}^{1}\left(z^{k}-\left(I-V+\gamma/\mu
SS^{\top}\right)z^{k}+b_{z}\right)=\hat{\sigma}\left(\left(V-\gamma/\mu
SS^{\top}\right)z^{k}+b_{z}\right).\end{split}$

Thus, any solution $z^{\star}$ of the equilibrium network (12) is also an
equilibrium point of the above iteration.

###  C.5 Proof of Proposition 5.1

First, we have $\nabla f(x)=G^{\top}(x){\mathcal{G}}(x)$ where
$G(x)=\nabla{\mathcal{G}}(x)$ satisfies $\|G(x)\|\geq\mu$. Then, the PL
inequality holds for $f$ with $m=\mu^{2}$, i.e.,

$\frac{1}{2}\|\nabla
f(x)\|^{2}=\frac{1}{2}{\mathcal{G}}(x)^{\top}G(x)^{\top}G(x){\mathcal{G}}(x)\geq\frac{\mu^{2}}{2}\|{\mathcal{G}}(x)\|^{2}=\mu^{2}(f(x)-f^{\star}).$
(38)

##  Appendix D Experiments

###  D.1 Training details

For all experiments, we choose ReLU as our default activation and use Adam
Kingma and Ba, (2015) with one-cycle linear learning rate Coleman et al.,
(2017) except the NGP case which use a fixed rate. For the NGP case, we use
the cross entropy loss while the L2 loss is used for the rest of the examples.
For our proposed BiLipNet, we fixed the free parameter $F_{\star}=0$, which
implies that $Q^{\top}Q=I$ instead of $Q^{\top}Q\preceq I$. By enforcing all
singular values of $Q$ to be 1 (as suggested by Prach and Lampert, (2022),
this can improve the model training when there are multiple monotone and
Lipschitz layers. Dataset and model architectures are described as follows.

#### 1D Step function.

The target function is a step function

$f(x)=\begin{cases}2,&x>0\\\ -2,&x<0\end{cases}$

which is monotone and $0$-Lipschitz everywhere except the singularity point
$x=0$. We try to fit this curve with $(0.1,10)$-Lipschitz models. The optimal
fit is a linear piecewise continuous function with slope of 10 near $x=0$ and
slope of 0.1 near $x=\pm 2$. We take 1000 random samples from $[-2,2]$ for
training. Our model (BiLipNet) is an one-layer residual network
${\mathcal{F}}(x)=\mu x+{\mathcal{H}}(x)$ where ${\mathcal{H}}$ has 8 hidden
layers of width 32, giving the model 15.8K parameters. We compare to i-ResNet
Chen et al., (2019) and i-DenseNet Perugachi-Diaz et al., (2021), where the
nonlinear block ${\mathcal{H}}$ has 2 and 4 hidden layers, respectively. For
those two models, we test for depth from 2 to 8 with proper hidden width (so
that they has similar amount of parameters). And the empirical Lipschitz bound
is computed via finite difference over the test data. As shown in Figure 1,
our model achieves much tighter bounds than other models.

#### Neural Gaussian process.

We take 1000 two-moon data points as training data and 1000 Gaussian samples
with mean $(1.3,-1.8)$ and variance $(0.02,0.01)$ as OOD data. For all models,
we use fixed input weight to mapping the 2D input into 128D hidden space, then
perform hidden space transformation using bi-Lipschitz models, and finally add
a Gaussian process as the output layer. SNGP uses 3 residual layer
$x+{\mathcal{H}}(x)$ where the Lipschitz bound of ${\mathcal{H}}$ is $c<1$.
BiLipNet has one monotone and Lipschitz layer with two orthogonal layer, i.e.,
$K=1$ for (10). The nonlinear block ${\mathcal{H}}$ of our model has 6 hidden
layers with width of 32. Both models are chosen to have the same amount of
parameters, roughly $233K$.

#### 2D Rosenbrock function.

The true function is a Rosenbrock function defined by

$r(x,y)=\frac{1}{200}(x-1)^{2}+\frac{1}{2}\bigl{(}y-x^{2}\bigr{)}^{2}.$

Note that we use a scaling factor of $1/200$ for the classic Rosenbrock
function. The above function is non-convex but has one minimum at $(1,1)$. We
also consider the combination of the above Rosenbrock function with the
following 2D Sine function:

$s(x,y)=0.25(\sin(8(x-1)-\pi/2)+\sin(8(y-1)-\pi/2)+2).$

In this case $r(x,y)+s(x,y)$ still has a unique global minimum at $(1,1)$. But
there are many local minima. We take 5K random training samples from the
domain $[-2,-2]\times[-1,3]$. The proposed BiLipNet contains two monotone and
Lipschitz layers (i.e., $K=2$ for (10)). The nonlinear block ${\mathcal{H}}$
has 4 hidden layers of width 128. The model size is roughly 16K. The ICNN
model has 8 hidden layers with width of 180. The MLP has hidden units of
$[128,256,256,512]$. We trained i-ResNet and i-DenseNet with different depth
and width such that the total amount of parameters is comparable with
BiLipNet.

#### Parametric Rosenbrock function.

We consider the following parametric Rosenbrock function

$r(x,y;p)=\frac{1}{200}(x-a)^{2}+\frac{1}{2}(y-bx^{2})^{2},\quad
p=(a,b)\in[-1,1]^{2}.$

We take 10K random training data. The partially BiLipNet contains 3 orthogonal
layers, and 2 monotone and Lipschitz layers (the ${\mathcal{H}}$ block of each
layer has 4 hidden layer with width 128). The bias term of each orthogonal
layer is produced by an MLP with hidden units of $[64,128,2]$ while the bias
for those hidden units inside the ${\mathcal{H}}$ block is generated by an MLP
of $[64,128,256,512]$. The model’s bi-Lipschitz bound is chosen to be
$(0.04,16)$. The resulting model size is 604K.

#### ND Rosenbrock function.

We also consider the $N$-dimensional (with $N=20$) Rosenbrock function:

$R_{N}(x)=\frac{1}{N-1}\sum_{i=1}^{N-1}r(x_{i},x_{i+1})$

which is non-convex and has a unique global minimum at $(1,1,\ldots,1)$.
Besides it also has many local minima. We take 10K random samples over the
domain $[-2,2]^{20}$ and do training with batch size of 200. Note that the
data size is very small compared to the dimension. We then use 500K samples
for testing. BiLipNet has two monotone and Lipschitz layers (i.e., $K=2$ for
(10)) where each layer has a nonlinear block ${\mathcal{H}}$ with 8 hidden
layer of width 256 (model size $\sim$ 2.1M). For the i-ResNet/i-DenseNet, we
try different depths from 2 to 10 and observe that depth of 5 yields slightly
better results. The width of hidden layer is chosen so that it has a similar
amount of parameters as BiLipNet.

###  D.2 Extra results

Some extra results for uncertainty quantification via neural Gaussian process
and surrogate loss learning are shown in in Figure 10 and Figure 11,
respectively.

Figure 10: Uncertainty qualification via neural Gaussian process with
different bi-Lipschitz bound specifications on the hidden space mapping.
Figure 11: Additional results for Learning a surrogate loss for the Rosenbrock
and Rosenbrock + Sine functions. The first row contains the true functions
while the remaining rows show learned functions and errors for various
surrogate loss models. Our model has the flexibility of capturing the non-
convex sub-level sets, but can also fit smoothed representations by reducing
the distortion parameter.