12. Lipschitz Bounds via SDP: LipSDP and Beyond

Certified radii, slope-restricted activations as quadratic constraints, the LipSDP LMI, Pauli's training and estimation extensions, and the limits of SDP certificates

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

A two-number bridge to the LMI

A difference is just subtraction between two evaluations: $u=a-b$ and $w=\varphi(a)-\varphi(b)$. For ReLU, $a=2,b=-1$ gives $u=3,w=2$, a chord slope $2/3$. The allowed slope interval $[0,1]$ becomes $2w(u-w)=4\ge0$. This is one scalar fact that the matrix certificate will combine with the network equations.

Read the final inequality in two stages: $M\preceq0$ means $\xi^\top M\xi\le0$ for every vector $\xi$; the activation QC is nonnegative only for valid network increments. Subtracting that nonnegative QC from the matrix inequality leaves $\|\Delta f\|^2-\rho\|\Delta x\|^2\le0$. Hence $L=\sqrt\rho$. If matrix signs feel unfamiliar, review quadratic forms and PSD order and the S-procedure sign convention.

A Lipschitz bound is the simplest certificate a neural network can carry: one number $L$ with $\|f(x)-f(y)\|\le L\|x-y\|$ for all inputs. It turns a classifier's margin (see Primer E) into a deterministic robustness radius, bounds the incremental gain of a neural controller in small-gain arguments (see Primer D and, later, Module 14), and is the quantity that Lipschitz-by-design architectures enforce by construction (see Module 13). Computing the smallest such $L$ is NP-hard already for one hidden layer (Section 2; P, NP and NP-hardness are defined in Primer 0), while multiplying spectral norms (see Primer A) is cheap and often very loose. LipSDP sits in between: it describes every activation by the chord-slope inequality it satisfies, writes the Lipschitz question as a quadratic form, and lets the S-procedure of Module 2 turn it into a semidefinite program.

This page derives LipSDP line by line and then follows Patricia Pauli's work on it: the counterexample showing that only diagonal multipliers are valid, training under the LMI (ADMM, then log-det barriers), convolutions as 1-D and 2-D state-space systems (and the layer-wise GLipSDP), and quadratic constraints for activations that are not slope-restricted. It ends with the scalable reformulations and with what SDP certificates provably cannot do.

1. Lipschitz Constants and Certified Robustness

Why this one number certifies a classifier: every logit moves by at most $L\|\delta\|$ when the input moves by $\delta$, and a logit difference by at most $\sqrt2\,L\|\delta\|$. As long as the winning logit is ahead by more than that, nothing can overtake it.

Lipschitz-margin training. To make the certificate hold at radius $c$ on the training data, Tsuzuku et al. add $\sqrt2\,c\,L_F$ to every non-target logit before the cross-entropy loss (see Primer E), where $L_F$ is a cheap, differentiable upper bound of the network's Lipschitz constant (their LMT). For target class $y$ and logits $a$ the cross-entropy is $-a_y+\log\sum_je^{a_j}$; after the shift it is at least $\log2$ unless $a_y$ exceeds every other logit by more than $\sqrt2cL_F$. The network is thereby trained to have margins larger than $\sqrt2cL_F$, which is exactly the hypothesis of the theorem; training only encourages this, and the radius $c$ is certified on those inputs where the margin inequality holds after training. For context on what this buys at scale: the 2026 state of the art among deterministic $\ell_2$ certificates on CIFAR-10 is LipNeXt with $85.0\%$ clean and $73.2\%$ certified accuracy at $36/255$ without extra data (see Module 13). Those networks are $1$-Lipschitz by construction, so they need no estimation at all; LipSDP and its descendants matter whenever the architecture does not come with a certificate, and for control, where the Lipschitz bound of a trained controller is an input to a stability test.

2. The Product Bound and Why It Is Loose

The product bound pays for two things that never happen:

• Misaligned singular vectors. For nonzero $A,B$, equality $\|AB\|=\|A\|\|B\|$ requires a unit vector $v$ with $\|Bv\|=\|B\|$ whose image $Bv$ lies in the top right-singular subspace of $A$ (the directions $A$ stretches most). Diagonal slope matrices only rescale coordinate axes: they can change a direction ($\mathrm{diag}(1,0)$ maps $(1,1)$ to $(1,0)$) but cannot perform the arbitrary rotation that would realign two layers. Example (the one worked in the walkthrough): $W_0=\begin{bmatrix}1& 2\\1& -2\end{bmatrix}$, $W_1=\begin{bmatrix}1& 1\end{bmatrix}$, ReLU. $W_0$ stretches $e_2$ most, into the direction $(1,-1)$, but $W_1$ only responds to $(1,1)$. The product is $2\sqrt2\cdot\sqrt2=4$, while the four regional Jacobians $W_1DW_0=(d_1+d_2,\,2d_1-2d_2)$ have norms $0,\sqrt5,\sqrt5,2$: $L^\star=\sqrt5\approx2.24$.
• Impossible activation patterns. $f(x)=\mathrm{ReLU}(x)-\mathrm{ReLU}(-x)=x$, i.e. $W_0=\begin{bmatrix}1\\-1\end{bmatrix}$, $W_1=\begin{bmatrix}1& -1\end{bmatrix}$. The product is $\sqrt2\cdot\sqrt2=2$ but $L^\star=1$. Only the patterns $D=\mathrm{diag}(1,0)$ and $\mathrm{diag}(0,1)$ occur; $D=I$ would give slope $W_1W_0=2$ but needs $x\gt0$ and $-x\gt0$ at once. Every method that treats neurons independently, LipSDP included, pays for $D=I$ here.

