13. Lipschitz-by-Design Networks & Direct Parameterizations

Orthogonal and Cayley layers, SLL, Sandwich layers, Pauli's Cayley–Gramian CNNs and LipKernel, RENs, and the certified-robustness state of the art

Book contents · Apply this chapter to a real decision · Glossary

A scalar model of “by construction”

Suppose a scalar weight must obey $|w|\le1$. The free parameterization $w=\tanh\theta$ meets the constraint for every real $\theta$: training can move $\theta$ freely. It reaches only $(-1,1)$, so validity does not imply completeness at the boundary. In contrast, adding a penalty $\lambda w^2$ to a loss lets training trade a larger weight against a smaller loss; a finite penalty does not enforce $|w|\le1$.

The matrix constructions replace this interval fact by identities such as $Q^\top Q=I$ or $H=PP^\top\succeq0$. Before reading their formulas, check the dimensions, the domain of each inverse, and the set for which completeness is claimed. Matrix multiplication and transposes and Gram matrices and PSD factors supply the needed algebra.

Module 12 answered the analysis question: given the weights, LipSDP returns a certified Lipschitz bound, and constrained training (ADMM, barriers) can push it down. This module turns the question around and builds the certificate into the architecture: every weight is a function of free parameters (smooth in the Cayley-based constructions), chosen so that the Lipschitz certificate holds for every value of the free parameters in its stated domain. For a direct parameterization with domain $\mathbb R^N$, training is then unconstrained SGD (see Primer B), and every iterate, including the deployed one, is certified in exact arithmetic. The Cayley–Gramian construction below additionally requires nonsingular gain factors; approximate orthogonal updates need their numerical error included.

The algebra is small: a Schur complement turns a Lipschitz LMI into a norm bound, and a norm bound is met by construction with a matrix with orthonormal columns, which the Cayley transform produces from arbitrary matrices. The control-theoretic twist, due to Patricia Pauli and co-authors, is to treat convolutions as dynamical systems (with Dennis Gramlich and Frank Allgöwer), so that the storage matrix of their dissipativity certificate can be chosen as the inverse of a controllability Gramian (see Primer D) in closed form (with Ruigang Wang, Ian Manchester and Frank Allgöwer). The same ideas give recurrent models that are contracting by construction (RENs, R2DN), which Module 14 uses as certified controllers.

1. Constrain or Parameterize?

The training problem we would like to solve is

where $\mathcal L$ is the training loss, $\theta$ collects the trainable parameters and "s.t." means "subject to": a constrained optimization problem whose feasible points are the $\theta$ satisfying the constraint (see Primer B). Computing $\operatorname{Lip}(f_\theta)$ exactly is NP-hard already for two-layer ReLU networks (Virmaux & Scaman, NeurIPS 2018), so the constraint is always replaced by a certificate set $\Theta_{\rm cert}(\gamma)\subseteq\{\theta:\operatorname{Lip}(f_\theta)\le\gamma\}$, e.g. "there is a diagonal $\Lambda\succ0$ satisfying the LipSDP LMI" (Fazlyab et al., NeurIPS 2019; Module 12). There are three ways to live with that set.

1. Constrain: optimize over $\Theta_{\rm cert}$ with ADMM (Pauli et al., L-CSS 2022), log-det barriers (Pauli et al., CDC 2022) or projections (Module 12). The set is not jointly convex in $(W,\Lambda)$; the substitution $\hat W=\Lambda W$ makes it an LMI, but every step still touches a semidefinite constraint. Normalized training times on the CDC 2022 toy example: 1 (nominal), 7.35 (projected), 78.72 (ADMM), 1.15–1.30 (barrier); Wang & Manchester report that barriers and projections become the bottleneck around $10^3$ neurons.
   Derivation — why $\hat W=\Lambda W$ gives an LMI

   In the one-hidden-layer LMI of Section 3, $H=\begin{bmatrix}\gamma I&-W_0^{\top}\Lambda&0\\-\Lambda W_0&2\Lambda&-W_1^{\top}\\0&-W_1&\gamma I\end{bmatrix}\succeq0$, the only products of unknowns are $\Lambda W_0$ and $W_0^{\top}\Lambda$. With the new variable $\hat W_0=\Lambda W_0$ the off-diagonal blocks become $-\hat W_0^{\top},-\hat W_0,-W_1^{\top},-W_1$, so $H$ is affine in $(\hat W_0,\Lambda,W_1,\gamma)$: an LMI. The network is recovered as $W_0=\Lambda^{-1}\hat W_0$ ($\Lambda$ is positive diagonal, hence invertible). Only the constraint becomes convex; the training loss, written in the new variables, is still nonconvex.

2. Regularize: penalize a differentiable bound, e.g. the loop-transformation bound LipLT with a margin regularizer (Fazlyab, Entesari, Roy & Chellappa, NeurIPS 2023). Cheap, but nothing forces the certified bound below a prescribed $\gamma$.
   Going deeper — what LipLT computes, and why a penalty is not a constraint

   A bound routine returns a differentiable $U(\theta)\ge\operatorname{Lip}(f_\theta)$, and training minimizes $\mathcal L(f_\theta)+\lambda U(\theta)$ with $\lambda\ge0$ (or a margin loss built from $U$). For one hidden layer, $f(x)=W_1\varphi(W_0x+b_0)+b_1$ with slopes in $[0,1]$, the loop transformation splits $\varphi(t)=\tfrac t2+r(t)$: $r$ has slopes in $[-\tfrac12,\tfrac12]$, so it is $\tfrac12$-Lipschitz, and $f(x)=\tfrac12W_1W_0x+W_1r(W_0x+b_0)+\text{const}$. Bounding the two terms separately gives

   $$\operatorname{Lip}(f)\le U=\tfrac12\|W_1W_0\|_2+\tfrac12\|W_1\|_2\|W_0\|_2\ \le\ \|W_1\|_2\|W_0\|_2 ,$$

   never worse than the product bound because the first term sees the combined direction of $W_1W_0$; deeper networks repeat the transformation layer by layer (Fazlyab et al., NeurIPS 2023). A finite $\lambda$ lets training trade a larger $U$ for a smaller loss, so nothing enforces $U(\theta)\le\gamma$.