SeqLip as an intermediate. Virmaux & Scaman start from the all-patterns bound with $[\alpha,\beta]=[0,1]$ (their eq. 7), which is itself hard to evaluate, and split it layer by layer using SVDs $W_i=U_i\Sigma_iV_i^{\top}$ (box below): each factor involves only two consecutive layers and captures their misalignment. For small layers ($\le20$ neurons) each factor's maximum is found by enumerating $0/1$ patterns; otherwise a power-method heuristic ("Greedy SeqLip") is used, which is not guaranteed to be an upper bound, as Fazlyab et al. point out. For one hidden layer, exact SeqLip coincides with the all-patterns bound, at a cost exponential in the width.

3. Slope Restriction as an Incremental Quadratic Constraint

Incremental versus non-incremental. An incremental QC relates two arbitrary inputs; it is the right object for global Lipschitz bounds and contraction. A non-incremental QC relates $v$ and $\phi(v)$ at a single input (with $\phi(0)=0$, the classical sector condition), or $v-v_*$ and $\phi(v)-\phi(v_*)$ around a fixed point $v_*$. Those are what stability analysis of an equilibrium uses, often in a local, much tighter form (see Module 14; Yin, Seiler & Arcak, TAC 2022), and what the reachability and verification SDPs DeepSDP (Fazlyab, Morari & Pappas, TAC 2022) and Reach-SDP (Hu et al., CDC 2020) use. The distinction decides which multipliers are admissible (Section 5).

4. The LipSDP Derivation

The whole method fits in one sentence: the Lipschitz inequality is a quadratic inequality in the increments, every neuron contributes a quadratic inequality in the same increments, and the S-procedure turns "the neuron inequalities imply the Lipschitz inequality" into a linear matrix inequality (Fazlyab, Robey, Hassani, Morari & Pappas, NeurIPS 2019). We derive it for one hidden layer, $f(x)=W_1\phi(W_0x+b_0)+b_1$ with $W_0\in\mathbb R^{n\times n_0}$ and $W_1\in\mathbb R^{m\times n}$, and then read off the deep version.

Fix two inputs $x,y$ and write $\Delta x=x-y$, $\Delta v=W_0\Delta x$ (the bias cancels), $\Delta z=\phi(W_0x+b_0)-\phi(W_0y+b_0)$ and $\Delta f=f(x)-f(y)=W_1\Delta z$. Everything below is a statement about the stacked increment $\xi=(\Delta x,\Delta z)\in\mathbb R^{n_0+n}$.

Step 1: the goal as a quadratic form. With $\rho=L^2$,

because $\|\Delta f\|^2=\Delta z^{\top}W_1^{\top}W_1\Delta z$. The output layer is linear, so it only contributes the block $W_1^{\top}W_1$.

Step 2: what the activations guarantee. The lemma of Section 3 applies to the pair $(\Delta v,\Delta z)$. Since $(\Delta v,\Delta z)=\mathrm{blkdiag}(W_0,I)\,\xi$, a congruence moves the constraint onto $\xi$:

for every $T\in\mathcal T_n$ and every pair of inputs.

Step 3: the S-procedure. Call the two matrices $G(\rho)$ (goal) and $Q(T)$ (constraint). If $G(\rho)+Q(T)\preceq0$ for some $T\in\mathcal T_n$, then for every increment the network can produce,

The multipliers $\lambda_i$ on the diagonal of $T$ are exactly the S-procedure multipliers of Module 2, one per neuron. The implication runs one way only: the S-procedure with $n\ge2$ constraints is sufficient, not necessary, and this is one of the places where conservatism enters.

Step 4: assemble. Adding the blocks gives the LipSDP matrix

The deep network. For $l$ hidden layers, stack $\xi=(\Delta z_0,\Delta z_1,\dots,\Delta z_l)\in\mathbb R^{n_0+n}$ with $\Delta z_0=\Delta x$. Two selection matrices produce the pre-activation and activation increments of all hidden layers at once:

With $T=\mathrm{blkdiag}(T_1,\dots,T_l)$ diagonal, Steps 1–4 go through unchanged and give Theorem 2 of LipSDP:

Written out for two hidden layers, the matrix is block tridiagonal, because each neuron's constraint couples only its own input and output increments, which live in consecutive layers:

This banded structure is the single most important fact for everything that follows: Pauli et al. factor it layer by layer (Section 6), the dissipativity view reads it as a chain of layer-wise storage functions (Section 7), and chordal and compositional methods decompose it (Section 9).

Splitting is valid because the Lipschitz constant of a composition is at most the product of the constants of its parts, but it reintroduces a product bound at every cut. In their MNIST experiments (five hidden layers, test accuracy above $97\%$), LipSDP-Neuron was tighter than CPLip and SeqLip.

5. What Went Wrong With Coupled Multipliers

The NeurIPS 2019 version of LipSDP (arXiv v1) stated its Lemma 1 for a larger multiplier class than the diagonal one,

on the grounds that all neurons apply the same function $\varphi$ ("repeated nonlinearities"). The resulting variant with $O(n^2)$ multipliers (at most a constant times $n^2$ for large $n$; see Primer 0, which also covers the other use, $O(\varepsilon^2)$ as $\varepsilon\to0$, in Section 6), LipSDP-Network, was presented as the most accurate one. It is wrong. Pauli, Koch, Berberich, Kohler & Allgöwer (arXiv May 2020; IEEE L-CSS 2022 and ACC 2021) gave two counterexamples, and the current arXiv version of LipSDP (v2, January 2023) states Lemma 1 and Theorems 1–2 for diagonal $T$ only and lists only LipSDP-Neuron and LipSDP-Layer.

What a coupling term asserts. For $\alpha=0$, $\beta=1$ the term $\lambda_{12}(e_1-e_2)(e_1-e_2)^{\top}$ contributes

i.e. it treats the pair $(\Delta v_1-\Delta v_2,\ \Delta z_1-\Delta z_2)$ as if it were a chord of $\varphi$ with slope in $[0,1]$. But

is a difference of two chords taken at different places of the graph (between $\bar v_1,v_1$ and between $\bar v_2,v_2$), and nothing constrains its ratio to $\Delta v_1-\Delta v_2$.

Where the proof broke. Version 1 proved the incremental lemma by fixing $x$ and applying a non-incremental repeated-nonlinearity lemma to the shifted map $\tilde\phi(x,\delta)=\phi(x+\delta)-\phi(x)$. Its $i$-th component is $\delta_i\mapsto\varphi(x_i+\delta_i)-\varphi(x_i)$. Each component is slope-restricted and vanishes at $\delta_i=0$, so the diagonal terms are fine. But component $i$ and component $j$ are different functions unless $x_i=x_j$, so the repeated-nonlinearity lemma does not apply to them. In the words of Wang & Manchester (ICML 2023, Remark A.2), who credit the explanation to Revay, Wang & Manchester (2020): the nonlinearities are repeated as functions of $v_i$, but their increments are not repeated as functions of $\Delta v_i$.

What survives. Coupling is legitimate whenever the two compared quantities really are a chord of one function. At a single input, $\big(v_i-v_j,\ \varphi(v_i)-\varphi(v_j)\big)$ is a genuine chord of $\varphi$ between the points $v_i$ and $v_j$, so for ReLU

is valid. Such non-incremental repeated-nonlinearity terms are used, correctly, in DeepSDP (Fazlyab, Morari & Pappas, TAC 2022) and Reach-SDP (Hu, Fazlyab, Morari & Pappas, CDC 2020), where the input set is fixed and each constraint relates quantities at the same input, and in Appendix A.8 of LipSDP v2, which bounds $\|f(x)-f(y)\|\le\sqrt{\rho(y)}\|x-y\|$ for one fixed point $y$. Linear operations are also safe: in the 1-D CNN paper, average pooling admits full multiplier matrices while max pooling requires diagonal ones (Section 7).

6. Training Under Lipschitz Constraints: ADMM and Barriers (Pauli et al.)

Estimating $L$ after training is passive: if the bound is bad, one can only retrain and hope. Pauli and co-authors turned the certificate into a training constraint, $\min_W\ \mathcal L(W)$ subject to "the LMI holds". Two obstacles have to be removed first: $M(\rho,T)$ is not affine in the weights ($W_1^{\top}W_1$ is quadratic and $TW_0$ is bilinear in the unknowns; see Primer B), and stochastic gradient descent (see Primer B) does not handle semidefinite constraints.

Making the certificate convex in the weights. Following Pauli et al. (L-CSS 2022), fix $T$, set $\alpha=0$, and take a Schur complement to move the last layer out of the quadratic term. For one hidden layer (their Remark 4, with $L^2$ in place of $\rho$):

Because the corner block $-I$ is negative definite, the Schur complement of Module 2 gives

which is $M(L^2,T)$ with $\alpha=0$. Every entry of $M_{\rm tr}$ is affine in $(L^2,W_0,W_1)$, so for fixed $T$ the set of certified weights is convex. The multi-layer version (their eq. (7)) has the same form, with $[A;B]^{\top}\begin{bmatrix}0 & \beta T\\ \beta T & -2T\end{bmatrix}[A;B]$ in place of the first two block rows and the last layer moved into a $-I$ corner. Two choices make this work:

• $\alpha=0$ (their Remark 3). With $\alpha\gt0$ the block $-2\alpha\beta W_0^{\top}TW_0$ is matrix-concave in $W_0$, and a Schur complement cannot remove it because it enters with the wrong sign. $\alpha=0$ is the tight value for ReLU, tanh and sigmoid (slopes $[0,1]$, $[0,1]$, $[0,\tfrac14]$). It costs something only for activations such as leaky ReLU.
  Derivation — why $-W_0^{\top}TW_0$ is matrix-concave

  For $T\succeq0$, the map $W\mapsto W^{\top}TW$ is matrix-convex: for $0\le s\le1$ and any $W,V$,

  $$sW^{\top}TW+(1-s)V^{\top}TV-\big(sW+(1-s)V\big)^{\top}T\big(sW+(1-s)V\big)=s(1-s)(W-V)^{\top}T(W-V)\ \succeq\ 0 .$$

  So its negative is matrix-concave, and the set where a matrix-concave expression (plus affine terms) is $\preceq0$ need not be convex. With $\alpha=0$ the term is absent.