3. Parameterize: choose a smooth map whose image lies inside $\Theta_{\rm cert}(\gamma)$ and train over its unconstrained argument. This module.

Layer-wise versus network-wise certificates

The cheapest certificate multiplies layer bounds, $\operatorname{Lip}(f)\le\prod_k\|W_k\|$ for 1-Lipschitz activations (Module 12); making every layer 1-Lipschitz (Section 2) is the corresponding design principle. It ignores that the directions amplified by one layer need not be those amplified by the next. Network-level certificates such as LipSDP couple the layers through multipliers, and their parameterizations produce networks whose individual layers have norms well above 1 while the network is still $\gamma$-Lipschitz (Wang & Manchester, Remark C.2; see the explorer). LipKernel's Fig. 1 shows it geometrically: propagating ellipsoids through a two-layer network certifies bound 1 where propagating balls certifies almost 2.

Expressivity: why gradient norm has to be preserved

In words: backpropagation multiplies the gradient by $W_k^{\top}$ (norm at most 1) and by diagonal activation Jacobians with entries in $[0,1]$; to arrive with norm exactly 1, every unit that matters must sit on a slope-1 piece everywhere, so $f$ is affine. The norm bound makes each factor non-expansive; "elementwise and monotone" restricts the Jacobian to a diagonal with entries in $[0,1]$ (a sorting activation escapes because its Jacobian is a permutation). Hence a ReLU network with $\|W_k\|\le1$ cannot represent $|x|$ and can be too restrictive for Wasserstein critics, which optimize over all scalar 1-Lipschitz functions (Anil et al., ICML 2019). The remedy: weights may be taken orthonormal without changing the function (their Thm 2, under the unit-gradient hypothesis; stated below), and activations must be gradient-norm preserving: GroupSort sorts groups of pre-activations, and its group-size-2 case MaxMin maps $(a,b)\mapsto(\max(a,b),\min(a,b))$. With both, scalar-output networks approximate every 1-Lipschitz function (their Thm 3, below).

The table orders the approaches of this module by what holds for every parameter value; "layer" certificates make each layer 1-Lipschitz, "network" certificates couple layers through multipliers or gain matrices.

2. Spectral Normalization, Parseval, Cayley, SOC, AOL

The layer-wise route makes every affine layer 1-Lipschitz and uses $\operatorname{Lip}(g\circ h)\le\operatorname{Lip}(g)\operatorname{Lip}(h)$. A linear layer $x\mapsto Wx$ is 1-Lipschitz iff $\|W\|\le1$. It is gradient-norm preserving, i.e. backpropagation $g\mapsto W^{\top}g$ keeps $\|g\|$ for every $g$ (Li et al., NeurIPS 2019), iff $WW^{\top}=I$ (orthonormal rows); the forward map keeps every $\|x\|$ iff $W^{\top}W=I$ (orthonormal columns), and a square orthogonal $W$ does both. A tall $W$ with orthonormal columns has all singular values 1 but is not gradient-norm preserving: $W=(1,0)^{\top}$ sends $g=(0,1)^{\top}$ to $0$. The methods differ in how exactly and cheaply they enforce this, and whether they work for convolutions.

Spectral normalization and Parseval networks

Parseval networks keep a weight $W\in\mathbb R^{d_{\rm out}\times d_{\rm in}}$, $d_{\rm out}\le d_{\rm in}$, close to a Parseval tight frame, $WW^{\top}=I$ (orthonormal rows), with the Frobenius-norm regularizer $R_\beta(W)=\tfrac\beta4\|WW^{\top}-I\|_F^2$, whose gradient is $\beta(WW^{\top}-I)W$ (see Primer B), and take one unit gradient step on it after every main update: $W\leftarrow(1+\beta)W-\beta WW^{\top}W$ (Cisse et al., ICML 2017; the paper prints the penalty as $\tfrac\beta2\|W^{\top}W-I\|_2^2$, which does not have this gradient). In the SVD $W=U\Sigma V^{\top}$ this maps each singular value $s\mapsto(1+\beta)s-\beta s^3$, which has an attracting fixed point at $s=1$ for $0\lt\beta\lt1$ (see Primer D). The result is only approximately orthogonal, so there is no certificate without a post-hoc norm computation.

The Cayley transform