• $T$ fixed. $TW_0$ is bilinear, so optimising $T$ and $W$ jointly gives a bilinear matrix inequality. Pauli et al. fix $T$, e.g. from a LipSDP run on a network trained without the constraint or with L2 regularisation, and accept the conservatism; their Remark 5 notes that $T$ could be updated in an additional ADMM block at extra cost. The CDC 2022 paper below instead trains $T$ together with the weights and lets the barrier handle the bilinear inequality directly; its Remark 9 adds that the substitution $\tilde W_{k-1}=T_kW_{k-1}$ would make the constraint linear in $(\tilde W,T)$, at the price of a nonlinear reparameterisation of the network.

ADMM (L-CSS 2022). Duplicate the weights, $W=\bar W$, so that the loss sees $W$ and the certificate sees $\bar W$:

where the indicator $\iota(\cdot)$ is $0$ when its condition holds and $+\infty$ otherwise, and $\mu\gt0$ trades accuracy against the bound. With the augmented Lagrangian (see Primer A) $\mathcal L_\eta=\mathcal L(W)+\mu L^2+\iota(\cdot)+\mathrm{tr}\big(Y^{\top}(W-\bar W)\big)+\tfrac\eta2\|W-\bar W\|_F^2$ (the paper calls the penalty $\rho$; we write $\eta$ because $\rho=L^2$ here), ADMM alternates a loss step (backpropagation), a Lipschitz step (an SDP) and a dual update:

Here $W$ stands for all layers $(W_0,\dots,W_l)$, which have different shapes: $\|W-\bar W\|_F^2=\sum_j\|W_j-\bar W_j\|_F^2$ and $\mathrm{tr}\big(Y^{\top}(W-\bar W)\big)=\sum_j\mathrm{tr}\big(Y_j^{\top}(W_j-\bar W_j)\big)$ with one dual matrix $Y_j$ per layer, while the subscript $k$ above counts iterations. The loss step is solved only approximately, by stochastic gradient steps. The Lipschitz step follows from completing the square, $\mathrm{tr}(Y^{\top}(W-\bar W))+\tfrac\eta2\|W-\bar W\|_F^2=\tfrac\eta2\|\bar W-W-Y/\eta\|_F^2-\tfrac1{2\eta}\|Y\|_F^2$. It projects $W_{k+1}+Y_k/\eta$ onto the set of weights with a certified bound, traded against $\mu L^2$. That set is described by an LMI in $\bar W$, not by $\bar W\succeq0$, so there is no eigenvalue-clipping shortcut: every iteration calls an SDP solver. To enforce a prescribed bound, drop $\mu L^2$ and fix $L=L_{\rm des}$. Their Theorem 6 then guarantees that $L_{\rm des}$ bounds the network with weights $\bar W$ at every iteration, because every $\bar W$ is produced by a feasible SDP. The guarantee is about $\bar W$, so $\bar W$ is the network to deploy. Convergence of ADMM is not guaranteed for the non-convex loss; each subproblem is well behaved (SGD for the first, a convex SDP with a unique minimiser for the second).

An eleven times smaller certified bound for $0.2$ percentage points of accuracy; weight decay reaches a similar bound only by giving up six points. The price is a factor of about $45$ in training time, and the SDP step dominates as the network grows.

Log-det barriers (CDC 2022). Pauli, Funcke, Gramlich, Msalmi & Allgöwer (CDC 2022) replace the SDP subproblem by an interior-point penalty that ordinary backpropagation can handle. They write the certificate in positive-definite form (the negative of $M_{\rm tr}$, with $\beta=1$ and $\Lambda$ for $T$). For one hidden layer:

with barrier weight $\mu\gt0$ (their $\rho_j$), and the barrier set to $+\infty$ where $N\not\succ0$. As $\mu\to0$ the barrier vanishes at each fixed strictly feasible point, so optima on the boundary can be approached from inside; as $N$ approaches singularity, $-\log\det N\to+\infty$ and pushes the iterate back. A finite step can still jump across the boundary, which is why Algorithm 1 below tests every step.

Treating $\Lambda$ as fixed gives the "linear" barrier method (convex constraint); training it as well gives the "bilinear" one. The template is not limited to Lipschitz bounds. Their Theorem 4 certifies any incremental property

with an LMI of the same form. $Q=L^2I$, $S=0$, $R=-I$ is Lipschitz continuity; note that the vector is ordered (input, output) here, the reverse of the QSR convention in Module 2. Their Remark 5 points to closed-loop stability of feedback loops containing a network as another semidefinite constraint (see Module 14).

On MNIST ($14\times14$ inputs, hidden layers of $100$ and $30$ tanh neurons, $L\le20$) the bilinear barrier method reached $98.2\%$ with mean bound $16.7$ at $1.74$ times the nominal training time, where ADMM needed $104.4$ times (normalised to its own framework). The same machinery trains Wasserstein-GAN critics with a hard 1-Lipschitz constraint. There, weight clipping kept the bound far below $1$ (overly conservative), while with a gradient penalty the estimated bound was far above $1$, which suggests that the constraint was violated.

7. CNNs as Dynamical Systems: 1-D, Roesser and GLipSDP

A convolutional layer (see Primer E) acting on a signal of length $N$ can be unrolled into a Toeplitz matrix (see Primer A) and handed to LipSDP; the CDC 2022 paper did exactly this for its WGAN critic. The unrolled matrix has $cN$ rows, most of its entries are zero, LipSDP needs one multiplier per channel and position, and the resulting bound is valid only for that one input size. Control theory offers a better description: a convolution is a shift-invariant linear system, so it has a state-space realisation whose size depends on the channels and the kernel, not on $N$. Lipschitz estimation then becomes a dissipativity problem with a storage function (Module 2; Willems 1972).

Layer-wise incremental dissipativity. Run the layer on two input signals and write $\Delta$ for the differences. Give the layer an input weight $X_{\rm in}\succeq0$, an output weight $X_{\rm out}\succeq0$, a storage matrix $P\succeq0$ and a diagonal $\Lambda\succeq0$ for the activation. The layer LMI of the 1-D CNN paper (Corollary 10 there, written with $X=-\tilde Q$ and slopes $[0,1]$) is

Fully connected layers use the same chain with $A=B=C=0$, $D=W$: $\begin{bmatrix}X_{k-1} & -W^{\top}\Lambda\\ -\Lambda W & 2\Lambda-X_k\end{bmatrix}\succeq0$. Pooling layers need their own incremental QCs. The paper takes the stride equal to the window length $\ell$, so the windows do not overlap. Averaging over one window is then $1/\sqrt\ell$-Lipschitz (Proposition 7: the averaging row $\tfrac1\ell\mathbf 1^{\top}$ has norm $\sqrt\ell/\ell$), and so is the whole layer, which acts on disjoint windows: square the per-window bounds and add. Being linear, average pooling admits full multiplier matrices; max pooling, $1$-Lipschitz under the same assumption (per window, $|\max_ia_i-\max_ib_i|\le\max_i|a_i-b_i|\le\|a-b\|$), admits only diagonal ones, the same distinction as in Section 5. Overlapping windows need a bound on the whole pooling operator: averaging the adjacent pairs of $(u_1,u_2,u_3)$ is the matrix $\tfrac12\begin{bmatrix}1&1&0\\0&1&1\end{bmatrix}$, with norm $\sqrt3/2\gt1/\sqrt2$. For fully connected networks the layer-wise LMIs are equivalent to LipSDP (Remark 12, via chordal sparsity). For CNNs they are slightly more conservative than LipSDP on the unrolled network, but they hold for every input length and use multipliers of size $c$ instead of $cN$. In the experiments, LipSDP on the unrolled network was a little tighter for short signals, but its computation time grew quickly with the signal length, while the state-space SDP of a fully convolutional network does not grow at all. A fully connected head after flattening still sees all $c_pN_p$ features.

Two-dimensional convolutions. Images are indexed by two coordinates, so the natural state space is two-dimensional. Gramlich, Pauli, Scherer & Ebenbauer (Automatica 2026) represent a 2-D convolutional layer as a Roesser model

with a horizontal state $x_1$ and a vertical state $x_2$. A fully convolutional CNN is then a 2-D Lur'e system (see Primer D): a Roesser system in feedback with the activations. The dissipation inequality becomes $V_1(x_1(i_1{+}1,i_2))+V_2(x_2(i_1,i_2{+}1))\le V_1(x_1(i_1,i_2))+V_2(x_2(i_1,i_2))+s(u,y)$ with $V_j(x)=x^{\top}P_jx$. Summed over an image with zero initial states, the $V_1$ terms telescope along every row and the $V_2$ terms along every column; only nonnegative terms at the far edges remain, so the supply summed over the image is nonnegative. This yields an LMI in $(P_1,P_2)$ and a Lipschitz bound valid for images of any size. The arXiv version of the paper also lists Allgöwer as an author.

GLipSDP: one LMI per layer, any architecture. Pauli, Gramlich & Allgöwer (arXiv 2405.01125, v2 2024; a preprint without a journal record so far) read the network itself as a time-varying system whose "time" is the layer index, $u_{k+1}=y_k=\mathcal L_k(u_k)$. The exact Lipschitz constant then satisfies a dynamic-programming recursion on value functions (see Primer D):

Relax the equalities to inequalities $V_{k-1}\ge V_k\circ\mathcal L_k$, which still carry the output bound backwards, and the value functions to quadratics $V_k=\Delta^{\top}X_k\Delta$. The layer's QC and the S-procedure then turn each inequality into one LMI, $G_k(X_k,X_{k-1},\nu_k)\succeq0$, where $\nu_k$ collects the layer's other certificate variables (activation multipliers and, for state-space layers, storage matrices), and the bound comes from the SDP $\min\gamma^2$ s.t. $X_0=\gamma^2I$, $G_k\succeq0$, $X_l=I$. Lemmas supply $G_k$ for fully connected, 1-D/2-D/N-D convolutional (Roesser), state-space-model, residual, average and max pooling, slope-restricted and GroupSort layers. For GroupSort, $X_k$ is constrained to the group structure of Section 8 with $0\preceq X_k\preceq X_{k-1}$. Theorem 1 of the paper shows that for fully connected networks the layer-wise LMIs are equivalent to LipSDP, so the splitting adds no conservatism, and Remark 8 identifies them with LipSDP's chordal decomposition. The remaining sources of conservatism are the QCs themselves, structure restrictions on $X_k$, boundary effects of treating convolutions on infinite signals, and quadratic value functions. The authors report bounds tighter than LipLT and the weight-norm product on MNIST and CIFAR-10 CNNs and 18-layer ResNets, including networks for which LipSDP on the unrolled convolutions runs into memory issues.