Two dimensions. For $A=\begin{bmatrix}0&a\\-a&0\end{bmatrix}$, $Q=\frac{1}{1+a^2}\begin{bmatrix}1-a^2&-2a\\2a&1-a^2\end{bmatrix}$, a rotation by $\theta=2\arctan a$ (put $a=\tan(\theta/2)$ in the double-angle formulas). As $a\to\pm\infty$ the angle tends to $\pm\pi$ without reaching it: $Q=-I$ is missing, and so are all reflections. Trockman & Kolter note that a fixed diagonal $\pm1$ factor recovers them, but that choice is discrete and cannot be learned by gradient descent. Explorer (a) below animates this. Rectangular version: for free $X\in\mathbb R^{n\times n}$, $Y\in\mathbb R^{m\times n}$ and $Z=X-X^{\top}+Y^{\top}Y$, the matrix $\big[(I+Z)^{-1}(I-Z);\,-2Y(I+Z)^{-1}\big]$ has orthonormal columns (walkthrough, Step 3). Pauli et al. write the same map as $U=(I+M)^{-1}(I-M)$, $V=2Z(I+M)^{-1}$ with $M=Y-Y^{\top}+Z^{\top}Z$; the sign of the lower block is irrelevant.

Cayley convolutions. The 2-D DFT block-diagonalizes a circular convolution with $c$ channels on $n\times n$ images, $\mathrm{FFT}(\mathrm{conv}_W(X))[:,i,j]=\tilde W[:,:,i,j]\,\tilde X[:,i,j]$, and transposes or inverts it frequency by frequency. Trockman & Kolter skew-symmetrize per frequency ($\tilde A=\tilde W-\tilde W^{*}$), apply the Cayley transform to each $c\times c$ block and return with an exactly real inverse FFT. Price: $n^2$ inverses of $c\times c$ matrices, $O(n^2c^3)$ operations (see Primer 0), per layer and training forward pass (with frozen weights the transformed blocks can be cached; see "Inference cost" below), circular padding only, strides emulated by invertible downsampling.

SOC and BCOP

For skew-symmetric $J$, $\exp(J)^{\top}\exp(J)=\exp(-J)\exp(J)=I$. Skew Orthogonal Convolutions build a filter whose Jacobian is skew-symmetric ($L=M-\mathrm{conv\_transpose}(M)$ for an arbitrary filter $M$, their Theorem 2) and apply the convolution exponential truncated after $k$ terms (Singla & Feizi, ICML 2021).

BCOP builds orthogonal convolutions as block convolutions of an orthogonal matrix with symmetric projectors, $W=H\,\Box\,[P_1\ \ I-P_1]\,\Box\cdots\Box\,[P_{K-1}\ \ I-P_{K-1}]$, $P_i=P_i^2=P_i^{\top}$ (1-D case; see the fact below). It also shows why such parameterizations are incomplete: 1-D orthogonal convolutions with kernel size $K$ and $n$ channels form $2(K-1)n+2$ connected components, while a continuous map from $\mathbb R^N$ reaches only one (Li et al., NeurIPS 2019): $\mathbb R^N$ is path-connected and continuous images of path-connected sets are path-connected, so a continuous parameterization cannot jump between components, however many parameters it has (see Primer 0; the simplest case is $O(1)=\{-1,1\}$, which no continuous map from $\mathbb R$ covers).

AOL: almost-orthogonal layers

AOL needs no inverse, FFT or iteration; for convolutions the rescaling is channel-wise, so the trained layer is an ordinary convolution, and the learned $PD$ ends up almost orthogonal (Prach & Lampert, ECCV 2022). Section 3 shows it is one analytic solution of a single matrix inequality.

The cost of orthogonality

Expressivity: 1-Lipschitz layers compose to networks whose true constant is often far below the certified 1. Wang & Manchester's toy regression of a discontinuous function (their Table 1) measures tightness as empirical lower bound over certified bound:

Inference cost: Fourier-domain layers recompute their per-frequency Cayley inverses whenever the weights change, i.e. in training. With frozen weights and a fixed image size the transformed blocks can be computed once and cached (Wang & Manchester's code does this in evaluation mode), so an inference pass needs FFTs, per-frequency matrix products and inverse FFTs but no inverses, and no image-size kernel has to be materialized. LipKernel (Section 4) reports standard-form kernels two to three orders of magnitude faster in its benchmark, which matters inside a feedback loop.

3. SLL and Sandwich Layers: LMIs Solved by Construction

Section 2 bounded weights. The two constructions here take a LipSDP-type LMI, which also involves the activation through its multiplier, and solve it analytically: SLL layer by layer for residual blocks, the Sandwich parameterization for the whole network.

One inequality behind SN, orthogonal layers, AOL and CPL

Statement 1 is one line: $\|g(x)-g(y)\|^2=(x-y)^{\top}T^{-1/2}W^{\top}WT^{-1/2}(x-y)\le\|x-y\|^2$. Every known 1-Lipschitz layer is a choice of $T$: SN is $T=\|W\|^2I$ in Statement 1; orthogonal layers are $T=I$ with $W^{\top}W=I$; AOL is $T=\operatorname{diag}\big(\sum_j|W^{\top}W|_{ij}\big)$, for which $T-W^{\top}W$ is diagonally dominant with nonnegative diagonal, hence PSD; CPL (Meunier et al. 2022) is the SN choice in Statement 2; SLL is Statement 2 with $T_{ii}=\sum_j|W^{\top}W|_{ij}\,q_j/q_i$ for learnable $q_i\gt0$ (e.g. $q_i=e^{r_i}$ with free $r_i$). Here $|W^{\top}W|_{ij}$ is the absolute value of the $(i,j)$ entry. Since $T_{ii}\ge(W^{\top}W)_{ii}=\|w_i\|^2$ for the $i$-th column $w_i$ of $W$, this $T$ is nonsingular unless a column of $W$ vanishes; adding a fixed $\epsilon\gt0$ to every $T_{ii}$ removes that exception and keeps $T\succeq W^{\top}W$, because $\epsilon I\succeq0$.

Why: $T-QW^{\top}WQ^{-1}=Q(T-W^{\top}W)Q^{-1}$ is similar to the symmetric $T-W^{\top}W$ (a similarity $QMQ^{-1}$ keeps the eigenvalues themselves; a congruence $Q^{\top}MQ$ would keep only their signs), so it has the same real eigenvalues. Gershgorin's theorem puts every eigenvalue of a square matrix $M$ in one of the discs $|z-M_{ii}|\le\sum_{j\ne i}|M_{ij}|$; here the discs are centred at nonnegative diagonal entries with radii no larger than those entries, so the eigenvalues are $\ge0$. The $q_i$ are trained with the weights; $q_i\equiv1$ recovers AOL's choice of $T$, here used in the residual layer of Statement 2.

The Sandwich layer: a complete parameterization of LipSDP

Wang & Manchester consider $z_0=x$, $z_{k+1}=\varphi(W_kz_k+b_k)$ ($k=0,\dots,L-1$), $y=W_Lz_L+b_L$, with $W_k\in\mathbb R^{n_{k+1}\times n_k}$ and $\varphi$ piecewise differentiable and slope-restricted in $[0,1]$ (their Assumption 2.1), and write LipSDP with the bound $\gamma$ itself.

A block-bidiagonal factor $H=PP^{\top}$ automatically has the tridiagonal pattern; the difficulty is to make the middle diagonal blocks diagonal. Wang & Manchester build the blocks of $P$ as $\Psi_kA_k$ and $\Psi_kB_k$ with $[A_k^{\top};B_k^{\top}]$ from a Cayley transform, so that $A_kA_k^{\top}+B_kB_k^{\top}=I$ collapses each diagonal block to a diagonal matrix. Reading off the weights gives:

The converse is constructive ($\Psi_k=\Lambda_k^{1/2}$, $B_k=\tfrac12\Psi_kW_k\Psi_{k-1}^{-1}A_{k-1}^{-\top}$, $A_k$ from a Cholesky factor of $I-B_kB_k^{\top}$ times an orthogonal matrix chosen to avoid eigenvalue $-1$, then the inverse Cayley map). This recursion is the nonsingular case: on the boundary of the LMI a block $A_{k-1}$ can be singular (for one hidden layer, $\gamma=\Lambda_0=1$, $W_0=\sqrt2$, $W_1=0$ gives $B_0=1$, $A_0=0$), and the paper's Appendix D.3 then replaces $A_{k-1}^{-\top}$ by the pseudoinverse $(A_{k-1}^{\top})^{+}$, or by $I$ when $A_{k-1}=0$, which forces $W_k=0$; its sufficiency proof uses the matching generalized Schur complements. (For a singular pivot block $C$, $\begin{bmatrix}C&B^{\top}\\B&E\end{bmatrix}\succeq0$ iff $C\succeq0$, $(I-CC^{+})B^{\top}=0$ and $E-BC^{+}B^{\top}\succeq0$; the middle range condition, which says that the coupling only acts in directions where $C$ has energy, is what a pseudoinverse alone would miss; see Primer A.) So completeness is relative to LipSDP: every network LipSDP can certify is reached, not every $\gamma$-Lipschitz network. With $h_k=\sqrt2A_{k-1}^{\top}\Psi_{k-1}z_k$ the network becomes a composition of identical modules with decoupled parameters:

The proof is one identity (Exercise 13.3): with $U=\sqrt2\Psi^{-1}B$, $Y=\sqrt2A^{\top}\Psi$ and $\Lambda=\Psi^2$, $2\Lambda-Y^{\top}Y-\Lambda UU^{\top}\Lambda=2\Psi(I-AA^{\top}-BB^{\top})\Psi=0$, so the certificate sits exactly on the boundary of the PSD cone. $\Psi$ reshapes each activation channel ($\Psi\varphi(\Psi^{-1}v+b)$ still has slopes in $[0,1]$), while $B$ and $A^{\top}$ are the input and output halves of one matrix with orthonormal columns, "sandwiching" the nonlinearity.

Network-level coupling. $W_k$ depends on the parameters of layers $k$ and $k-1$ ("interlacing"). The paper's Proposition C.1 shows the network satisfies weighted layer-wise norm bounds whose product is at most $\gamma$, while the plain norms $\|W_k\|$ and their product can exceed 1 even for $\gamma=1$ (Remark C.2); on random draws we found products far above 1, and the explorer lets you find your own.

"Backpropagate through the solves" needs only the derivative of the defining equation: if a layer computes $v$ from $Mv=r$, then $(dM)\,v+M\,dv=dr$, so $dv=M^{-1}\big(dr-(dM)\,v\big)$, one more solve with the same matrix (see Primer B). Standard automatic differentiation uses this rule for linear solves rather than differentiating the elimination steps.

Convolutions and results. Circular, unstrided convolutions are evaluated in the Fourier domain ($s(\lfloor s/2\rfloor+1)$ parallel inverses of $q\times q$ complex matrices for $q$ channels and $s\times s$ images). SLL was slightly better on CIFAR-10 with a much larger model (41M versus 3M parameters), Sandwich about 4% better on CIFAR-100 (48M versus 118M) and, with last-layer normalization, 33.4% clean / 24.7% certified (see Primer E) at $\varepsilon=36/255$ on Tiny-ImageNet with 39M parameters versus 32.1% / 23.0% for the 1.1B-parameter SLL X-Large as quoted there (Wang & Manchester, ICML 2023); the SLL paper's own table lists 23.2% certified.

4. Lipschitz-Bounded CNNs: Cayley–Gramian and LipKernel (Pauli et al.)

Convolutions are where LipSDP-style certificates get expensive: as dense layers they are Toeplitz matrices growing with the signal length, in the Fourier domain they need circular padding and an inverse per frequency. Pauli et al. take a systems-theory route: a convolution is a finite impulse response (FIR) filter, a linear state-space system whose state is a window of past inputs. Its certificate is a dissipation inequality with quadratic storage, and the LMI size depends on kernel size and channels only, not on the signal length (Pauli, Gramlich & Allgöwer, L4DC 2023; Module 12).

Convolutional layers as FIR systems

A 1-D convolutional layer with $c_{i-1}$ input channels, $c_i$ output channels and kernel size $\ell_i$ computes $w^i_k=\varphi\big(b_i+\sum_{j=0}^{\ell_i-1}K^i_j\,w^{i-1}_{k-j}\big)$, $K^i_j\in\mathbb R^{c_i\times c_{i-1}}$. With the state $x^i_k\in\mathbb R^{(\ell_i-1)c_{i-1}}$ holding the last $\ell_i-1$ inputs (Pauli et al. CDC 2023, eqs. (5)–(6)):

$A_i$ is a block shift, hence nilpotent ($A_i^{\ell_i-1}=0$), and all learnable kernel entries sit in $\hat C_i$.

Here $\mathbb S^{c}$ is the set of real symmetric $c\times c$ matrices, and $\|v\|^2_{Q}:=v^{\top}Qv$; this is a squared norm, and $\{v:\|v\|_Q\le1\}$ an ellipsoid, only when $Q\succ0$, which the construction below guarantees through $Q_i=L_i^{\top}L_i$. At the flattening step, $N$ is the number of spatial positions: stacking the channel vectors $u_1,\dots,u_N$ of the positions one after another, $I_N\otimes Q=\operatorname{blkdiag}(Q,\dots,Q)$ and $\operatorname{vec}(u)^{\top}(I_N\otimes Q)\operatorname{vec}(u)=\sum_{k=1}^Nu_k^{\top}Qu_k$ (see Primer A; another flattening order permutes this matrix).

Fully connected layers use $W_i=\sqrt2\,\Gamma_i^{-1}V_i^{\top}L_{i-1}$, $L_i=\sqrt2\,U_i\Gamma_i$, $Q_i=L_i^{\top}L_i$, $\Lambda_i=\Gamma_i^{\top}\Gamma_i$, with $\Gamma_i=\operatorname{diag}(\gamma_i)$ and $[U_i;V_i]=\operatorname{Cayley}(Y_i,Z_i)$ (Thm 4): $U_i^{\top}U_i+V_i^{\top}V_i=I$, multiplied by $\sqrt2\,\Gamma_i$ on both sides, becomes $Q_i+\Lambda_iW_iQ_{i-1}^{-1}W_i^{\top}\Lambda_i=2\Lambda_i$, and a Schur complement gives the FC LMI. This is equivalent to the Sandwich parameterization (their Remark 6); what is new is the gain factor $L_i$, passed on to the next layer. Convolutions need one more idea because of the dynamics in the upper-left block.

Step 1: the storage matrix from a controllability Gramian

With $F_i:=\begin{bmatrix}P_i-A_i^{\top}P_iA_i&-A_i^{\top}P_iB_i\\-B_i^{\top}P_iA_i&Q_{i-1}-B_i^{\top}P_iB_i\end{bmatrix}$ the conv LMI reads $\begin{bmatrix}F_i&-\hat C_i^{\top}\Lambda_i\\-\Lambda_i\hat C_i&2\Lambda_i-Q_i\end{bmatrix}\succeq0$, the FC shape with $F_i$ in place of $Q_{i-1}$. We need $F_i\succ0$ for every parameter value.

Step 2: the kernel from a Cayley transform

Proof in one identity: with the multiplier $\Lambda_i=\Gamma_i^2$, the two terms of the Schur complement are $Q_i=L_i^{\top}L_i=2\Gamma_iU_i^{\top}U_i\Gamma_i$ and, because $L^F_iF_i^{-1}L_i^{F\top}=I$ for $F_i=L_i^{F\top}L^F_i$, $\Lambda_i\hat C_iF_i^{-1}\hat C_i^{\top}\Lambda_i=2\Gamma_iV_i^{\top}V_i\Gamma_i$. Multiplying $U_i^{\top}U_i+V_i^{\top}V_i=I$ by $\sqrt2\,\Gamma_i$ on both sides therefore gives $Q_i+\Lambda_i\hat C_iF_i^{-1}\hat C_i^{\top}\Lambda_i=2\Gamma_i^2=2\Lambda_i$, so $2\Lambda_i-Q_i-\Lambda_i\hat C_iF_i^{-1}\hat C_i^{\top}\Lambda_i=0$ and a Schur complement with respect to $F_i\succ0$ gives the LMI. Compare the sandwich layer: $\Psi$ became $\Gamma$, and the Cholesky factor of the Gramian-based $F_i$ plays the role of the input-side scaling.

All sizes are set by $c_{i-1},c_i,\ell_i$; the signal length never enters. On the MIT-BIH arrhythmia ECG data (5 classes, mean of 5 CNNs) the Lipschitz-bounded CNNs reached 84.8%, 90.0% and 94.7% test accuracy for $\rho=5,10,50$ (certified bounds 4.99, 9.85, 45.3), against 94.9% for a vanilla CNN with bound 147 and 92.4% / 90.6% for L2-regularized CNNs (training loss plus a weight penalty $\lambda\sum_k\|W_k\|_F^2$) with bounds 45.7 / 33.9; they degraded more slowly under $\ell_2$ PGD attacks (projected gradient ascent on the loss over an $\ell_2$ ball around the input: a successful attack disproves robustness, a failed one certifies nothing; see Primer E) and could be trained at bounds where L2-regularized training failed (Pauli, Wang, Manchester & Allgöwer, CDC 2023).

LipKernel: 2-D convolutions parameterized in standard form

LipKernel extends the construction to 2-D and 1-D CNNs with pooling, strides, dilation and zero padding, with the goal $\|\mathrm{NN}_\theta(u_a)-\mathrm{NN}_\theta(u_b)\|_Q\le\|u_a-u_b\|_R$ ($(Q,R)$-Lipschitz; $Q=I$, $R=\rho^2I$ is $\rho$-Lipschitz, other choices weight input directions and output classes) (Pauli, Wang, Manchester & Allgöwer, Automatica 2026).

Each layer type has its own LMI: $\varphi\circ$conv (the $3\times3$ block LMI above with $X_{k-1},X_k$ for $Q_{i-1},Q_i$), pooling$\circ\varphi\circ$conv (the $(3,3)$ block becomes $2\Lambda-\rho_p^2X$ with pooling constant $\rho_p$; $X$ diagonal for max pooling), $\varphi\circ$FC, and the last layer $X_{l-1}-W_l^{\top}QW_l\succeq0$. A 2-D convolution is realized as a Roesser model with a horizontal and a vertical state, $$\begin{bmatrix}x_1[i+1,j]\\x_2[i,j+1]\end{bmatrix}=\begin{bmatrix}A_{11}&A_{12}\\A_{21}&A_{22}\end{bmatrix}\begin{bmatrix}x_1[i,j]\\x_2[i,j]\end{bmatrix}+\begin{bmatrix}B_1\\B_2\end{bmatrix}u[i,j],$$ $$y[i,j]=C_1x_1[i,j]+C_2x_2[i,j]+Du[i,j]+b,$$ where $A_{21}=0$, $A_{11},C_1,A_{22},B_2$ are fixed shift structures and the kernel blocks fill $A_{12},B_1,C_2,D$ (Pauli, Gramlich & Allgöwer, IFAC-PapersOnLine 2024; the 2-D Lur'e view of CNNs is Gramlich et al., Automatica 2026), with block-diagonal storage $P=\operatorname{blkdiag}(P_1,P_2)$. 2-D Lyapunov equations have no general solution, so LipKernel constructs one for this FIR structure: since $A_{21}=0$ the $x_2$ dynamics are decoupled, $P_2^{-1}$ is a nilpotent Gramian-type sum, and $P_1^{-1}$ a second sum absorbing the coupling through the now free $A_{12},B_1$ (Lemma 16). The remaining blocks $[C_2\ D]$ come from a Cayley transform (Theorem 18), where the needed $2\Gamma-C_1F_1^{-1}C_1^{\top}\succ0$ is secured by $\gamma_i=\epsilon+\delta_i^2+\tfrac12\sum_j|C_1F_1^{-1}C_1^{\top}|_{ij}\,q_j/q_i$: SLL's Gershgorin trick. After training, the bijection $(A,B,C,D)\leftrightarrow K$ returns an ordinary kernel.

Accuracy as reported on MNIST (32×32; 20 epochs; mean of three initializations; certified if the margin exceeds $\sqrt2\rho\varepsilon$):

LipKernel is slightly less expressive than Sandwich but more flexible and much faster at inference, and it beats Orthogon and AOL on clean accuracy.

5. Recurrent Equilibrium Networks and R2DN

For dynamic models (system identification, i.e. fitting a dynamical model to input–output data; observers, which estimate an unmeasured state; feedback policies) the model itself must be stable for every parameter value, including those SGD visits on the way. A recurrent equilibrium network (REN) is an LTI system in feedback with an implicit (equilibrium) layer (Revay, Wang & Manchester, IEEE TAC 2024; recurrent and equilibrium layers: Primer E):

with $\sigma$ elementwise, piecewise differentiable and slope-restricted in $[0,1]$ (Assumption 1). With $D_{11}$ strictly lower triangular the equilibrium is computed row by row (acyclic REN); a deep feedforward network is the case of block sub-diagonal $D_{11}$, so RENs contain the networks of Section 3, all stable LTI systems and previously known contracting RNNs.

Why: take two trajectories with the same input, so $\Delta x_{t+1}=A\Delta x_t+B_1\Delta w_t$ and $\Delta v_t=C_1\Delta x_t+D_{11}\Delta w_t$. The quadratic form of the LMI matrix at $[\Delta x_t;\Delta w_t]$ equals $\bar\alpha^2|\Delta x_t|_P^2-|\Delta x_{t+1}|_P^2-2\Delta w_t^{\top}\Lambda(\Delta v_t-\Delta w_t)$, and the last term is $\le0$ by the diagonal incremental QC. The matrix is $\succeq\eta I$ for some $\eta\gt0$ (its smallest eigenvalue), so $|\Delta x_{t+1}|_P^2\le\bar\alpha^2|\Delta x_t|_P^2-\eta|\Delta x_t|^2\le\alpha^2|\Delta x_t|_P^2$ with $\alpha^2=\bar\alpha^2-\eta/\lambda_{\max}(P)$ (shrink $\eta$ if needed so that $\alpha\gt0$): an incremental Lyapunov function decreasing at a rate $\alpha\lt\bar\alpha$. Iterating and using $\lambda_{\min}(P)|z|^2\le|z|_P^2\le\lambda_{\max}(P)|z|^2$ gives $|\Delta x_t|\le\sqrt{\lambda_{\max}(P)/\lambda_{\min}(P)}\,\alpha^t|\Delta x_0|$: contraction with a constant $K$ that depends only on $P$, not on the initial states or the input.

From LMI to free parameters. The LMI is convex in $(P,\Lambda)$ for a fixed model but not jointly in the weights, so RENs use an implicit form with extra invertible $E$ and $\Lambda$, $Ex_{t+1}=Fx_t+\mathcal B_1w_t+\mathcal B_2u_t+\dots$ and $\Lambda v_t=\mathcal C_1x_t+\mathcal D_{11}w_t+\mathcal D_{12}u_t+\dots$ (so $A=E^{-1}F$, $C_1=\Lambda^{-1}\mathcal C_1$, and so on), in which contraction is an LMI $\mathcal H(\theta_{\rm cvx})\succ0$ jointly convex in weights, certificate and multiplier. The map $\theta_{\rm cvx}\mapsto\mathcal H$ is onto the positive definite cone and has an explicit right inverse, so writing $\mathcal H=X^{\top}X+\epsilon I$ with a free square $X$ and reading the model off its blocks gives a direct parameterization; for a fixed $\epsilon$ it reaches exactly the matrices $\mathcal H\succeq\epsilon I$, not every $\mathcal H\succ0$. Random free parameters give random contracting systems, a strictly larger family of stable echo state networks than previously known.

R2DN: removing the equilibrium solve

The implicit layer must be solved at every time step (e.g. by operator splitting), which is slow on GPUs for long sequences. R2DN (Barbara, Wang & Manchester, accepted to CDC 2026) replaces it by an explicit network, $w_t=\phi_g(C_1x_t+D_{12}u_t+b_v)$ with the same linear equations for $x_{t+1}$ and $y_t$, where $\phi_g$ is 1-Lipschitz for every value of its own parameters $g$, e.g. a stack of sandwich layers. The network enters the certificate only through $|\Delta w_t|\le|\Delta v_t|$, and the linear part is directly parameterized so that the interconnection is contracting and, if desired, $\gamma$-Lipschitz from $u$ to $y$ (their Prop. 1 covers general incremental IQCs; the direct parameterization of the Lipschitz case keeps $D_{22}=0$). Training and inference become up to an order of magnitude faster than for RENs with similar performance on system identification, observer design, feedback control and sequential image classification, scaling better with model size.

6. Certified Robust Accuracy: State of the Art 2026

For logits $f(x)\in\mathbb R^N$ and prediction $y=\arg\max_if_i(x)$, a Lipschitz certificate turns the margin into a radius: the deterministic guarantee type of Module 1.

Why $\sqrt2$: $f_y-f_j=(e_y-e_j)^{\top}f$, so $|\Delta(f_y-f_j)|\le\|e_y-e_j\|\,\|\Delta f\|\le\sqrt2L\|\delta\|$, and a label flip needs some margin to reach zero (Tsuzuku et al., NeurIPS 2018). Three conventions are in use: the global rule $M_f/(\sqrt2L)$ (BRONet, LipKernel); a constant $L_{yj}$ per logit difference (GloRo nets, Leino, Wang & Fredrikson, ICML 2021; LiResNet, Hu, Zou, Wang, Leino & Fredrikson, NeurIPS 2023), which for a 1-Lipschitz feature map and a linear last layer with rows $w_i$ can be $\|w_y-w_j\|$; and LipNeXt's decomposition into $N-1$ binary classifiers ("one-vs-rest" in the paper's words, which with $N-1$ classifiers we read as $f_y-f_j$), each certified with radius $|f_y-f_j|/K_j$ for a Lipschitz bound $K_j$, the overall radius being the minimum. The last two coincide when $K_j=L_{yj}$, and with a linear last layer $\|w_y-w_j\|$ is never worse than the $\sqrt2L$ rule based on the product certificate (Exercise 13.5).

Standard $\ell_2$ radii for images in $[0,1]$: $36/255\approx0.141$, $72/255\approx0.282$, $108/255\approx0.424$. CIFAR-10 certified robust accuracy (CRA) as tabulated by LipNeXt (Tables 1–2), no extra data unless marked:

BRONet uses block-reflector orthogonal layers $W=I-2V(V^{\top}V)^{-1}V^{\top}$ for full-column-rank $V\in\mathbb R^{m\times n}$: symmetric, orthogonal, with $n$ eigenvalues $-1$ and $m-n$ eigenvalues $+1$ (precisely the eigenvalue the Cayley transform cannot produce), without iterative approximation (their Prop. 1). The reason: $P_V=V(V^{\top}V)^{-1}V^{\top}$ is the orthogonal projector onto the column space of $V$, $P_V^{\top}=P_V=P_V^2$, so $W^{\top}W=I-4P_V+4P_V^2=I$, $Wx=-x$ on the $n$-dimensional column space and $Wx=x$ on its orthogonal complement. BRONet also trains with a logit-annealing loss: for logits $z$, true class $t$ and $p=\operatorname{softmax}\big((z-\xi e_t)/T\big)$ (temperature $T$, offset $\xi$), $\mathcal L_{\rm LA}=-T(1-p_t)^{\beta}\log p_t$ (their eq. (8); this $\beta$ is unrelated to $\beta$-Abs below). The factor $(1-p_t)^{\beta}$ shrinks the loss of examples that are already classified with a large margin ($0.01$ for $p_t=0.9$, $\beta=2$), so capacity goes to the others; the loss only changes which network training finds, the certificate still comes from the orthogonal layers (Lai, Huang, Kung & Chen, ICML 2025).

LipNeXt drops parameterizations: orthogonal weights are updated directly on the orthogonal manifold by a manifold Adam whose exponential map is replaced by a norm-adaptive truncated Taylor series ("FastExp"); the small loss of orthogonality this causes is removed by a polar retraction (SVD) once per epoch, plus a Lookahead step taken in the tangent space. Convolutions are replaced by a spatial-shift module (each selected circular shift is a norm-preserving permutation), the activation $\beta$-Abs (absolute value on a fraction $\beta$ of the coordinates, identity on the rest) contains MinMax up to an orthogonal change of basis, and scale does the rest (Hu, Hu & Fredrikson, ICLR 2026). On ImageNet at $36/255$ its 2B-parameter model reaches 57.0% clean / 41.2% CRA without generated data, against 49.3% / 37.6% for BRONet (86M) and 45.6% / 35.0% for LiResNet (51M), and 41.0% / 22.4% at $\varepsilon=1$; BRONet with 2M generated images reaches 40.7% CRA, so this comparison mixes data regimes.

Deterministic versus randomized-smoothing certificates

Randomized smoothing certifies a different object. With the base label classifier $C(x)=\arg\max_i f_i(x)$ (fixed tie-breaking) and a noise level $\sigma\gt0$, the smoothed classifier is $g(x)=\arg\max_c\mathbb P_{\varepsilon\sim\mathcal N(0,\sigma^2I)}(C(x+\varepsilon)=c)$, and its certificate is the theorem below (Cohen, Rosenfeld & Kolter, ICML 2019; see Module 15). With an off-the-shelf diffusion denoiser and a large pretrained classifier, Carlini et al. (ICLR 2023) report certified top-1 accuracy of 71.1% / 54.3% / 38.1% / 29.5% on ImageNet at $\ell_2$ radius 0.5 / 1.0 / 1.5 / 2.0 (best noise level per radius).

Walkthrough: Deriving the Sandwich Layer From the LMI

The central derivation of this module, for one hidden layer so that every block is visible: start from the LipSDP condition, turn it into a norm bound, satisfy the norm bound with Cayley factors, and read off the weights. The result is exactly Wang & Manchester's formula (8) for a single hidden layer, and the same three moves reappear in the Cayley–Gramian layers of Section 4. The Schur complement and congruence steps use only Module 2.

Interactive: Cayley Transform and a 1-Lipschitz Layer

Panel (a): the Cayley transform of $A=\begin{bmatrix}s&a\\-a&-s\end{bmatrix}$; for $s=0$, $A$ is skew-symmetric and $W=(I-A)(I+A)^{-1}$ is a rotation by $2\arctan a$, while any $s\ne0$ turns the unit circle into an ellipse. Panel (b): the 2-2-1 network $f(x)=W_1\,\mathrm{ReLU}(W_0x+b_0)$ of the walkthrough built from arbitrary free parameters (for $X_0$ only the skew part $X_0-X_0^{\top}$ enters the Cayley map, so the slider sets its entry $(X_0-X_0^{\top})_{12}$). On each region where the ReLU activation pattern $D=\operatorname{diag}(d_1,d_2)$, $d_i\in\{0,1\}$, is fixed, $f$ is affine with gradient $W_0^{\top}DW_1^{\top}$, so its exact Lipschitz constant is the largest of these gradient norms over the patterns that occur on open sets of inputs: all four when $W_0$ is invertible ($x\mapsto W_0x+b_0$ is then onto), those met along a line when $W_0$ has rank one (the code sorts the two kinks), one when $W_0=0$. The panel compares it with the largest ratio $|f(x)-f(x')|/\|x-x'\|$ over random input pairs, which can only underestimate it, and reports the smallest eigenvalue of the certificate $H(\gamma,\Lambda=\Psi^2)$: a floating-point check of the identity proved in the walkthrough, so values of order $-10^{-16}$ are round-off; the guarantee for all parameters is the proof, not this number. The comparison network skips the Cayley step: $W_0=\sqrt{2\gamma}\,\Psi^{-1}Y_0^{\top}$, $W_1=\sqrt{2\gamma}\,Y_1^{\top}\Psi$.

Exercises

Graded practice: build the calculation, then audit the claim

These twelve new questions each include an independent hint and a fully worked solution. The difficulty measures the amount of reasoning, not the amount of notation. The original research exercises follow below.

Easy: read definitions and compute

Medium: connect two or three steps

Hard: combine calculations with assumptions

Original research exercises

Continue here when the core calculations and the certificate assumptions are clear. Use the graded questions above as a warm-up; the original derivations below remain available in full.

Key Papers

Flashcards