8. Beyond Slope Restriction: GroupSort, MaxMin, Householder

Many Lipschitz-by-design networks avoid elementwise activations. Anil, Lucas & Grosse (ICML 2019) showed that 2-norm-constrained networks with ReLU, tanh or sigmoid activations cannot even represent the absolute value function, because these activations cannot preserve the gradient norm. They proposed GroupSort: split the $n$ pre-activations into $N=n/n_g$ groups of size $n_g$ and sort each group. MaxMin is the case $n_g=2$, $(x_1,x_2)\mapsto(\max,\min)$. GroupSort is 1-Lipschitz, its Jacobian is a permutation matrix almost everywhere ("gradient norm preserving"), and with $\infty$-norm-constrained weights it gives universal approximators of Lipschitz functions; whether 2-norm-constrained GroupSort networks are universal was left open by Anil et al. (see Module 13). The Householder activation (see Primer A) applies, per group, $x\mapsto x$ if $v^{\top}x\gt0$ and $x\mapsto(I-2vv^{\top})x$ otherwise, with $\|v\|=1$ learnable; $v=(1,-1)/\sqrt2$ recovers MaxMin. None of these is elementwise, so assumption (A1) of LipSDP fails.

Why the obvious workaround fails. MaxMin is a residual ReLU network,

so LipSDP applies to it. Pauli, Havens, Araujo, Garg, Khorrami, Allgöwer & Hu (ICLR 2024) show that the resulting SDP returns $\rho=2$: it certifies $\sqrt2$ for a map that is exactly 1-Lipschitz. The reason is the impossible-pattern cost from Section 2 in its purest form. The two ReLUs receive $u$ and $-u$ and can never be active together, but the slope QC allows $D=I$, for which the Jacobian is $H+GW=\begin{bmatrix}1 & 0\\ 1 & 0\end{bmatrix}$ with norm $\sqrt2$. What is needed is a constraint that describes the group as a whole.

LipSDP-NSR. Replacing the slope multiplier $\begin{bmatrix}0 & T\\ T & -2T\end{bmatrix}$ by $\begin{bmatrix}T-2S & P+S\\ P+S & -T-2P\end{bmatrix}$ in the LipSDP matrix of Section 4 gives the paper's Theorem 1: if $M(\rho,P,S,T)\preceq0$, the GroupSort (or Householder) network is $\sqrt\rho$-Lipschitz. In practice the authors set $S=P=0$, which in their tests gave the same bounds. For one hidden layer, $f(x)=W_1\phi(W_0x+b_0)+b_1$, the matrix then becomes block diagonal,

This is a weighted product bound: $T=\|W_1\|^2I$ recovers $\|W_0\|^2\|W_1\|^2$, and a better $T$ within the group structure gives less. For MaxMin alone ($W_0=W_1=I$), $T=I$ gives $\rho=1$, the correct constant that the residual-ReLU route missed. The multipliers $\gamma_i$ are a free adjustment along $\mathbf 1$ inside each group: sum preservation makes input and output agree in that direction, so weight placed there costs nothing.

FGL ("formal global Lipschitz") is the largest Jacobian norm over all combinations of within-group sorting permutations, including combinations that no input realises; it is an upper bound, computable only for small networks, while the sampled column is a lower bound, so the true constant lies between the two.

These networks were trained with $\ell_2$ adversarial training (see Primer E) and reach more than $96\%$ accuracy. Three observations. The residual-ReLU route can be worse than the naive product: the both-active pattern it cannot exclude costs up to a factor $\sqrt2$ at each MaxMin. The group-level constraint beats the product at every depth. And the gap between certified upper bounds and the sampled lower bound grows with depth, to a factor of about $65$ at 18 layers; Section 9 explains why no polynomial-time method can be expected to close it in general. The paper also extends the constraint to residual networks (generalising the residual LipSDP condition behind SLL layers; Araujo et al., ICLR 2023), to implicit models, to $\ell_\infty\to\ell_1$ bounds and, via the state-space view, to CNNs. GLipSDP uses the same structure: for a GroupSort layer it restricts $X_k$ to the group form and requires $0\preceq X_k\preceq X_{k-1}$.

9. Scalable Variants and Hardness Results

A generic SDP solver treats $M(\rho,T)$ as a dense matrix. That is what limits the plain method to a few thousand neurons (Section 4). Every scalable variant exploits the block-tridiagonal structure, and each one is exact in a precise sense or loses something identifiable.

Chordal decomposition. Xue, Lindemann, Robey, Hassani, Pappas & Alur (CDC 2022) use a chordal block sparsity pattern containing all nonzeros of $M(\rho,T)$, with maximal cliques given by pairs of consecutive layers, and the decomposition theorem below. Chordal-LipSDP replaces the one large LMI by one small LMI per clique and is equivalent to LipSDP (their Theorem 3); on deep networks it is much faster. Its tunable sparsity parameter $\tau\gt0$ adds neuron couplings and is invalid after Section 5; only $\tau=0$ certifies. The layer-wise chain of Section 7 (1-D CNN paper, Remark 12; GLipSDP, Remark 8) is exactly this decomposition, with the clique blocks written as value functions $X_k$.

ECLipsE: eliminate the value functions in closed form. In the chain of Section 7, the largest admissible $X_k$ for given multipliers is the Schur complement $X_k=2T_k-T_kW_{k-1}X_{k-1}^{-1}W_{k-1}^{\top}T_k$, and a larger $X_k$ only makes the next layer's LMI easier. Xu & Sivaranjani (NeurIPS 2024) run this recursion forward. Normalising $X_k=\rho M_k$ and $T_k=\rho\Lambda_k/2$ turns it into their recursion, in our indexing

Their Theorem 3 states that, for slopes in $[0,1]$ and given diagonal $\Lambda_k$, the strict LipSDP certificate holds if and only if every $M_k\succ0$ and $M_l-FW_l^{\top}W_l\succ0$. So the decomposition itself is exact, and the only approximation is how the $\Lambda_k$ are chosen: one layer at a time instead of jointly. ECLipsE picks $\Lambda_k$ by a small SDP of layer size that looks one layer ahead. ECLipsE-Fast uses $\Lambda_k=\lambda_kI$ with the closed form $\lambda_k=2/\lambda_{\max}\big(W_{k-1}M_{k-1}^{-1}W_{k-1}^{\top}\big)$ and needs no SDP at all. The authors report speed-ups of up to several thousand times on deep networks with estimates comparable to LipSDP; ECLipsE-Fast is a feasible point of LipSDP-Layer, hence never below it. On the walkthrough network, ECLipsE-Fast gives $\lambda_1=1/4$ and $L\approx3.02$, between LipSDP ($2.31$) and the product ($4$). The same Schur recursion is the block Cholesky factorisation of Pauli et al.'s barrier method (Section 6).

EP-LipSDP: an eigenvalue problem for first-order methods. Wang, Hu, Havens, Araujo, Zheng, Chen & Jha (ICLR 2024) move the LMI into the objective. The naive penalty $\zeta+c\,\lambda_{\max}^+(\cdot)$ is exact only if $c$ is large enough, and for LipSDP it was not known how large. Their trick is a redundant constraint. Slope restriction and $\|\Delta x\|\le1$ imply $\Delta z_i^2\le\|w_i\|^2=:a_i$ ($w_i^{\top}$ the $i$-th row of $W_0$), which gives the Shor relaxation an explicit trace bound $1+\|\Delta x\|^2+\|\Delta z\|^2\le2+\sum_ja_j$. With a trace bound, the penalty weight equal to that bound is exact. For a scalar-output one-hidden-layer network $f(x)=v^{\top}\phi(W_0x+b_0)+b_1$ with slopes in $[0,1]$:

The letters are the paper's: $\mathrm{diag}(\tau)$ is the slope multiplier ($T$ of Section 4), $\lambda$ multiplies the redundant constraints $\Delta z_i^2\le a_i$ and $\gamma$ the constraint $\|\Delta x\|^2\le1$; neither is the $\lambda_i$ of $T$ or GLipSDP's bound $\gamma$. Their Theorem 1: the optimal value equals that of the SDP with the redundant constraints and equals $2\sqrt{\rho^\star}$, twice the LipSDP bound. The problem has only sign constraints, so subgradient or bundle methods with autodiff and a Lanczos routine for $\lambda_{\max}$ apply. The implementation, LipDiff, scaled SDP-based Lipschitz estimation to networks on ImageNet, far beyond the memory limits of interior-point solvers. Theorem 2 of the paper covers multi-layer networks.

Cheaper ingredients. The product bound itself can be improved per layer: for convolutions, Gram iteration gives fast, differentiable and certified upper bounds on the spectral norm (Delattre, Barthélemy, Araujo & Allauzen, ICML 2023). For local certificates around one input, SDP-CROWN imports SDP-style coupling between neurons into bound propagation (Chiu, Chen, Zhang & Zhang, ICML 2025; see Module 15). Bound propagation pushes guaranteed enclosures of every neuron's possible values, for one specified input set, through the network layer by layer; SDP-CROWN tightens those enclosures on $\ell_2$ balls. These are per-input certificates and play no role in the global LipSDP derivation.

What no polynomial-time method can do. Exact computation of $L^\star$ is NP-hard already for one hidden layer (Section 2). Kuang, Xu, Sivaranjani & Lin (arXiv 2603.28113, v2 May 2026) locate the difficulty. Estimating a network's Lipschitz constant requires knowing which hidden states are reachable, and reachability is NP-hard. Through explicit constructions they show that this blindness can force SDP-based bounds to inherit the qualitative failures of the product bound, including polynomial conservatism per layer (box below). They argue that the difficulty afflicts every instance of the verification problem, not only worst-case reductions, so SDP is not sufficient: their Theorem 2.7 (assuming P $\ne$ NP) gives every ReLU network with more than three hidden layers a twin of width one larger that either computes the same function or has an arbitrarily large Lipschitz constant, depending on an infinitesimal change of two biases, and that no polynomial-time verifier can tell apart from the original. Read it carefully: it says that each network has a hard-to-distinguish twin, not that computing $L^\star(f)$ is hard for each fixed $f$. They also argue SDP is not necessary: several failures of the product bound come from removable parameterisation pathologies, which optimising or regularising the product bound itself can mitigate, and for a new bias-free trigonometric layer $x\mapsto A\cos(Wx)+B\sin(Wx)$ with $W$ of full row rank the product bound is tight up to a factor $\sqrt2\,\kappa(W)$, the condition number (their Theorem 3.1). This page has shown the mechanism several times: the slope QC forgets which patterns occur together (Section 3); the certificate ignores the biases that decide reachability (Section 4); MaxMin's $u$ and $-u$ pair (Section 8); and the growing gap in the 18-layer table.

Walkthrough: LipSDP for One Hidden Layer

The derivation of Section 4, carried out by hand for the network $f(x)=W_1\mathrm{ReLU}(W_0x+b_0)+b_1$ with $W_0=\begin{bmatrix}1 & 2\\ 1 & -2\end{bmatrix}$ and $W_1=\begin{bmatrix}1 & 1\end{bmatrix}$ (any biases), from the neuron constraints to the optimal multiplier and a comparison with the product bound and the true constant.

Interactive: Naive vs LipSDP vs Empirical

A 2-3-1 ReLU network with Gaussian weights and biases (see Primer C) from a seeded generator. For the same network the explorer computes six numbers. The product bound uses spectral norms from a Jacobi eigen-solver. LipSDP-Layer and LipSDP-Neuron are approximated numerically on the $5\times5$ LMI by an interior-point method: the log-det barrier of Section 6, minimised by Newton steps while $\mu\to0$ (Section 6 uses gradient steps because it trains the weights). The numerical LipSDP value is a point where $-M(\rho,T)\succ0$ passes a Cholesky test. The all-patterns bound is the largest $\|W_1DW_0\|$ over all eight $0/1$ patterns. The enumerated global estimate of $L^\star$ is the largest gradient over the patterns that some input realises, found by solving one small linear program per pattern (see Primer B). The sampled lower bound is the largest gradient norm at $4000$ random inputs in the square $[-3,3]^2$. The slider $\kappa$ blends neuron 3 into the negative of neuron 1 (weights, bias and output weight); at $\kappa=1$ the pair computes a linear function, like $\mathrm{ReLU}(u)-\mathrm{ReLU}(-u)$ in Section 2.

All six numbers are floating-point computations. In exact arithmetic any feasible $(\rho,T)$ is an upper bound and enumerating the full-dimensional patterns gives $L^\star$ exactly; here the barrier stops at a finite tolerance and the Cholesky test and the pattern LPs run with rounding and thresholds, so the labels say "numerical" and "estimate". A rigorous certificate would need validated arithmetic (background box in Section 6). If a Cholesky check fails, the explorer says so instead of labelling the number a bound.

What to look for. (1) With $\kappa=0$ LipSDP-Neuron equals the all-patterns bound for $280$ of the first $300$ seeds, and the largest relaxation gap among the rest is $3.4\%$. The product bound is typically $24$–$86\%$ above $L^\star$ (quartiles over the same seeds, see Primer C; median $45\%$): LipSDP removes the alignment cost. (2) "Walkthrough network" loads the example above: product $4$, LipSDP $2.309$, $L^\star=\sqrt5\approx2.236$. The small gap between LipSDP and the all-patterns bound is the relaxation gap of the S-procedure step, and LipSDP-Layer equals LipSDP-Neuron by the symmetry argument of Step 6, to the digits shown. The explorer pads the network with a third neuron that has output weight $0$ and input row $(0.01,0)$. LipSDP-Neuron gives it multiplier $0$, but the shared Layer multiplier cannot, which raises LipSDP-Layer by about $5\cdot10^{-6}$. (3) Move $\kappa$ to $1$. At exactly $\kappa=1$ neurons 1 and 3 can no longer be active together; for $80\%$ of the seeds $L^\star$ drops by more than $10\%$ between $\kappa=0.99$ and $\kappa=1$ (median drop $34\%$), while the all-patterns bound and both LipSDP bounds, which depend continuously on the weights, change much less: the median change is below $0.7\%$, all three stay within $1\%$ for $191$ of the $300$ seeds, and the largest change is $6.5\%$. This is the impossible-pattern cost that no bias-independent certificate can avoid (Section 4). For $\kappa$ slightly below $1$ the both-active region still exists as a thin wedge. The enumerated value counts it, and sampling often misses it: at $\kappa=0.99$ the sampled value reaches $L^\star$ for only $64\%$ of the seeds. (4) Lower the scale $s$. The product, both LipSDP bounds, the all-patterns bound and $L^\star$ scale like $s^2$, so their ratios do not change: the biases are not scaled, so $f_s(x)=s\,f_1(sx)$ and $\nabla f_s(x)=s^2\nabla f_1(sx)$. The sampled value does not scale like $s^2$, because the square stays fixed while the boundaries move away from the origin relative to it, so sampling sees fewer regions: at $s=0.4$ it reaches $L^\star$ for $68\%$ of the seeds, against $85\%$ at $s=1$.

Why can $L^\star$ jump at $\kappa=1$ although $f$ changes continuously at every input? $L^\star$ is a supremum over all inputs, including distant ones. For $\kappa\lt1$ the steep both-active region can be a thin wedge, or lie far away, and still determine that supremum; at $\kappa=1$ its two boundaries coincide, the region loses its interior and its Jacobian no longer counts. Sampling in one fixed square can miss the region long before it disappears.

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