1 Introduction↩︎

Deep learning has achieved remarkable success across various domains, including computer vision, natural language processing, and graph-based learning. Advances in state-of-the-art architectures, such as convolutional neural networks (CNNs) [1], transformers [2], [3], and graph neural networks (GNNs) [4], as well as their large-scale applications, including large language models (LLMs) [3], [5]–[8] and large vision-language models (LVLMs) [9]–[11], have significantly expanded the frontiers of deep learning, enabling powerful capabilities across a wide range of tasks.

Beyond their remarkable success, ensuring safety [12], robustness [13] and generalization, remains a critical challenge. These properties are intrinsically linked to the sensitivity of neural networks to input perturbations [13]–[16], which can significantly impact their reliability in real-world applications. For instance, adversarial vulnerabilities [17], [18], and out-of-distribution (OOD) generalization [19]–[23] are major concerns in developing reliable deep learning models. A fundamental mathematical concept that characterizes network sensitivity and governs these properties is Lipschitz continuity, which measures the upper bound of the outputs of neural networks to inputs.

Although numerous studies have explored Lipschitz continuity from both theoretical [21], [24] and applied [25]–[27] perspectives, a unified survey that consolidates these insights is still lacking. Existing works primarily focus on isolated aspects, such as adversarial robustness [13], [28], [29], regularization [25], [26], or theoretical bounds without practical considerations [21]. The absence of a comprehensive survey integrating these perspectives leaves a critical gap in the literature. To address this, we provide a systematic review of Lipschitz continuity in deep learning, covering its theoretical foundations, estimation methods, regularization approaches in optimization and architecture, and certifiable robustness. By synthesizing key findings across these domains, this survey serves as a valuable resource for researchers and practitioners seeking a deeper understanding of Lipschitz continuity and its role in trustworthy learning systems. Beyond conducting a survey, we also critically distill the existing knowledge, present it in a more accessible manner, and provide the results or proofs missing or incorrect in the literature. The message of this survey is that Lipschitz continuity is not just a niche mathematical property; it is a fundamental principle for building and analyzing trustworthy neural networks.

1.1 Scope↩︎

This survey aims to consolidate existing knowledge across these domains, providing a comprehensive resource for researchers and practitioners seeking to understand and leverage Lipschitz continuity in deep learning. To ensure a high-quality review, we focus primarily on the literature published in prestigious and other well-regarded sources. The scope of this survey is organized as below:

1. Section 2:  presents the formal definitions and fundamental properties of Lipschitz continuity, followed by key theoretical developments in deep learning. It covers the mathematical foundations, Lipschitz properties of activation functions and neural networks, Lipschitz analysis of multi-head self-attention, complexity-theoretic generalization bounds, and advances in training dynamics.

2. Section 3:  surveys the principal methods for estimating the Lipschitz constants of neural networks, including derivative bound propagation, spectral alignment, convex optimization relaxations, and relaxed integer programming approaches.

3. Section 4:  reviews approaches for enforcing Lipschitz continuity in neural networks through both optimization-based and architectural approaches, covering weight regularization, gradient regularization, activation regularization, and class-margin regularization.

4. Section 5:  examines methods for achieving certifiable robustness through Lipschitz bounds, including the theoretical foundations of Lipschitz-margin robustness radii, certification via global Lipschitz bounds, certification via local Lipschitz bounds and certificated robustness in LLMs.

1.2 Contributions↩︎

Our key contributions are as follows:

1. A Unified Systematic Review. We close a critical gap in the literature by conducting a comprehensive survey of Lipschitz continuity in deep learning. Our review unites disparate research threads — spanning theoretical foundations, estimation methods, regularization approaches, architectural designs, and certified robustness — into a coherent framework that highlights their interconnections from the lens of Lipschitz continuity.

2. Comprehensive Analysis. We provide in-depth examinations of existing techniques organized into four categories:

   • Theoretical Foundations. We review the formal definitions and properties of Lipschitz continuity, analyze Lipschitz behavior of activations and network architectures, and synthesize key results on generalization bounds, optimization convergence, and training-induced Lipschitz dynamics.

   • Estimation Methods. From power-iteration spectral-norm approximations and extreme-value theory to randomized and interval-analysis approaches, we compare each method’s tightness, computational overhead, and scalability to modern architectures.

   • Regularization Approaches. We critically assess gradient-based penalties, spectral normalization, orthogonal and Parseval constraints, and novel Lipschitz-bounded activations, examining their theoretical guarantees, implementation trade-offs, and empirical effects on robustness, generalization, and training dynamics.

   • Certifiable Robustness. We survey Lipschitz-based certification frameworks — covering global and local Lipschitz bounds, path-based and convex-relaxation certificates, randomized smoothing, and Lipschitz-constrained architectures — analyzing their provable guarantees, computational requirements, and applicability to real-world deep models.

3. Results and Proofs. We contribute new and complementary results, including corrected Lipschitz constants for common activation functions (validated by our numerical experiments) and a Lipschitz continuity bound for arbitrary neural networks represented as directed acyclic graphs (DAGs). Furthermore, we provide theoretical results missing or incorrect from prior work, such as rigorous proofs of the closed-form expressions for the Lipschitz constants of common activation functions and a general perturbation bound for certifiable robustness based on Hölder’s conjugacy relation.

2 Theoretical Foundations↩︎

Lipschitz continuity is a foundational concept in deep learning theory, offering a rigorous framework for quantifying a neural network’s sensitivity to input perturbations. It plays a central role in understanding robustness, generalization, and optimization dynamics. Throughout, we use operator \(\mathrm{Lip}_p\left[ \bullet \right]\) to denote the \(p\)-norm Lipschitz constant for \(\bullet\), and use operator \(\mathrm{Lip}\left[ \bullet \right]\) to denote the \(2\)-norm Lipschitz constant for \(\bullet\).

2.1 Preliminaries↩︎

2.2 \(1\)-Lipschitz Orthogonal Matrix↩︎

2.2.1 Matrix Orthogonalization↩︎

We also survey useful matrix orthogonalization algorithms below.

2.2.1.1 Lie Group Exponential Map.

Let: \[\begin{align} U(n) = \Bigl\{ A \in \mathbb{C}^{n \times n} \mid A^\ast A = I \Bigr\} \end{align}\] be a unitary group and: \[\begin{align} \mathfrak{u}(n)=\{B\in\mathbb{R}^{n\times n}\;:\;B=-B^\ast \} \end{align}\] be a Skew-Hermitian group. Then there exists an Lie Exponential Map: \[\begin{align} \exp(A): \mathfrak{u}(n) \to U(n) , \end{align}\] connecting the groups \(\mathfrak{u}(n)\) and \(U(n)\), defined as: \[\begin{align} \exp(A)=I+A+\frac{1}{2}A^2+\cdots \label{eq:Lie95Exponential95Map} , \end{align}\tag{1}\] which is surjective. For a matrix \(W\), its corresponding skew-Hermitian matrix can be constructed by taking: \[\begin{align} B=W-W^\ast \in \mathfrak{u}(n) \label{eq:skew95Hermitian95conversion} \end{align}\tag{2}\] [32].

2.2.1.2 Cayley Transform.

If \(B \in \mathbb{R}^{n \times n}\) is a skew-symmetric matrix, the map: \[\begin{align} \phi(B)=(I+\frac{1}{2}B)(I-\frac{1}{2}B)^{-1} \end{align}\] is referred to as the Cayley map which transforms \(B\) to an orthogonal matrix \(\phi(B)\) [33], [34]. Then the transformed orthogonal matrix \(\phi(B)\) has a constant Lipschitz one. The Cayley map is not surjective onto the orthogonal group: its image excludes orthogonal matrices with \(-1\) in their spectrum.

2.2.1.3 Björck Orthogonalization Approximation.

Suppose \(W_0 \in \mathbb{R}^{n \times n}\) is not an orthogonal matrix. Let \(W_k\) be an iterative sequence such that \(W_k\) is the closest orthogonal matrix to \(W_0\) as \(k \to \infty\). Björck Orthogonalization [35] is a differentiable method to iteratively solve the closest orthogonal matrix to a given matrix \(W_0\). [36] use Björck Orthogonalization in an optimization problem for constraining the orthogonality of parameter matrices [36]. Björck Orthogonalization [35] is defined by: \[\begin{align} W_{k+1} = W_k \Bigl( I + \frac{1}{2} Q_k + \cdots + (-1)^p \begin{pmatrix} -\frac{1}{2} \\ p \end{pmatrix} Q_k^p \Bigr) , \end{align}\] where \(W_{k}\) is the \(k\)-th iterative result, \(p\) is a chosen hyper-parameter, \(\begin{pmatrix} -\frac{1}{2}\\ p \end{pmatrix}\) is binomial coefficient, and \(Q_k\) is: \[\begin{align} Q_k = I - W_k^\top W_k . \end{align}\] As \(k \to \infty\), \(W_k^\top W_k\) converges to \(I_n\). For a rectangular matrix \(W_0 \in \mathbb{R}^{m \times n}\), if \(\rank(W_0) = n\), the convergence still holds: \[\begin{align} W_k^\top W_k \to I_n \notag , \end{align}\] as \(k \to \infty\).

2.3 Functional Inequalities↩︎

2.4 Sub-Differential Convex Functions↩︎

For a function that is not differentiable but is locally Lipschitz continuous — such as the ReLU activation — the Clarke sub-differential provides a generalized notion of gradient that allows one to analyze and bound local Lipschitz constants [37].

2.5 Lipschitz Continuity of Activation Functions↩︎

An activation function: \[\begin{align} \rho: \mathbb{R}^d \to \mathbb{R}^d \notag \end{align}\] is a nonlinear mapping applied after neuron’s output at a layer, enabling neural networks to learn complex patterns. For example, if \(\sigma\) is piecewise defined as: \[\begin{align} \rho(x) = \begin{cases} x, & \text{if } x > 0, \\ 0, & \text{otherwise} \end{cases} , \end{align}\] then the \(\rho\) is referred to as ReLU activation function [41].

The Lipschitz constants of commonly used activation functions are provided in Table 1. The constants are derived as the supremum of the \(2\)-norm of the derivatives of activation functions. Proofs are provided in the Appendix 7. To numerically validate the results, we have conducted experiments by maximizing the gradient norm using Lemma 1 for a further sanity check.

Notes. The Lipschitz constant for the softmax function reported in the literature [42] and Sigmoid function reported in the literature [43] as 1. We show in Appendix 7.6 that the exact Lipschitz constant of softmax is indeed \(\frac{1}{2}\), validated by our numerical experiment. This result has been reported in recent concurrent literature [44]. We also show in Appendix 7.1 that exact Lipschitz constant of sigmoid is indeed \(\frac{1}{4}\), validated by our numerical experiment. For Leaky ReLU and ELU, the parameter \(\alpha\) is typically chosen with \(\alpha < 1\), so their effective Lipschitz constants are usually \(1\) in practice.

Auto-Differentiation based Numerical Method. Suppose \(f: \mathbb{R}^{d} \to \mathbb{R}\) is a first- and second-order differentiable neural network under domain convex condition (Remark 2). Leveraging the autograd mechanism of PyTorch, we use gradient descent — e.g., Adam [45] — to numerically approximate the supremum: \[\begin{align} \mathrm{Lip}[f] = \sup_x \|f'(x)\|_2 \notag , \end{align}\] using Lemma 1.

2.6 Lipschitz Continuity of Dot-Product Self-Attention↩︎

LLMs [3], [5]–[8], built on transformer architectures [2], rely heavily on dot-product self-attention mechanisms, whose Lipschitz continuity governs their sensitivity to input perturbations, such as token-level adversarial attacks. Dot-product self-attention mechanisms form the core of transformer architectures, enabling effective modeling of long-range dependencies in sequences for applications such as natural language processing and vision-language tasks.

2.6.0.1 Single-Head Dot-Product Self-Attention.

For an input sequence matrix \(x \in \mathbb{R}^{n \times d}\), where \(n\) is the sequence length and \(d\) is the input dimension, a self-attention layer computes the output as \[\begin{align} \text{Attention}(x) = \underbrace{ \mathrm{softmax}\left( \frac{Q K^{\top}}{\sqrt{d_K}} \right) }_{\text{row-wise softmax}} V \in \mathbb{R}^{n \times m} ,\label{eq:shdp95attention} \end{align}\tag{3}\] [2] where \[\begin{align} Q = x W_Q \in \mathbb{R}^{n \times m}, \quad K = x W_K \in \mathbb{R}^{n \times m}, \quad V = x W_V \in \mathbb{R}^{n \times m} \label{eq:KQV95def} , \end{align}\tag{4}\] are the query, key, and value matrices, with weight matrices: \[\begin{align} W_Q, W_K, W_V \in \mathbb{R}^{d \times d_K}\notag , \end{align}\] and \(d_K\) denotes the key (and query) vector dimension. The softmax is applied row-wise to normalize attention scores. The scaling factor \(\frac{1}{\sqrt{d_K}}\) stabilizes the computation numerically.

2.6.0.2 Multi-Head Dot-Product Self-Attention.

In practice, transformers employ \(h\) parallel self-attention heads to jointly attend to information from different representation subspaces [2]. Each head \(i\) has its own projection matrices: \[\begin{align} W_{Q,i}, W_{K,i}, W_{V,i} \in \mathbb{R}^{d \times m_i}\notag , \end{align}\] where \(m_i = \frac{m}{h}\), \(h\) denotes the number of heads, and produces an output for head \(i\): \[\begin{align} \mathrm{Attention}_i(x) \in \mathbb{R}^{n \times m_i} . \end{align}\] The outputs are concatenated: \[\begin{align} \mathrm{Attention}(x) =\bigl[ \mathrm{Attention}_1(x),\mathrm{Attention}_2(x), \cdots, \mathrm{Attention}_h(x) \bigr] \in \mathbb{R}^{n \times m} , \end{align}\] produces the final output.

2.6.0.3 Locally Lipschitz Continuous.

Standard dot-product self-attention is locally Lipschitz continuous, and several works have derived explicit upper bounds for its local Lipschitz constant. Let \(\mathcal{B}(x_0,\delta)=\{x \mid \|x-x_0\|_2 < \delta \}\) denote a ball centered \(x_0\) with a radius \(\delta\). [52] show that for the dot-product self-attention layer within a ball \(\mathcal{B}_2(x,\delta)\) is locally Lipschitz continuous bounded by: \[\begin{align} \mathrm{Lip}\Bigl[\mathrm{Attention}; \mathcal{B}_2(x,\delta)\Bigr] \leq n(n+1)\Bigl(\|x\|_2 + \delta\Bigr)^2 \Bigl[ \|W_V\|_2 \|W_Q\|_2 \|W_K^\top\|_2 + \|W_V\|_2 \Bigr] , \end{align}\] where \(x \in \mathbb{R}^{n \times d}\), \(n\) is the sequence length, \(d\) is the input dimension, \(W_K, W_Q, W_V \in \mathbb{R}^{d \times d_K}\) are the key, query and value parameter matrices, and \(\delta\) is a \(2\)-norm radius [52]. [53] provide a tighter bound by \(\sqrt{n}\) up to a constant factor, showing that within a ball \(\mathcal{B}_2(0, \delta)\), the Lipschitz constant of dot-product self-attention is bounded by: \[\begin{align} \mathrm{Lip}\Bigl[\mathrm{Attention}; \mathcal{B}_2(0,\delta)\Bigr] \leq \sqrt{3} \|W_V\|_2 \Biggl[ \Bigl\|\frac{W_K^\top\,W_Q}{\sqrt{d_K}}\Bigr\|_2^2 \delta^4 (4n+1)+n \Biggr]^{\frac{1}{2}} , \end{align}\] where \(A\) is the attention score matrix [53].

Most recently, [54] refine the bound by explicitly incorporating the Jacobian of the softmax function, given by: \[\begin{align} \mathcal{M}(P_x) := \nabla \, \mathrm{softmax}\left(x\right) = \mathrm{diag}(P_x) - P_x P_x^\top , \end{align}\] where: \[\begin{align} P_x = \frac{\mathrm{e}^{x_i}}{\sum_j \mathrm{e}^{x_j}} \notag , \end{align}\] showing that the Lipschitz constant for a single head is bounded by: \[\begin{align} \mathrm{Lip}\Bigl[\mathrm{Attention}\Bigr] \leq \|W_V\|_2 \Biggl( \|P\|_2 + 2 \|x\|_2^2 \; \|A\|_2 \; \max_i \|\mathcal{M}(P_{i, :})\|_2 \Biggr) \end{align}\] with: \[\begin{align} P = \mathrm{softmax}\Bigl( x A x^\top \Bigr) \in \mathbb{R}^{n \times n} \qquad \text{and} \qquad A = \frac{W_QW_K^\top}{\sqrt{d_K}} \in \mathbb{R}^{d \times d} \notag \end{align}\] [54].

2.6.0.4 Globally Non-Lipschitz Continuous.

It is worth noting that the Lipschitz constants of both single-head and multi-head dot-product self-attention are usually not globally bounded for arbitrary inputs. This is because the Lipschitz bound for dot-product self-attention grows with the sequence length and input norm itself, implying that the dot-product self-attention is not globally Lipschitz continuous [55].

2.7 Lipschitz Continuity of Neural Networks↩︎

2.7.1 Lipschitz Continuity for a Directed Acyclic Graph (DAG)↩︎

To the best of our knowledge, the existing literature does not present an explicit Lipschitz bound for a general feedforward network with skip connections. Prior literature such as [43], [57] focuses strictly on sequential architectures, while norm-based capacity measures such as the path norm [58] implicitly involve path enumeration but are not formulated as direct Lipschitz bounds. More recent path-metric bounds [59] apply to arbitrary DAGs but are expressed in a rescaling-invariant metric form rather than the simple sum-over-paths structure we consider.

Figure 1 provides an illustrative example of a DAG network. To complement the survey, we derive an explicit bound in Theorem 11 using graph-theoretic method. To the best of our knowledge, this result has not previously appeared in the literature.

2.7.2 Lipschitz Constant of Non-Biconnected DAG↩︎

If a DAG is separated into two disconnected sub-DAGs by the removal of a vertex, then the removed vertex is referred to as cut vertex or articulation point. Modern deep learning architectures often give rise to computation DAGs that are not biconnected: the DAGs can be separated into multiple subgraphs by removing a small set of vertices.

2.7.3 Lipschitz Constant of Residual Network↩︎

Residual networks [60] address the optimization challenges of very deep neural architectures by introducing residual modules, which are illustrated in Figure 3. Each module \(m_{\mathrm{res}}\) has the structure: \[\begin{align} m_{\mathrm{res}}(x) = x + \phi(x), \notag \end{align}\] which consists of a non-identity unit \(\phi\) and an identity skip connection. The resulting residual mapping \(m_{\mathrm{res}}(x)\) maintains stable gradient flow, mitigates vanishing-gradient effects, and enables the training of substantially deeper networks without degradation. This simple architectural principle has proven highly effective and forms the backbone of many state-of-the-art deep learning models. Using Theorem 11 (), it is not difficult to show that the Lipschitz constant bound of this structure is: \[\begin{align} \mathrm{Lip}\left[m_{\mathrm{res}}\right] \leq \mathrm{Lip}\left[s\right] +\mathrm{Lip}\left[\phi\right] = 1 + \mathrm{Lip}\left[\phi\right] \notag , \end{align}\] which recovers the same results in [61], [42] and Proposition 5 () by setting \(s\) to an identity map.

2.8 Complexity-Theoretic Generalization Bound↩︎

2.8.0.1 Notations.

Let \(\mathcal{H}\) be a hypothesis family and \(\mathcal{H} \ni h: X \to Y\) be a hypothesis. Let \(\ell: Y \times Y \to \mathbb{R}\) be a loss function. For each \(\ell(h(x), y)\), we can associate \(\ell\) and \(h\) with a map \(G \ni g: X \times Y \to \mathbb{R}\), where \(G\) is the family of loss functions associated with the hypothesis family \(\mathcal{H}\).

Lipschitz constants can enter complexity–theoretic generalization bounds and often interact with capacity, robustness, and smoothness [62]. To quantify this capacity for a hypothesis, the concept of generalization gap is often used to measure the error between true risk and empirical risk (Definition 8). This gap for a model \(h \in \mathcal{H}\) is proportional to its empirical Rademacher complexity of the hypothesis family \(\mathcal{H}\), and the empirical Rademacher complexity is proportional to the Lipschitz constant \(\mathrm{Lip}\left[\mathcal{H}\right]:=\sup_{h \in \mathcal{H}} \mathrm{Lip}\left[h\right]\) of the hypothesis family: \[\begin{align} \text{Generalization gap of h} \sim \text{Rademacher complexity} \sim O(\mathrm{Lip}\left[\mathcal{H}\right]) \notag \end{align}\] ([62, Sec. 3]; [63, Pt. IV]). A smaller Lipschitz constant enforces smoother mappings, reducing the complexity of the function class and mitigating overfitting, thereby enhancing robustness to noise and adversarial perturbations [42], [53], [58].

2.8.1 Rademacher Complexity and its Bounds↩︎

In learning theory, Rademacher complexity quantifies the expressiveness of a hypothesis family \(G:=\{g: X \to \mathbb{R}\}\) by measuring how well the \(G\) can fit data associated with random labels [62]. Let \(S = (x_1, x_2, \dots, x_m)\) be a dataset of size \(m\) with each \(x_i \in X\) drawn i.i.d. from an unknown distribution \(\mathcal{D}_X\). Let \(h : X \to \mathbb{R}\) be a hypothesis that assigns a real-valued score to each input. Let \(\sigma_1, \dots, \sigma_m\) be \(m\) independent Rademacher variables, i.e., random variables uniformly distributed over \(\{-1, +1\}\) [62]. For random labels \(\sigma_1, \dots, \sigma_m\) assigned to \(S\), the maximal correlation between the labels and the predictions of hypotheses in \(\mathcal{H}\) is given by \[\begin{align} \sup_{g \in G} \frac{1}{m} \sum_{i=1}^m \sigma_i g(x_i). \notag \end{align}\] By assigning all possible random labels to the dataset \(S\), we can thus define Rademacher Complexity for hypothesis family \(\mathcal{G}\), stated in Definition 9 ().

In the paradigm of machine learning, for the case of the hypothesis family \(\ell\circ \mathcal{H}\) is the composition of a model \(h \in \mathcal{H}:=\{h:\mathbb{R}^n \to \mathbb{R} \}\) and loss function \(\ell: \mathbb{R} \to \mathbb{R}\): \[\begin{align} \ell\circ \mathcal{H} := \{\ell\circ h \mid h \in \mathcal{H}\}, \notag \end{align}\] by Talagrand’s Contraction Lemma [62], the empirical Rademacher complexity of \(\ell\circ \mathcal{H}\) on \(S\) admits: \[\begin{align} \widehat{\mathfrak{R}}_S(\ell\circ \mathcal{H}) \leq \mathrm{Lip}\left[\ell\right] \;\widehat{\mathfrak{R}}_S(\mathcal{H}) . \notag \end{align}\]

For any \(\delta > 0\), with probability at least \(1-\delta\), the generalization gap for \(h \in \mathcal{H}\) on \(S\) is bounded by empirical Rademacher complexity of \(\ell\circ \mathcal{H}\): \[\begin{align} R(h) - R_S(h) \leq 2 \;\widehat{\mathfrak{R}}_S(\ell\circ\mathcal{H}) + 3 \sqrt{\frac{\log{\frac{2}{\delta}}}{2m}} \leq 2\; \mathrm{Lip}\left[\ell\right] \;\widehat{\mathfrak{R}}_S(\mathcal{H}) + 3 \sqrt{\frac{\log{\frac{2}{\delta}}}{2m}} \end{align}\] [62].

For general vectorized function, let: \[\begin{align} \Phi:=\{\phi: \mathbb{R}^n \to \mathbb{R}^n \mid \phi~\text{is a coordinate-wise identity map}\} \notag \end{align}\] be a coordinate-wise identity family. For vectorized function \(h: \mathbb{R}^n \to \mathbb{R}^k\) and compositional hypothesis \(h \circ \Phi\), suppose \(h_i\) is \(K\)-Lipschitz continuous, using vector contraction theorem [64] yields: \[\begin{align} \widehat{\mathfrak{R}}_S(h \circ \Phi) &\leq \sqrt{2} \; K \; \frac{1}{m} \; \mathbb{E}\Bigl[ \sup_{\phi \in \Phi} \sum_{i=1}^m \sum_{j=1}^n \sigma_{ij} \phi_j(x_i)\Bigr] \notag\\ &\leq \sqrt{2} \; K \; \frac{1}{m} \; \sum_{i=1}^m \sum_{j=1}^n \bigl|\phi_j(x_i)\bigr| = \sqrt{2} \; K \; \frac{1}{m} \; \sum_{i=1}^m \Biggl(\sum_{j=1}^n \bigl|x_{i:j}\bigr| \Biggr) = \sqrt{2} \; K \; \frac{1}{m} \; \sum_{i=1}^m \bigl\|x_i\bigr\|_1 \notag \\ &\leq \sqrt{2} \; K \; \frac{1}{m} \; m \; \sup_{x \in S} \bigl\| x \bigr\|_1 = \sqrt{2} \; K \; \sup_{x \in S} \bigl\| x \bigr\|_1 , \end{align}\] where \(\sigma_{ij}\) are an \(m \times n\) matrix of independent Rademacher variables, \(x_{i:j}\) denotes the \(j\)-th component of the \(i\)-th data point, and \(\sup_{x \in S} \bigl\| x \bigr\|_1\) is the \(1\)-norm dataset diameter. This result implies that Lipschitz constant of \(h\) controls its generalization bound: \[\begin{align} R(h) - R_S(h) \leq \mathcal{O}(K) . \end{align}\]

2.9 Training Dynamics of Lipschitz Continuity↩︎

The Lipschitz continuity of a parameterized network evolves as optimization updates its parameters. Consequently, the optimization process induces the dynamics of the network’s Lipschitz continuity. While existing work in deep learning theory has largely focused on bounding the Lipschitz constants of neural networks, understanding how Lipschitz continuity evolves throughout training remains an open problem. Only a few studies have begun to explore this aspect. In particular, [24] introduces the concept of optimization-induced dynamics and establishes a continuous-time stochastic framework that explicitly links the optimization process to the time-varying Lipschitz continuity bound of deep networks. The literature [24] uses the operator-theoretical results regarding how the largest singular values of parameter matrices vary with respect to perturbations from the literature [65], along with the theory from high-dimensional stochastic differential equations (SDEs) from stochastic analysis, for establishing a mathematical framework modeling the evolution of a spectral-norm based Lipschitz continuity upper bound. They further validate their theoretical framework with experiments across datasets and regularization scenarios.

Consider an \(L\)-layer feed-forward neural network \(f: \mathbb{R}^m \to \mathbb{R}\) with 1-Lipschitz activation functions (e.g., ReLU). Let \(\boldsymbol{\theta}^{(\ell)}(t) \in \mathbb{R}^{m_\ell \times n_\ell}\) be the \(\ell\)-layer parameter matrix at time \(t\) and \(\boldsymbol{\theta}(t)\) be the collection of \(\boldsymbol{\theta}^{(\ell)}(t)\) for all layers. Let \(\mathcal{L}_f(\boldsymbol{\theta}^{(\ell)}(t))\) be the loss expectation for the \(f\) with parameters \(\boldsymbol{\theta}(t)\) at time \(t\). Suppose that \(\boldsymbol{\Sigma}_t^{(\ell)} \in \mathbb{R}^{(m_\ell n_\ell) \times (m_\ell n_\ell)}\) is the layer-wise covariance matrix of stochastic gradient noise, arising from mini-batch sampling. Let \(K^{(\ell)}(t)\) be the Lipschitz constant at layer \(\ell\). Let \(K(t)\) be the network Lipschitz spectral-norm bound. [24] show that the continuous-time dynamics of the parameters under stochastic gradient descent (SGD) with a sufficiently small learning rate \(\eta > 0\) are given by a system of SDEs [24]:

\[\require{physics} \begin{align} \begin{cases} \dd \mathrm{vec}(\boldsymbol{\theta}^{(\ell)}(t)) = - \mathrm{vec}\left[\nabla^{(\ell)} \mathcal{L}_f(\boldsymbol{\theta}(t))\right] \, \dd t + \sqrt{\eta} \left[\boldsymbol{\Sigma}_t^{(\ell)} \right]^{\frac{1}{2}} \, \dd \boldsymbol{B}_t^{(\ell)} \\ K^{(\ell)}(t)=\|\boldsymbol{\theta}^{(\ell)}(t)\|_{op} \\ Z(t) = \sum_{l=1}^L \log K^{(\ell)}(t) \\ K(t) = \mathrm{e}^{Z(t)} \end{cases} \end{align}\] where:

• \(\mathrm{vec}: \mathbb{R}^{m \times n} \to \mathbb{R}^{mn}\) represents major-column vectorization operator.

• \(\|\cdot\|_{op}\) the matrix operator (spectral) norm.

• \(\nabla^{(\ell)} \mathcal{L}_f(\boldsymbol{\theta}(t))\) is the gradient of the loss with respect to the \(\ell\)-th layer parameters.

• \(\boldsymbol{B}_t^{(\ell)}\) is a standard \((m_\ell n_\ell)\)-dimensional Wiener process adapted to the filtration induced by mini-batch sampling.

This stochastic dynamical system characterizes, in continuous time, the evolution of both layer-wise and network-level Lipschitz continuity bounds induced by the optimization process. [24] show that the layer-wise dynamics decompose into three driving forces: \[\require{physics} \begin{align} \frac{\dd K^{(\ell)}(t)}{K^{(\ell)}(t)} = \left( \mu^{(\ell)}(t) + \kappa^{(\ell)}(t) \right) \, \dd t + \boldsymbol{\lambda}^{(\ell)}(t)^\top \, \dd \boldsymbol{B}_t^{(\ell)} , \end{align}\] where \(\require{physics} \dd\,\boldsymbol{B}_t^{(\ell)} \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{\mathrm{I}}_{m_{\ell}n_{\ell}} \dd t)\) represents the increment of a standard Wiener process in \(\mathbb{R}^{m_{\ell}n_{\ell}}\), and the three driving forces are:

1. Optimization-induced drift: \[\begin{align} \mu^{(\ell)}(t) &= \frac{\left\langle \boldsymbol{J}_{op}^{(\ell)}(t), -\mathrm{vec}\!\left[\nabla^{(\ell)} \mathcal{L}_f(\boldsymbol{\theta}(t))\right] \right\rangle}{\sigma_{1}^{(\ell)}(t)}, \end{align}\] where \(\boldsymbol{J}_{op}^{(\ell)}(t)\) is the Jacobian of the operator norm with respect to \(\mathrm{vec}(\boldsymbol{\theta}^{(\ell)}(t))\) and \(\sigma_{1}^{(\ell)}(t)\) is its largest singular value. This term captures the deterministic component of Lipschitz evolution driven by the mean gradient flow. This shows that the alignment of the gradient \(\nabla^{(\ell)} \mathcal{L}_f(\boldsymbol{\theta}(t))\) and principal direction of the \(\ell\)-th layer parameter matrix \(\boldsymbol{J}_{op}^{(\ell)}(t)\) at time \(t\) determines the contribution of the optimization to the increment of the layer-wise Lipschitz constant.

2. Noise–curvature entropy production: \[\begin{align} \kappa^{(\ell)}(t) &= \frac{\eta}{2\,\sigma_{1}^{(\ell)}(t)} \left\langle \boldsymbol{H}_{op}^{(\ell)}(t), \boldsymbol{\Sigma}_t^{(\ell)} \right\rangle \ge 0, \end{align}\] where \(\boldsymbol{H}_{op}^{(\ell)}(t)\) is the Hessian of the operator norm and \(\boldsymbol{\Sigma}_t^{(\ell)}\) is the gradient-noise covariance. This non-negative term quantifies irreversible growth of the Lipschitz bound due to the interaction between stochastic gradient noise and curvature.

3. Diffusion-modulation intensity: \[\begin{align} \boldsymbol{\lambda}^{(\ell)}(t) &= \frac{\sqrt{\eta}}{\sigma_{1}^{(\ell)}(t)} \left[\boldsymbol{\Sigma}_t^{(\ell)}\right]^{1/2\top} \boldsymbol{J}_{op}^{(\ell)}(t), \end{align}\] which controls the variance of Lipschitz evolution by scaling the Wiener noise term from mini-batch sampling.

[24] show that, at the network level, these quantities aggregate as: \[\begin{align} \mu_Z(t) = \sum_{\ell=1}^L \mu^{(\ell)}(t) ,\qquad \kappa_Z(t) = \sum_{\ell=1}^L \kappa^{(\ell)}(t) , \qquad \lambda_Z(t) = \left[\sum_{\ell=1}^L \|\boldsymbol{\lambda}^{(\ell)}(t)\|_2^2\right]^{\frac{1}{2}} , \end{align}\] governing the deterministic trend, irreversible growth, and stochastic fluctuations of the overall Lipschitz bound \(K(t)\).

Their framework is particularly useful for interpreting the behaviors of neural networks, such as the near-convergence behavior, noisy supervision and mini-batch sampling trajectories. [24] show that:

1. The Lipschitz constant bound irreversibly increases since the term noise-curvature entropy production \(\boldsymbol{\kappa}_Z(t)\) is non-negative, which increases the system entropy [24].

2. The magnitude of uniform label noise affects the Lipschitz bound [24]. In particular, larger supervision noise leads to a lower Lipschitz bound for the network. This is because the supervision noise shrinks the optimization-induced drift \(\boldsymbol{\mu}_Z(t)\).

3. Mini-batch trajectories do not affect the variance of the Lipschitz bound if batch size is sufficiently large [24].

3 Estimation Methods↩︎

Proposition 10 gives an upper bound of the Lipschitz constant of a network [24], [27]. However, exact computation of the Lipschitz constant \[\begin{align} K = \sup_{x \neq y} \frac{\|f(x) - f(y)\|_2}{\|x - y\|_2} = \sup_{x} \|\nabla f(x)\|_2 \notag \end{align}\] is generally NP-hard [43], and even inapproximable under the Exponential Time Hypothesis [39]. This section surveys major estimation and bounding techniques, ranging from statistical sampling to certified convex relaxations. We group them into Power Iteration, Extreme Value Theory, Derivative Bound Propagation, Spectral Alignment, Convex Optimization Relaxations, and Exact/Relaxed MILP formulations.

3.1 Power Iteration for Single Linear Unit↩︎

Power Iteration [66] is also known as the Von Mises iteration [66]. It is the standard way to estimate the Lipschitz constant for a linear layer — fully-connected or convolutional, by approximating the spectral-norm of weight matrices, which bounds the layer Lipschitz constant (Proposition 8) [24], [27], [55], [57], [67]. Another desirable property of power iteration in optimization is that power iteration is differentiable. For example, [57] and [27] use power iteration for computing the largest singular values of parameter matrices, in which the differentiability of power iteration allows the gradient propagation for gradient based optimization.

Suppose that \(W \in \mathbb{R}^{m \times m}\) is the parameter matrix of a fully-connected layer or a convolutional layer, then the largest singular value \(\sigma_1(W)\) is its operator norm: \[\begin{align} \mathrm{Lip}[W] = \sigma_1(W) . \end{align}\] Computing the exact value of \(\sigma_1(W)\) for a large \(m \times m\) matrix \(W\) is computationally expensive, as classical algorithms require \(O(m^3)\) time. However, the largest singular value can be approximated using power iteration by assuming the existence of a uniquely dominant singular value [68]. Starting with random unit vectors \(u_0, v_0\), one step of power iteration updates \[\begin{align} v_{k+1} \;\leftarrow\; \frac{W^\top u_k}{\|W^\top u_k\|_2}, \qquad u_{k+1} \;\leftarrow\; \frac{W v_k}{\|W v_k\|_2} , \end{align}\] and uses \[\begin{align} \hat{\sigma}_1^{(k+1)} = u_{k+1}^\top W v_{k+1} \end{align}\] as a differentiable estimate of \(\sigma_1(W)\) at step \(k+1\). At the iteration number \(T\), the estimated Lipschitz constant of \(W\) is approximated as: \[\begin{align} \mathrm{Lip}[W] \approx u_{T}^\top W v_T , \end{align}\] where \(T\) is the number of iterations. The estimation error is proportional to: \[\begin{align} |\sigma_1 - \hat{\sigma}_1^{(T)} | = O\Bigl( \Bigl[ \frac{\sigma_2}{\sigma_1} \Bigr]^T \Bigr) , \end{align}\] where \(\sigma_1\) and \(\sigma_2\) are the largest, and second largest singular values, respectively [68].

3.2 Extreme Value Theory↩︎

CLEVER [29] — Cross-Lipschitz Extreme Value for nEtwork Robustness — converts the problem of estimating the attack-independent robustness into the problem of estimating the local Lipschitz constant by sampling gradient norms in a neighborhood of the input [29]. This is because the \(q\)-norm Lipschitz constant of a network \(g: \mathbb{R}^m \to \mathbb{R}\) is determined by: \[\begin{align} \mathrm{Lip}_q\left[g\right] = \sup_x \|\nabla\; g(x) \|_q , \end{align}\] so that, for any inputs \(x\) and \(y\), the inequality holds: \[\begin{align} |g(x) - g(y)| \leq \mathrm{Lip}_q\left[f\right] \|x-y\|_p \notag , \end{align}\] [29] where \(1 \leq p, q \leq \infty\) satisfy the Hölder conjugacy relation [31]: \[\begin{align} \frac{1}{p} + \frac{1}{q} = 1 . \notag \end{align}\]

3.2.0.1 Adversarial Perturbation Bound.

[29] further show that the local Lipschitz constant can be used for deriving the adversarial perturbation bound. Let \[\begin{align} f: \mathbb{R}^m \to \mathbb{R}^k \notag \end{align}\] be a \(k\)-class classifier and \(f_i\) be the \(i\)-th prediction. Let: \[\begin{align} c(x) = \mathop{\mathrm{arg\,max}}_{1 \leq i \leq k} f_i(x) \notag \end{align}\] be the prediction. Then the \(p\)-norm adversarial perturbation \(\|\delta\|_p\) to an input \(x_0\) is bounded above by: \[\begin{align} \|\delta\|_p \leq \min_{j \neq c} \frac{f_c(x_0) - f_j(x_0)}{\mathrm{Lip}_q\left[f_c - f_j\right]} ,\notag \end{align}\] where the local Lipschitz constant \(\mathrm{Lip}_q\left[f_c - f_j\right]\) can be estimated through: \[\begin{align} \mathrm{Lip}_q\left[f_c - f_j\right] = \max_{x \in B_p(x_0, \delta)} \|\nabla\; (f_c - f_j)(x) \|_q \end{align}\] and \(B_p(x_0, \delta)\) is a \(p\)-norm ball centered at \(x_0\) with a radius \(\delta\) [29].

3.2.0.2 Distribution of Extreme Value \(\max_{x \in B_p(x_0, \delta)} \|\nabla\; (f_c - f_j)(x) \|_q\).

Estimating: \[\begin{align} \max_{x \in B_p(x_0, \delta)} \|\nabla_x\; (f_c - f_j)(x) \|_q \end{align}\] can be through sampling \(x \in B_p(x_0, \delta)\). Suppose the \(n\) samples are \(\{x^{(1)}, x^{(2)}, \cdots, x^{(n)} \}\). Their gradient \(q\)-norms are \(\{\|\nabla\; (f_c - f_j)(x^{(1)}) \|_q, \|\nabla\; (f_c - f_j)(x^{(2)}) \|_q, \cdots, \|\nabla\; (f_c - f_j)(x^{(n)}) \|_q\}\). [29] show that the extreme value of the gradient \(q\)-norm: \[\begin{align} Y := \lim_{n \to \infty} \; \max_{x^{(i)} \in B_p(x_0, \delta)} \Bigl\{\|\nabla\; (f_c - f_j)(x^{(1)}) \|_q, \|\nabla\; (f_c - f_j)(x^{(2)}) \|_q, \cdots, \|\nabla\; (f_c - f_j)(x^{(n)}) \|_q \Bigr\} \end{align}\] can only be a distribution \(P_Y\) of three distribution classes:

• Type I: Gumbel class,

• Type II: Fréchet class,

• Type III: Reverse Weibull class,

according to Fisher-Tippett-Gnedenko Theorem [29]. The extreme value \(\max\|\nabla\; (f_c - f_j)(x^{(1)}) \|_q\) is thus converted to estimating the right end-point (the location parameter \(a_W\)) of the reverse Weibull distribution \(P_Y\): \[\begin{align} \mathrm{Lip}_q\left[f_c - f_j\right] \approx a_W . \end{align}\]

3.3 Coordinate-Wise Gradient↩︎

Fast-Lip [69] derives an upper bound of local Lipschitz constant in the form of coordinate-wise gradients by analyzing the activation patterns of ReLU networks. Let \(n_k\) denote the number of neurons at the \(k\)-th layer of an \(m\)-layer network and Let \(n_0\) be the input dimension. Let \(\phi_k: \mathbb{R}^{n_0} \to \mathbb{R}^{n_k}\) be the map from input to the output of \(k\)-th layer. Let \(\rho\) be the coordinate-wise activation function. Then the relation between the \((k-1)\)-th layer and the \(k\)-th layer can be written as: \[\begin{align} \phi_k(x) = \rho\Bigl( W^{(k)} \phi_{k-1}(x) + b^{(k)} \Bigr) , \end{align}\] where \(W^{(k)} \in \mathbb{R}^{n_k \times n_{k-1}}\) is the \(k\)-th layer parameter matrix, and \(b^{(k)} \in \mathbb{R}^{n_k}\) is the bias. Set \(f(x) = \phi_m(x)\) and let \(f_j(x)\) denote the \(j\)-th output of \(f\). [69] start from analyzing the ReLU activation patterns and then deriving gradient bounds under this setting for bounding local Lipschitz constant.

3.3.0.1 Activation Patterns of ReLU Networks.

Let \(l_r^{(k)}\) and \(u_r^{(k)}\) denote the lower and upper bound for the \(r\)-th neuron in the \(k\)-th layer and let \(z_r^{(k)}\) be the pre-activation at the \(k\)-th layer, given as: \[\begin{align} z_r^{(k)} = W_{r,:}^{(k)} \phi_{k-1}(x) + b_r^{(k)} , \end{align}\] where \(W_{r,:}^{(k)}\) denotes the \(r\)-th row of \(W^{(k)}\) [69]. For the neurons indexed by \([n_k]:=\{1, 2, \cdots, n_k\}\), then there are only three activation patterns:

1. Always Activated. Neurons are always activated: \(\mathcal{I}_k^+ := \{ r \in [n_k] \mid u_r^{(k)} \geq l_r^{(k)} \geq 0\}\).

2. Always Inactivated. Neurons are always inactivated: \(\mathcal{I}_k^- := \{ r \in [n_k] \mid l_r^{(k)} \leq u_r^{(k)} \leq 0\}\).

3. Either Activated or Inactivated. Neurons are either activated or inactivated: \(\mathcal{I}_k := \{ r \in [n_k] \mid l_r^{(k)} \leq 0 \leq u_r^{(k)} \}\).

3.3.0.2 Gradient Analysis of ReLU Networks.

Gradient norm carries the information for local Lipschitz constant. [69] then analyze the coordinate-wise gradients of ReLU networks with a manner of layer-wise. Using the activation patterns of ReLU networks, the \(k\)-th layer output can be rewritten into: \[\begin{align} \phi_k(x) = \Lambda^{(k)} \Bigl( W^{(k)}\phi_{k-1}(x) + b^{(k)} \Bigr) , \end{align}\] where \(\Lambda^{(k)}\) is the activation pattern matrix: \[\begin{align} \Lambda^{(k)}_{r,r} = \begin{cases} 1~\text{or}~0, &\text{if}~r \in \mathcal{I}_k \\ 1, &\text{if}~r \in \mathcal{I}_k^+ \\ 0, &\text{if}~r \in \mathcal{I}_k^- \end{cases} . \end{align}\] Let \(\Lambda^{(k)}_{a}\) be the diagonal activation matrix for neurons in the \(k\)-th layer that are always activated and set \(\Lambda^{(k)}_{u}: = \Lambda^{(k)} - \Lambda^{(k)}_{a}\). Starting by analyzing the bound of the gradient \(\nabla \phi_k(x)\) coordinate-wisely in a \(2\)-layer ReLU network, [69] show that an inequality for gradient \(q\)-norm holds: \[\begin{align} \max_{x \in B_p(x_0, \epsilon)} \; \Bigl| \left[ \nabla f_j(x) \right]_k \Bigr| \leq \max \Bigl( C_{j,k}^{(1)} + L_{j,k}^{(1)} ,\; C_{j,k}^{(1)} + U_{j,k}^{(1)} \Bigr) , \end{align}\] where: \[\begin{align} C_{j,k}^{(1)} = W_{j,:}^{(2)} \Lambda^{(1)}_{a} W_{:,k}^{(1)} , \quad L_{j,k}^{(1)} = \sum_{i \in \mathcal{I}_1, W_{j,i}^{(2)}W_{i,k}^{(2)} < 0} W_{j,i}^{(2)}W_{i,k}^{(2)}, \quad \text{and} \quad U_{j,k}^{(1)} = \sum_{i \in \mathcal{I}_1, W_{j,i}^{(2)}W_{i,k}^{(2)} > 0} W_{j,i}^{(2)}W_{i,k}^{(2)} . \end{align}\] Then the gradient \(q\)-norm can be bounded above by: \[\begin{align} \mathrm{Lip}_q\left[f_j; B_p(x_0, \epsilon)\right] = \max_{x \in B_p(x_0, \epsilon)} \| \nabla f_j(x) \|_q \leq \Biggl[ \sum_k \Bigl( \max_{x \in B_p(x_0,\epsilon)} \bigl| [\nabla f_j(x)]_k \bigr| \Bigr)^{q} \Biggr]^{\frac{1}{q}} , \end{align}\] yielding a local Lipschitz constant upper bound. This method is referred to as Fast-Lip in [69].

3.4 Spectral Alignment↩︎

SeqLip [43] refines the upper bound of an \(L\)-layer sequential network \(f \in \mathbb{R}^{n_1} \to \mathbb{R}^{n_L}\) with \(1\)-Lipschitz activation functions: \[\begin{align} \mathrm{Lip}\left[f\right]^+ = \prod_{\ell=1}^{L} \|W^{(\ell)}\|_2 , \end{align}\] where \(W^{(\ell)} \in \mathbb{R}^{n_{\ell} \times n_{\ell+1}}\) is the \(\ell\)-th layer parameter matrix, and \(\mathrm{Lip}\left[f\right]^+\) is referred to as AutoLip bound in [43]. By taking into account the spectral alignment between the parameter matrices of two consecutive layers, the bound is then refined as: \[\begin{align} \mathrm{Lip}\left[f\right] \leq \mathrm{Lip}\left[f\right]^+ \Biggl( \underbrace{ \prod_{\ell}^{L-1} \sqrt{ (1 - r_{\ell} - r_{\ell+1}) \max_{t^{(\ell)} \in [0,1]^{n_\ell}} \langle t^{(\ell)} \cdot v_1^{(\ell+1)} ,\; u_1^{(\ell)} \rangle^2 +r_{\ell} + r_{\ell+1} + r_{\ell}r_{\ell+1} } }_{\text{spectral alignment}} \Biggr) , \end{align}\] where \(t^{(\ell)} \in [0,1]^{n_\ell}\) represents the \(\ell\)-th layer gradient patterns from activation functions for \(n_{\ell}\) neurons — the coordinate-wise gradient of an activation function such as ReLU falls in \([0, 1]\), \(u_1^{(\ell)}\) is the first left-singular vector for \(W^{(\ell)}\), \(v_1^{(\ell+1)}\) is the first right-singular vector for \(W^{(\ell+1)}\), and \(r_{\ell}=\frac{\sigma_2^{(\ell)}}{\sigma_1^{(\ell)}}\) is the ratio of the second largest singular value to the first largest singular value for \(W^{(\ell)}\) [43].

[43] also show that, if \(u, v \in \mathbb{R}^n\) are two independent random vectors taken uniformly from \(\mathbb{S}^{n-1}=\{ x \in \mathbb{R}^n \mid \|x\|_2=1\}\), the following limit: \[\begin{align} \lim_{n \to \infty}\; \max_{t \in [0,1]^n} \Bigl | \langle t \cdot u \;, v \rangle \Bigr | =\frac{1}{\pi} \end{align}\] holds almost surely. Under this independent assumption, the bound reduces to: \[\begin{align} \mathrm{Lip}\left[f\right] \approx \frac{\mathrm{Lip}\left[f\right]^+}{\pi^{L-1}} \end{align}\] [43]. In real setting, the singular vectors are not independent across layers.

3.5 Convex Optimization Relaxation↩︎

LipSDP [70] interprets activation functions as gradients of convex potential functions, satisfying certain properties described by quadratic constraints. Therefore, Lipschitz constant estimation problem is treated as a semi-definite program (SDP) problem.

Let \(f: \mathbb{R}^{n_0} \to \mathbb{R}^{n_L}\) be an \(L\)-layer feed-forward network, recursively defined as: \[\begin{align} \begin{cases} x^{(0)} = x \\ x^{(k+1)} = \rho \Bigl( W^{(k)} x^{(k)} + b^{(k)} \Bigr) \end{cases} , \end{align}\] where \(x^{(\ell)}\) is the output at the \(\ell\)-th layer, \(W^{(\ell)} \in \mathbb{R}^{n_{\ell+1} \times n_{\ell}}\) is the parameter matrix at the \(\ell\)-th layer, \(b^{(\ell)} \in \mathbb{R}^{n_{\ell+1}}\) is the bias at the \(\ell\)-th layer, and \(\rho\) is the coordinate-wise activation function. For a single layer network: \[\begin{align} f(x) = W^{(1)} \rho \bigl( W^{(0)}x + b^{(0)} + b^{(1)} \bigr) , \end{align}\] [70] shows that the Lipschitz constant of \(f\) is given by a SDP problem: \[\begin{align} \mathrm{Lip}\left[f\right] \leq \sqrt{t} , \end{align}\] defined by: \[\begin{align} \text{minimize} \qquad &t \\ \text{subject to} \qquad &M(t, T) \preceq 0 \quad T \in T_n, \end{align}\] where \(\rho\) is slope-restricted on \([\alpha, \beta]\) [70]: \[\begin{align} \alpha \leq \frac{\rho_j(x) - \rho_j(y)}{x-y} \leq \beta \qquad \forall x,y \in \mathbb{R} , \end{align}\] \(T_n\) is a convex set [70]: \[\begin{align} T_n:= \bigl\{ T \in \mathbb{S}^n \mid T = \sum_{i=1}^n \lambda_{ii}e_ie_i^\top, \lambda_{ii} \geq 0 \bigr\} , \end{align}\] and: \[\begin{align} M(t, T) := \begin{bmatrix} -2\alpha\beta (W^{(0)})^\top T W^{(0)} - t I_{n_0} &(\alpha + \beta)(W^{(0)})^\top T \\ (\alpha + \beta) T W^{(0)} & -2T + (W^{(1)})^\top W^{(1)} \end{bmatrix} \preceq 0 \end{align}\] holds true for \(T \in T_n\) [70]. For more general \(L\)-layer fully-connected network, [70] show that \(M(t, T)\) is: \[\begin{align} M(t, T) = \begin{bmatrix} A \\ B \end{bmatrix}^\top \begin{bmatrix} -2\alpha\beta\,T & (\alpha + \beta)T \\ (\alpha + \beta)T & -2T \end{bmatrix} \begin{bmatrix} A \\ B \end{bmatrix} + \begin{bmatrix} -t I_{n_0} &0 &\cdots &0 \\ 0 &0 &\cdots &0 \\ \vdots &\vdots &\ddots &\vdots \\ 0 &0 &\cdots &(W^{(L)})^\top\,W^{(L)} \end{bmatrix} \preceq 0 , \end{align}\] where: \[\begin{align} A=\begin{bmatrix} W^{(0)} &0 &\cdots &0 &0 \\ 0 &W^{(1)} &\cdots &0 &0 \\ \vdots &\vdots &\ddots &\vdots &\vdots \\ 0 &0 &\cdots &W^{(L-1)} & 0 \end{bmatrix} , \qquad B=\begin{bmatrix} 0 &I_{n_1} &0 &\cdots &0 \\ 0 &0 &I_{n_2} &\cdots &0 \\ \vdots &\vdots &\vdots &\ddots &\vdots \\ 0 &0 &0 &\cdots &I_{n_L} \end{bmatrix} , \end{align}\] and: \[\begin{align} C=\begin{bmatrix} 0 &\cdots &0 &W^{(L)} \end{bmatrix} , \qquad b=\begin{bmatrix} (b^{(0)})^\top \cdots (b^{(L-1)})^\top \end{bmatrix}^\top \end{align}\] [70].

3.6 Integer Programming for ReLU Networks↩︎

It is well known that computing the Lipschitz constant of an arbitrary scalar- or vector-valued function is NP-hard [43], and even inapproximable under the Exponential Time Hypothesis [39]. A common approximation is to estimate the Lipschitz constant using gradient norms, where the Jacobians can be explicitly derived via the chain rule. However, nondifferentiabilities, e.g., those arising in ReLU networks, can introduce inaccuracies into these estimates. An alternative approach is to formulate the Lipschitz bounding problem as an integer programming problem, which avoids such inaccuracies by directly bounding the Jacobians.

LipMIP [39] is a method that provably exactly compute \(\ell_1\)- and \(\ell_\infty\)-norm Lipschitz constants of non-smooth ReLU networks. We summarize their method by simplifying their notations. Let \(f: \mathbb{R}^d \to \mathbb{R}\) be a ReLU network. Let \(\rho\) be activation function. Then, an \(L\)-layer ReLU network is recursively defined as: \[\begin{align} \begin{cases} f(x) = \rho(Z_L(x)) \\ Z^{(\ell)}(x) = W^{(\ell)}\;\rho\Bigl(Z^{(\ell-1)}(x)\Bigr) + b^{(\ell)} \\ Z^{(0)} = x \end{cases} , \end{align}\] where \(W^{(\ell)} \in \mathbb{R}^{n_{\ell} \times n_{\ell-1}}\) is the \(\ell\)-th layer parameter matrix, and \(Z^{(\ell)}(x)\) is the output of the \(\ell\)-th layer neurons. Let \(\partial f(x)\) be the Clarke sub-differential convex hull of \(f\) at point \(x\) — see Section 2.4 (). For example: \[\begin{align} \partial \rho(0)=[0,1] . \end{align}\]

3.6.0.1 Pushforward Sub-Differential Hull Set.

However, in practice, PyTorch implementations often set \(\partial \rho(0) = 1\), which can lead to inaccuracies in estimating the Lipschitz constant. To correctly estimate the Lipschitz constant of a ReLU network, the analysis must be carried out in the sense of sub-differentials. Let \(\partial f(X)\): \[\begin{align} \partial f(X) = \Bigl\{ g \mid\; g \in \partial f(x), x \in X \Bigr\} \end{align}\] be the set-valued Clarke sub-differentials at set \(X\). Let \(\nabla^{\sharp} f(\bullet)\) be the pushforward set that the sub-gradients of ReLU activation at value \(0\) come from the hull set \(\partial \rho(0)\): \[\begin{align} \nabla^{\sharp} f(\bullet) =\Bigl\{ \nabla f(\bullet) \mid \text{the sub-differentials of activations come from the hull set}~\partial\rho(0) \Bigr\} . \end{align}\] Thus for a ReLU network, the feasible set for its gradients on \(X\) is given as: \[\begin{align} \nabla^{\sharp} f(X) = \Bigl\{ G \in \nabla^{\sharp} f(x) \mid x \in X \Bigr\} . \end{align}\]

Theoretically, [39] show that the push forward set sub-gradient \(\nabla^{\sharp} f(x)\) is the Clarke sub-gradient hull \(\partial f(x)\), that is: \[\begin{align} \nabla^{\sharp} f(x) = \partial f(x) \end{align}\] [39] with all gradients come from \(\partial \rho(0)\). Then the local \(p \to q\) Lipschitz constant of \(f\) on \(X\) is given as: \[\begin{align} \mathrm{Lip}_{p \to q}\left[f;X\right] =\sup_{G \in \nabla^{\sharp} f(X)} \sup_{\|v\|_p \leq 1} \frac{\|G^\top\,v\|_{q}}{\|v\|_p} = \sup_{G \in \nabla^{\sharp} f(X)} \|G^\top \|_q . \end{align}\] Finding the Lipschitz constant on \(X\) can be formulated as an integer programming problem on the feasible set \(\nabla^{\sharp} f(X)\).

3.6.0.2 Solving the Supremum of Sub-Differential Hull Set.

The problem: \[\begin{align} \sup_{G \in \nabla^{\sharp} f(X)} \; \Bigl\{ \|G^\top\|_q \mid G \in \nabla^{\sharp} f(X) \Bigr\} \end{align}\] is then represented as a mixed-integer polytope, which is a mixed-integer programming problem for maximizing \(\|G^\top\|_q\) by solving the combinations of input \(x \in X\) and ReLU’s sub-gradient \(a \in [0, 1]^n\) in \(X \times [0, 1]^n\) [39]. This problem can be solved by off-the-shelf MIP solvers.

3.7 Remarks↩︎

The estimation methods surveyed above involve an inherent trade-off between computational efficiency, tightness, and scalability. Power iteration [27], [57], [66] (see Section 3.1) is lightweight and differentiable, but it only provides a layer-wise upper-bound estimate of the spectral norm, and therefore may induce a loose network-level upper bound when combined across layers (see Proposition 8 and Proposition 10). Extreme-value-theoretic approaches such as CLEVER [29] (see Section 3.2) are scalable and often yield tight empirical estimates of local Lipschitz behavior by leveraging the extreme value distributions of observed Lipschitz constants (see Lemma 1), yet they do not provide formal certificates. Methods based on coordinate-wise gradient propagation and spectral alignment, such as Fast-Lip [69] (see Section 3.3) and SeqLip [43] (see Section 3.4), improve upon spectral-norm product bounds by exploiting architectural structure (see Proposition 10), but they remain approximate and are typically tailored to specific network classes. In contrast, convex relaxation methods such as LipSDP [70] (see Section 3.5) and mixed-integer formulations such as LipMIP [39] (see Section 3.6) offer substantially stronger guarantees, with the latter being exact for ReLU networks on bounded domains by working directly with the feasible Jacobian or sub-differential set (see Section 2.4); however, these methods are considerably more computationally expensive and generally require additional architectural information to remain tractable.

4 Regularization Approaches↩︎

A variety of techniques have been developed to explicitly or implicitly enforce Lipschitz continuity in deep neural networks. These approaches can be categorized into:

1. Section 4.1: — constraining or regularizing weight matrices during initialization and training;

2. Section 4.2: — penalizing or normalizing gradient norms during training;

3. Section 4.3: — designing or constraining activation functions to be Lipschitz-bounded;

4. Section 4.4: — implicitly enforcing Lipschitz constraints via maximizing class decision boundaries;

5. Section 4.5: — enforcing Lipschitz constraints by designing architectures, particularly for transformers.

4.1 Weight Regularization↩︎

This subsection reviews methods that control the Lipschitz constant by directly constraining or regularizing the network’s weight parameters, either at initialization or during training.

4.1.1 Weight Clipping↩︎

Weight clipping is a straightforward method to enforce Lipschitz continuity by limiting the magnitude of weight matrices. This is because for a matrix \(W\), the operator-norm variation \(\|W\|_{op}\) in \(\ell_2\) for \(W\) under perturbation \(\Delta\,W\) is bounded above by: \[\begin{align} \|W\|_{op} = \sigma_1(W) = \langle u_1 v_1^\top,\; \Delta W \rangle \leq \|u_1 v_1^\top\|_2\; \|\Delta\,W\|_2 , \end{align}\] with Cauchy–Schwarz inequality, where \(\sigma_1(W)\) is the largest singular value of \(W\), \(u_1, v_1\) are the left- and right-singular vectors corresponding to the largest singular value [65].

In the setting of generative adversarial networks (GANs), [25] showthat the discriminator \(f\) of a GAN evaluates the Wasserstein distance \(W\left(\mathbb{P}, \mathbb{Q}\right)\) between two distributions \(\mathbb{P}\) and \(\mathbb{Q}\). The Kantorovich-Rubinstein duality [71] tells that: \[\begin{align} W\left(\mathbb{P}, \mathbb{Q}\right) = \frac{1}{K} \; \sup_{\mathrm{Lip}\left[f\right] \leq K} \; \mathbb{E}_{x \sim \mathbb{P}}\bigl[f(x)\bigr] - \mathbb{E}_{y \sim \mathbb{Q}}\bigl[f(y)\bigr] \end{align}\] [25]. Therefore reducing the Lipschitz constant \(K\) can reduce the Wasserstein distance \(W\left(\mathbb{P}, \mathbb{Q}\right)\). [25] propose clipping weights to a fixed range: \[\begin{align} W \leftarrow \text{clip}(W, -c, +c) \end{align}\] [25], where \(c\) is a small positive constant (e.g., \(c=0.01\)). Weight clipping is computationally efficient, but can lead to exploding or vanishing gradients if \(c\) is too small, reducing model capacity, as criticized in later literature [26].

4.1.2 Spectral Normalization and Regularization↩︎

4.1.2.1 Explicit Regularization by Penalizing Spectral Norm.

Referring to Propositions 8 () and 10 (), the product of the spectral norms of the parameter matrices gives rise to an upper bound on the Lipschitz constant of a neural network. To penalize this Lipschitz bound during training and thereby improve generalization, [57] propose, for an \(L\)-layer network parameterized by: \[\begin{align} W:=(W^{(1)},W^{(2)},\cdots,W^{(L)}) ,\notag \end{align}\] the following regularization objective: \[\begin{align} \min_{W} \mathscr{L}_{\mathrm{task}}(\xi;W) + \frac{\lambda}{2}\sum_{\ell=1}^L \big(\sigma_1^{(\ell)}\big)^2, \end{align}\] where \(\mathscr{L}_{\mathrm{task}}(\xi;W)\) represents the task loss evaluated on mini batch \(\xi\), \(\lambda\) denotes regularization coefficient and \(\sigma_1^{(\ell)}\) denotes the largest singular value of the \(\ell\)-th parameter matrix [57]. In implementation, each \(\sigma_1^{(\ell)}\) is estimated using power iteration (see Section 3.1).

4.1.2.2 Implicit Regularization by Normalizing Spectral Norm.

A challenge for Generative Adversarial Networks (GANs) [72] is that the prediction of discriminator is often inaccurate and unstable during training [27], [73]. [27] propose normalizing the weight matrix \(W\) of a linear layer by: \[\begin{align} W \leftarrow \frac{W}{\|W\|_2} = \frac{W}{\sigma_1} , \end{align}\] where \(\sigma_1\) is the largest singular value of parameter matrix \(W\), for ensuring a Lipschitz constant of at most \(1\) for the layer [27]. Readers can also refer to Proposition 9. This method was proposed by [27] for training generative adversarial networks (GANs) [72], improving the generations [27].

4.1.3 Orthogonal Weight↩︎

4.1.3.1 Parseval Networks [74].

Parseval networks enforce approximate orthonormality on weight matrices to bound their spectral-norms [74]. For a weight matrix \(W^{(\ell)}\) in layer \(\ell\), [74] minimize the \(2\)-norm of the deviation from orthonormality: \[\begin{align} \|W^{(\ell)\top} W^{(\ell)} - \mathrm{I}\|_F, \end{align}\] ensuring \(\|W^{(\ell)}\|_2 \leq 1\) [74]. For a feedforward network with \(L\) layers and 1-Lipschitz activations (e.g., ReLU), the global Lipschitz constant is bounded by: \[\begin{align} K \leq \prod_{\ell=1}^L \|W^{(\ell)}\|_2 \leq 1 , \end{align}\] (see Proposition 10 ()). [74] argue that the robustness to adversarial attacks is enhanced by the reduced sensitivity to perturbations.

4.2 Gradient Regularization↩︎

Gradient-based regularization methods enforce Lipschitz continuity by constraining the gradient norm of the loss function or network output with respect to inputs, ensuring smooth decision boundaries and robustness. Note that, for a network \(f\) on convex \(X\), its \(2\)-norm Lipschitz constant on \(X\) is given by: \[\begin{align} \mathrm{Lip}\left[f\right] = \sup_{x \in X} \|\nabla\,f(x)\|_2 , \end{align}\] (see Lemma 1 ()). Theoretically, constraining \(\|\nabla\,f(x)\|_2\) can bound Lipschitz constant of \(f\).

For example, [26] argue that bounding Lipschitz constant by weight clipping can lead to exploding or vanishing gradients, and reduce capacity of the critic in a GAN. They then introduce a regularization method by directly bounding the gradients of the critic network \(D\) through: \[\begin{align} \mathbb{E}_{x\,\sim\,\mathbb{P}_X} \Bigl[ \bigl( \|\nabla\,D(x)\|_2 -1 \bigr)^2 \Bigr] , \end{align}\] where \(\mathbb{P}_X\) is the distribution of training sample set \(X\) [26].

4.3 Activation Regularization↩︎

Activation regularization methods modulate Lipschitz continuity by constraining the properties of activation functions or their outputs. Such approaches may impose explicit norm bounds, design activations that are inherently norm-preserving, or apply penalties to activation magnitudes, thereby limiting the contribution of nonlinearities to the overall Lipschitz constant of the network.

4.3.1 Group-Sorting Activation↩︎

Bounding the Lipschitz constant of a network by ensuring \(1\)-Lipschitz for affine units through spectral-norm constraints can improve robustness [27], [57]. However, this often comes at the cost of reduced expressiveness [75]. [75] first show that, to remain expressive in a spectral-norm constrained network, the network must preserve the gradient norms at each layer [75]. ReLU networks can only satisfy this for positive values, this is because: \[\begin{align} \frac{\partial\,\mathrm{ReLU}(x)}{\partial\,x} = \begin{cases} 1, &x > 0\\ (0,1), &\text{in sub-differential sense} \\ 0, &x < 0 \end{cases} . \end{align}\] To preserve the gradient norms at activation layer, [75] use a general purpose \(1\)-Lipschitz activation function \(\boldsymbol{GroupSort}(x)\) which is homogeneous: \[\begin{align} \boldsymbol{GroupSort}(\alpha\,x)=\alpha\,\boldsymbol{GroupSort}(x) \end{align}\] [75]. GroupSort activation separates a pre-activation \(x\) into \(k\) groups: \[\begin{align} x = ( \underbrace{ x_1, \cdots, x_{g_1} }_{\text{group 1}} ,\, \underbrace{ x_{g_1+1}, \cdots, x_{g_2} }_{\text{group 2}} ,\cdots\, \underbrace{ x_{g_{k-1}+1}, \cdots, x_{g_k} }_{\text{group k}} ) , \end{align}\] and sorts each group into ascending order: \[\begin{align} \boldsymbol{GroupSort}(x) = ( \underbrace{ x_{s_1}, \cdots, x_{s_{g_1}} }_{\text{group 1}} ,\, \underbrace{ x_{s_{g_1}+1}, \cdots, x_{s_{g_2}} }_{\text{group 2}} ,\cdots\, \underbrace{ x_{s_{g_{k-1}}+1}, \cdots, x_{s_{g_k}} }_{\text{group k}} ) , \end{align}\] where the \(i\)-th sorted group satisfies: \[\begin{align} x_{s_{g_{i-1}}+1} \leq x_{s_{g_{i-1}}+2} \leq \cdots \leq x_{s_{g_{i}}} \end{align}\] [75]. For the pre-activation in a linear unit: \[\begin{align} Wx \end{align}\] where \(W \in \mathbb{R}^{m \times n}\) and \(x \in \mathbb{R}^n\), if the linear unit is with spectral-norm constraint such that: \[\begin{align} \|W\|_2 = 1 , \end{align}\] GroupSort activation directly gives rise to a \(1\)-Lipschitz constant activation: \[\begin{align} \sup_{\|x\|_2=1} \| \boldsymbol{GroupSort}(Wx) \|_2 = \sup_{\|x\|_2=1} \| Wx \|_2 = \|W\|_2. \end{align}\] When the group size is \(2\), [75] refers to this special case as MaxMin, which is equivalent to the Orthogonal Permutation Linear Unit (OPLU) [36]; when the group size is the entire input, this is referred to as FullSort in the literature [36].

To train a network with GroupSort activations and guarantee that all linear units are exactly \(1\)-Lipschitz during optimization instead of bounded by \(1\) [26], [74], [36] approximate the updated parameter matrix \(W_0\) with a closest orthogonal matrix through a differentiable, iterative algorithm, referred to as Björck Orthogonalization [35]. This algorithm is defined by: \[\begin{align} W_{k+1} = W_k \Bigl( I + \frac{1}{2} Q_k + \cdots + (-1)^p \begin{pmatrix} -\frac{1}{2} \\ p \end{pmatrix} Q_k^p \Bigr) , \end{align}\] where \(W_{k}\) is the \(k\)-th iterative result, \(p\) is a chosen hyper-parameter and \(Q_k\) is: \[\begin{align} Q_k = I - W_k^\top W_k \end{align}\] [36]. As a result, the iteration at step \(T\) leads to: \[\begin{align} W_T^\top W_T \approx I , \end{align}\] which provably ensures the linear unit with parameter matrix \(W\) is \(1\)-Lipschitz continuous: \[\begin{align} \|W_T\|_2 \approx 1. \end{align}\]

4.3.2 Contraction Activation and Invertible Residual Map↩︎

Flow-based models or flows learn to transform a source distribution \(p_X\) to target distribution \(p_Z\), consisting of invertible networks [76]. Ensuring that the transformation is invertible and its Jacobian is computable is a straightforward way for ensuring invertibility: \[\begin{align} \log\,p_X(x) = \log\,p_Z(z) + \log\, \left| \det J_F(x) \right| \end{align}\] where \(F: \mathbb{R}^m \to \mathbb{R}^m\) is an invertible map and \(J_F(x)\) is the Jacobian of \(F\) at \(x\) [77], [78]. Different from explicit computation of Jacobian, [61] propose a method for inverting residual networks [60] by ensuring the Lipschitz constants be less than \(1\), so that the residual networks are contraction maps [61], [79].

4.3.2.1 Inverting Residual Layer.

Let \[\begin{align} F(x) = x + g(x) \end{align}\] be a residual layer where \(F: \mathbb{R}^m \to \mathbb{R}^m\), \(g: \mathbb{R}^m \to \mathbb{R}^m\) and \(x \in \mathbb{R}^m\). Let \(g\) be a contraction map so that: \[\begin{align} \mathrm{Lip}\left[g\right] < 1 . \end{align}\] Suppose that the inversion \(F^{-1}\) exists and set \(y =F(x)\), so that: \[\begin{align} x = F^{-1}(y) = y - g(x). \end{align}\] A sufficient condition for inverting \(F(x)\) is that \(g\) is a contraction map: \[\begin{align} \mathrm{Lip}\left[g\right] < 1 \end{align}\] [61]. Set: \[\begin{align} T(x) = F^{-1}(y) = y - g(x) , \end{align}\] then the inversion of \(y\) is a fixed-point problem for solving: \[\begin{align} x = T(x) . \end{align}\] Using Banach contraction principle or Banach fixed-point theorem, set \(x_0=y\), the \(F^{-1}(y)\) can be approximated through: \[\begin{align} x_{k+1} = y - g(x_k) \end{align}\] iteratively.

4.3.2.2 Contraction Activation.

[80] first introduce LipSwish activation function for inverting residual networks, defined as: \[\begin{align} \boldsymbol{LipSwish}(x) = \frac{x\sigma(\beta\,x)}{1.1} , \end{align}\] where \(\sigma\) is the sigmoid function, and \(\beta > 0\) is a learnable positive constant (kept positive via softplus), initialized with \(0.5\). \(\boldsymbol{LipSwish}\) has a Lipschitz constant at most \(1\). However, the negative axis of \(\boldsymbol{LipSwish}\) has zero gradients almost everywhere. In a \(1\)-Lipschitz continuous network, the Jacobian norm across two consecutive layers is reduced to be at most one, which limits the network’s expressiveness and referred to as gradient norm attenuation problem [75], [81]. To mitigate this problem, [79] use specially designed activation function CLipSwish by concatenating two LipSwish functions. The CLipSwish is therefore given as: \[\begin{align} \Phi(x)= \begin{bmatrix} \boldsymbol{LipSwish}(x) \\ \boldsymbol{LipSwish}(-x) \end{bmatrix} , \qquad \boldsymbol{CLipSwish}(x)= \frac{\Phi(x)}{\mathrm{Lip}\left[\Phi\right]} \leq 1 , \end{align}\] where \(\mathrm{Lip}\left[\Phi\right] \approx 1.004\) [79]. This concatenation overcomes the gradient attenuation problem introduced by \(\boldsymbol{LipSwish}\) since the negative axis has zero gradients almost everywhere while ensuring the activation is a contraction map.

4.4 Class-Margin Regularization↩︎

4.4.0.1 Input Margin.

The goal of adversarial defense in classification task is to ensure that, for a \(k\)-class classifier \(f: \mathbb{R}^d \to \mathbb{R}^k\), and an input \(x \in \mathbb{R}^d\), a bounded perturbation \(\|\delta_x\|_2 < c\) does not alter the prediction: \[\begin{align} \hat{y} = \mathop{\mathrm{arg\,max}}\limits_{i \in [k]=\{1,2,\cdots,k\}} \;f_i(x+\delta_x), \end{align}\] where \(f_j(x)\) denotes the output at the \(j\)-th coordinate, and the \(\delta_x\) giving rise to the minimal \(\|\delta_x\|_2\) is often referred to as \(2\)-norm input margin [82]. More generally, input margin can also be discussed with respect to \(p\)-norm (\(1\leq p \leq \infty\)). The \(p\)-norm input margin represents the decision boundary in the input space with respect to the topology induced by \(p\)-norm.

4.4.0.2 Class/Logit Margin.

Of the same setting, the class margin (i.e. logit margin) \(m(x)\) is defined as: \[\begin{align} m(x) = f_i(x) - \max_{j \neq i} f_j(x) , \label{eq:classmargin95mx} \end{align}\tag{8}\] where \(i\) is the ground-truth [83]. By writing equation (equation 8 ) into: \[\begin{align} m(x) = (e_i - e_j)^\top \begin{bmatrix} f_i(x) \\ \max_{j \neq i} f_j(x) \end{bmatrix} , \end{align}\] where \(e_i\) and \(e_j\) are one-hot basis, for \(2\)-norm, applying Proposition 7 () immediately yields: \[\begin{align} \mathrm{Lip}\left[m\right]^2 \leq \mathrm{Lip}\left[f\right]^2 + \mathrm{Lip}\left[f\right]^2 \Longrightarrow \mathrm{Lip}\left[m\right] \leq \sqrt{2} \; \mathrm{Lip}\left[f\right] . \end{align}\] Therefore, a guaranteed \(2\)-norm perturbation \(\delta_x\) does not alter the prediction, if: \[\begin{align} \label{eq:2norm95perturbation95deltax} \|\delta_x\|_2 \leq \frac{m(x)}{\mathrm{Lip}\left[m\right]} = \frac{m(x)}{\sqrt{2} \;\mathrm{Lip}\left[f\right]} \end{align}\tag{9}\] [83], which implies that the class margin must be at least: \[\begin{align} m(x) = f_i(x) - \max_{j \neq i} f_j(x) \geq \sqrt{2} \,c\, \mathrm{Lip}\left[f\right] . \end{align}\]

4.4.0.3 Lipschitz-Margin Training.

Of the same setting, in training stage, [83] adjust each class margin \(j \neq i\) by: \[\begin{align} f_j(x) \leftarrow f_j(x) + \sqrt{2}\,c\,\mathrm{Lip}\left[f\right] , \end{align}\] where \(c\) is a guaranteed perturbation bound [83]. For each linear unit \(\phi: \mathbb{R}^m \to \mathbb{R}^n\) in the network, and let \(u \in \mathbb{R}^m\) be drawn i.i.d. from \(\mathcal{N}(0,1)\). [83] use a general method for estimating the Lipschitz constant of \(\phi\) iteratively by: \[\begin{align} u \leftarrow \frac{u}{\|u\|_2}, \qquad v \leftarrow \phi(u), \qquad \sigma \leftarrow \|v\|_2, \qquad u \leftarrow \frac{1}{2}\frac{\partial\,\|v\|_2^2}{\partial\,u} \end{align}\] [83]. At the end of the iteration \(T\), the Lipschitz constant of \(\phi\) is given by: \[\begin{align} \mathrm{Lip}\left[\phi\right] \approx \sigma , \end{align}\] almost surely.

4.5 Architectural Regularization↩︎

Lipschitz continuity can also be constrained through architectural design choices, including activation functions, parameter structures, initialization schemes, and optimization. As a complement to regularization approaches that do not fall into one single category, we survey these approaches in this section.

4.5.1 Orthogonalizing Convolution↩︎

Parameter orthogonality in neural networks can directly yield \(1\)-Lipschitz continuous (Proposition 3) ()). This is because, for a matrix \(W \in \mathbb{R}^{m \times n}\), a semi-orthogonal matrix \(W\) preserves the \(2\)-norm: \[\begin{align} W^\top\,W = I_n \quad \text{or} \quad WW^\top = I_m \quad \Longrightarrow \quad \mathrm{Lip}\left[ W \right]=1 , \end{align}\] if \(m \geq n\) and \(W^\top W=I_n\), \(W\) is an isometric map; if \(m < n\) and \(WW^\top = I_m\), \(W\) is non-expansive, vanishing on the kernel \(\ker(W)\) (Proposition 3) ()).

4.5.1.1 Extending Semi-Orthogonality to Convolution.

Let \(W \in \mathbb{R}^{c_{out} \times c_{in} \times n \times n}\) be a convolutional filter where \(c_{out}\) is the number of output channels, \(c_{in}\) is the number of input channels, and \(n \times n\) are the input dimensions. The action on an input \(x \in \mathbb{R}^{c_{in} \times n \times n}\) is denoted by: \[\begin{align} W \circledast x: \mathbb{R}^{c_{in} \times n \times n} \to \mathbb{R}^{c_{out} \times n \times n}. \end{align}\] Using the norm preserving concept for \(1\)-Lipschitz maps on \(\ell_2\), for \(\forall\;x \in \mathbb{R}^{c_{in} \times n \times n}\), if the \(2\)-norms are preserved: \[\begin{align} \|W \circledast x\|_2 = \|x\|_2 , \end{align}\] then the convolutional filter \(W\) is referred to as semi-orthogonal [34].

4.5.1.2 Filter Orthogonalization via Cayley Transform.

However, convolutional operation is not matrix multiplication, to apply the existing results from matrix theory, [34] harness the Convolutional Theorem [84] for converting the convolutional operation in spatial domain to multiplication in frequency domain by: \[\begin{align} \mathscr{F}\left[W \circledast x\right][:,i,j] =\mathscr{F}\left[W\right][:,:,i,j] \; \mathscr{F}\left[x\right][:,:,i,j] = \widehat{W}[:,i,j]\widehat{x}[:,i,j] , \end{align}\] where \(\mathscr{F}\) represents Fourier transform, \(\widehat{W}=\mathscr{F}\left[W\right]\), \(\widehat{x}=\mathscr{F}\left[x\right]\), \([:,:,i,j]\) and \([:,i,j]\) are slicing operations by fixing indices \((i,j)\)) [34]. Note that \(\widehat{W}[:,:,i,j]\) is not orthogonal or semi-orthogonal, [34] use a bijective Cayley Transform for deriving an orthogonal matrix \(Q[:,:,i,j]\) from \(\widehat{W}[:,:,i,j]\) by: \[\begin{align} Q[:,:,i,j] = (I - A[:,:,i,j])(I + A[:,:,i,j])^{-1}, \end{align}\] where: \[\begin{align} A[:,:,i,j] = \widehat{W}[:,:,i,j] - \widehat{W}[:,:,i,j]^\ast , \end{align}\] and \(A[:,:,i,j]\) is referred to as skew-Hermitian matrix [34]. Conceptually, the convolution on frequency domain is then computed through: \[\begin{align} \mathscr{F}\left[W \circledast x \right][:,i,j] &= Q[:,:,i,j] \; \widehat{x}[:,i,j] \\ &=(I - A[:,:,i,j])(I + A[:,:,i,j])^{-1} \;\widehat{x}[:,i,j] , \end{align}\] then \(Q[:,:,i,j] \; \widehat{x}[:,i,j]\) is inverted back to spatial domain by: \[\begin{align} \mathscr{F}^{-1}\left[ Q[:,:,i,j] \; \widehat{x}[:,i,j]\right], \end{align}\] which leads to the convolution \(1\)-Lipschitz continuous.

4.5.2 Orthogonality with Lie Unitary Group↩︎

[32] introduce a method, stemming from Lie group theory, through the exponential map, for ensuring \(1\)-Lipschitz by construction. Orthogonal constraints networks can improve the robustness and generalization capabilities [85], [86]. For a matrix \(A\), if \(AA^\ast =I\), then the matrix \(A\) is referred to as a unitary matrix. Unitary matrices allow to solve the exploding or vanish gradients [25], [87]. In the sense of Lie algebra, these unitary matrices form a special orthogonal group on field \(\mathbb{R}\): \[\begin{align} SO(n) = \Bigl\{ A \in \mathbb{R}^{n \times n} \mid A^\top A = I \Bigr\} \end{align}\] and unitary group on field \(\mathbb{C}\): \[\begin{align} U(n) = \Bigl\{ A \in \mathbb{C}^{n \times n} \mid A^\ast A = I \Bigr\} \end{align}\] [32]. We are interested in converting a parameter matrix \(W \in \mathbb{R}^{n \times n}\) to \(SO(n)\) or \(U(n)\). [32] harness the concept of skew-symmetric matrix group (Definition 4 ()) as a bridge: \[\begin{align} \mathfrak{so}(n)=\{B\in\mathbb{R}^{n\times n}\;:\;B=-B^\top\} \end{align}\] and: \[\begin{align} \mathfrak{u}(n)=\{B\in\mathbb{R}^{n\times n}\;:\;B=-B^\ast \} . \end{align}\] Taking: \[\begin{align} B=W - W^\top \end{align}\] can easily send a parameter matrix \(W\) to \(\mathfrak{u}(n)\).

4.5.2.1 Connecting to Lie Orthogonal Group.

Then there exists an exponential surjective map between two groups: \[\begin{align} \exp: \mathfrak{so}(n) \to SO(n) \end{align}\] and: \[\begin{align} \exp: \mathfrak{u}(n) \to U(n) , \end{align}\] which is defined as: \[\begin{align} \exp(B) = I + B + \frac{1}{2}B^2 + \cdots \end{align}\] [32]. Then a network \(f(x; W)\) is parameterized on the Lie group \(SO(n)\) by an algebraic mapping chain: \[\begin{align} W \in \mathbb{R}^{n \times n} \xlongrightarrow[]{B=W-W^\top} B \in \mathfrak{so}(n) \xlongrightarrow[]{A=\exp(B)} A \in SO(n) , \end{align}\] and the optimization problem becomes: \[\begin{align} \min_{A}~f(x; A) \Longleftrightarrow \min_{B}~f(x; \exp(B)) \end{align}\] [32].

4.5.2.2 Optimizing on Lie Orthogonal Group.

The Lie group \(SO(n)\) forms a Riemannian manifold, to optimize a function \(f\) on the Riemannian manifold \(SO(n)\): \[\begin{align} f: \mathcal{X} \times SO(n) \to \mathbb{R} , \end{align}\] where \(\mathcal{X}\) is input space, one may use Riemannian gradient descent [88]. Given a point on the manifold \(A \in SO(n)\), \(\mathcal{T}_A\,SO(n)\) is the tangent space at \(A\) with induced metric, and a gradient \(\Omega \in \mathcal{T}_A\,SO(n)\), then the geodesic update is: \[\begin{align} A \leftarrow A \, \exp\left(-\eta\,A^\ast \, \Omega\right) , \end{align}\] where \(\eta\) is a learning rate \(\eta > 0\) [32]. Computing the Riemannian exponential map is expensive; a cheaper first-order approximation, the Cayley map, may be used for the update instead: \[\begin{align} A \leftarrow A \, \phi\left(-\eta\,A^\ast \, \Omega\right) , \end{align}\] where \(\phi\) is the Cayley map: \[\begin{align} \phi(A) = (I+\frac{1}{2}A)(I-\frac{1}{2}A)^{-1} \end{align}\] [89], [90]. Finally, this induces an update rule on \(B \in \mathfrak{so}(n)\): \[\begin{align} \exp(B) \leftarrow \exp\Bigl(B-\eta\,\nabla\,f\left(x;\exp(B)\right)\Bigr)\, \end{align}\] where \(\nabla\,f\left(x;\exp(B)\right)\) is the usual gradient with Euclidean metric. In a specially designed recurrent neural network (RNN), [32] transform the matrix exponential maps skew-symmetric matrices to orthogonal matrices transforming an optimization problem with orthogonal constraints [32]. Particularly, [32] use Padé approximants and the scale-squaring trick to compute machine-precision approximations of the matrix exponential and its gradient [32]. As a result, the optimization remains in the group and the Lipschitz constants of the updated matrices remain \(1\) by design.

4.5.3 Lipschitz Continuous Transformers↩︎

4.5.3.1 \(L_2\) Self-Attention Transformer.

[55] have shown that dot-product self-attention is non-globally Lipschitz continuous [55]. Referring to equation 3 and equation 4 in Section 2.6, recall that the distance matrix for query \(Q\) and key \(K\) in multi-head dot-product attention matrix are computed through their row-wise dot-product: \[\begin{align} P:= \mathrm{Softmax}\Bigl( \frac{QK^\top}{\sqrt{d_K}} \Bigr) = \mathrm{Softmax}\Bigl( \frac{xW_Q(xW_K)^\top}{\sqrt{d_K}} \Bigr) \propto \exp{\Bigl(\frac{xW_Q(xW_K)^\top}{\sqrt{d_K}}\Bigr)} , \end{align}\] which leads to unbounded Jacobian norm. To solve the issue caused by dot-product attention computation, [55] propose to use \(L_2\) distance instead by: \[\begin{align} P_{ij} \propto \exp{\Bigl(\frac{-\|x_iW_Q - x_jW_K\|_2^2 }{\sqrt{d_K}}\Bigr)} , \end{align}\] where \(x_i\) and \(x_j\) are the \(i\)-th and \(j\)-th tokens of \(x\), and \(P_{ij}\) represents their attention score [55].

[55] further show that, for a sequence length \(n\), input dimension \(d\), and \(h\) heads, its \(2\)-norm Lipschitz bound is: \[\begin{align} \mathrm{Lip}\left[L_2\mathrm{-Attention}\right] \leq \frac{\sqrt{n}}{\sqrt{d_K}} \Bigl[ 4\phi^{-1}(n-1)+1 \Bigr] \Biggl( \sqrt{ \sum_{h_i} \|W_{Q,h_i}\|_2^2\; \|W_{V,h_i}\|_2^2\; } \Biggr) \|W_O\|_2 , \label{equ:kims95dotproduct95bound:l2} \end{align}\tag{10}\] and the \(\infty\)-norm Lipschitz bound: \[\begin{align} \mathrm{Lip}_{\infty}\left[L_2\mathrm{-Attention}\right] \leq \Biggl[ 4\phi^{-1}(n-1) + \frac{1}{\sqrt{d_K}} \Biggr] \; \|W_O^\top\|_{\infty} \; \max_{h_i} \|W_{Q,h_i}\|_{\infty} \|W_{Q,h_i}^\top\|_{\infty} \; \max_{h_j} \|W_{V,h_j}^\top\|_{\infty} , \label{equ:kims95dotproduct95bound:linfty} \end{align}\tag{11}\] where: \[\begin{align} \phi(x)=x\exp(x+1) \notag \end{align}\] is an invertible univariate function on \(x>0\), \(W_{Q,h_i}\) is the query parameter matrix for head \(h_i\), \(W_{V,h_i}\) is the value parameter matrix for head \(h_i\), and \(W_O \in \mathbb{R}^{d \times d}\) is the projection parameter matrix; the per-head key dimension is \(d_K = d/h\) [55]. They further conclude that their results can lead that the \(2\)-norm Lipschitz bound is at the scale: \[\begin{align} \mathrm{Lip}[\mathrm{Attention}] \sim O(\sqrt{n} \log n) , \end{align}\] and the \(\infty\)-norm Lipschitz bound is at the scale: \[\begin{align} \mathrm{Lip}_{\infty}[\mathrm{Attention}] \sim O(\log n) \end{align}\] [55]. The bounds in [55] are complemented by concurrent work [91] with a measure-theoretical framework. [91] show that the \(1\)-norm Lipschitz bound is at the scale: \[\begin{align} \mathrm{Lip}_1[\mathrm{Attention}] \sim O(\sqrt{\log{n}}) \end{align}\] [91].

4.5.3.2 LipsFormer.

[92] demonstrate that enforcing Lipschitz continuity in Transformer architecture [2] is more crucial for ensuring training stability in contrast to the tricks such as learning rate warmup, layer normalization, attention formulation, and weight initialization [92]. Their proposed method, referred to as LipsFormer, replaces the Transformer parts with their Lipschitz continuous counterparts: CenterNorm for LayerNorm, spectral initialization for Xavier initialization, scaled cosine similarity attention for dot-product attention, and weighted residual shortcut for residual connection [92]. These modifications result in a Transformer architecture with a provably bounded Lipschitz constant. Their key architectural design elements are summarized as follows:

1. CenterNorm instead of LayerNorm. LayerNorm [93] is widely used in Transformer, defined as: \[\begin{align} \mathrm{LN}(x) = \gamma \, \odot \, z + \beta , \end{align}\] and: \[\begin{align} z = \frac{y}{\mathrm{std}(y)}, \qquad y = \Bigl(I_d - \frac{1}{d}\boldsymbol{1}_d\boldsymbol{1}_d^\top\Bigr)x , \end{align}\] where \(x \in \mathbb{R}^d\) is input, \(I_d\) is an identity matrix in \(\mathbb{R}^{d \times d}\), \(\boldsymbol{1}_d\) is an all-ones vector in \(\mathbb{R}^d\), \(\mathrm{std}(y)\) is the standard deviation of \(y\), \(\odot\) is element-wise product, \(\gamma\) and \(\beta\) are learnable parameters initialized to \(1\) and \(0\) respectively [92]. The Jacobian of \(z\) with respect to \(x\) is given as: \[\begin{align} \frac{\partial\,z}{\partial\,x} = \frac{1}{\mathrm{std}(y)} \Bigl(I_d - \frac{1}{d}\boldsymbol{1}_d\boldsymbol{1}_d^\top\Bigr) \Bigl( I_d - \frac{yy^\top}{\|y\|_2^2} \Bigr) , \end{align}\] which shows that LayerNorm is not Lipschitz continuous when \(\mathrm{std}(y) \to 0\). To address this problem, [92] introduce CenterNorm, defined as: \[\begin{align} \mathrm{CN}(x) = \gamma\,\odot\,\frac{d}{d-1} \Bigl(I_d-\frac{1}{d}\boldsymbol{1}_d\boldsymbol{1}_d^\top\Bigr)x + \beta , \end{align}\] which admits a Lipschitz constant: \[\begin{align} \mathrm{Lip}\left[\mathrm{CN}\right] = \frac{d}{d-1} \approx 1 \end{align}\] for sufficient large \(d\), \(\gamma=1\), and \(\beta=0\) [92].

2. Scaled Cosine Similarity Attention (SCSA). [55] show that standard dot-product self-attention is not Lipschitz continuous since the Lipschitz constant depends on sequence length which is not bounded [55], see also Section 2.6 (). To address this problem, [92] replace standard dot-product attention part with SCSA, defined as: \[\begin{align} \mathrm{SCSA}(x;\nu,\tau)=\nu PV, \qquad P = \mathrm{softmax}(\tau\,QK^\top), \end{align}\] where \(\nu\) and \(\tau\) are learnable scalars, and \(K,Q,V\) are row-wise normalized in \(\ell_2\) [92].

   Remark 16. It is worth noting that their method reduces the Lipschitz constant of the attention module; however, the module remains not globally Lipschitz continuous since the Lipschitz constant still depends on sequence length \(n\): \[\begin{align} \mathrm{Lip}_2\left[SCSA\right] \leq O(n^2), \qquad \text{and} \qquad \mathrm{Lip}_{\infty}\left[SCSA\right] \leq O(n^2) \end{align}\] [92].

3. Weighted Residual Shortcut. For a residual module: \[\begin{align} h(x) = x + g(x) , \end{align}\] where \(g: \mathbb{R}^d \to \mathbb{R}^d\), \(h: \mathbb{R}^d \to \mathbb{R}^d\) and \(x \in \mathbb{R}^d\), the Lipschitz constant of \(h\) is bounded by: \[\begin{align} \mathrm{Lip}\left[h\right] \leq 1 + \mathrm{Lip}\left[g\right] , \end{align}\] see Section 2.7.3 (). To further reduce the Lipschitz constant introduced by residual module, [92] replace the standard residual module with: \[\begin{align} \mathrm{WRS}(x) = x + \alpha \, \odot\, g(x) , \end{align}\] where \(\alpha\) is a learnable parameter, initialized to a small value \([0.1, 0.2]\) [92]. Assuming the residual branch \(g\) is \(1\)-Lipschitz, the Lipschitz constant of \(\mathrm{WRS}\) in residual module is bounded by: \[\begin{align} \mathrm{Lip}\left[\mathrm{WRS}\right] \leq 1 + \max(\alpha) . \end{align}\]

4. Spectral-based Weight Initialization. [92] further normalize the parameters at initialization by spectral-norm. Let \(W\) be a parameter matrix, the \(W\) is normalized by: \[\begin{align} W \leftarrow \frac{W}{\|W\|_2} \end{align}\] [92].

4.5.4 Jacobian Norm Minimization↩︎

Constraining the Lipschitz constant of Transformer architecture for stable optimization and robustness remains an open research problem. [54] first derive an explicit tight local Lipschitz constant for self-attention module in Transformer architecture by analyzing the Jacobian of softmax [54]. Building on top of the Jacobian analysis, they present JaSMin (Jacobian Softmax norm Minimization) for enhancing Transformer’s robustness by constraining the local Lipschitz constant [54].

4.5.4.1 Bounding Lipschitz through Softmax Jacobian.

Consider the softmax function: \[\begin{align} \mathrm{softmax}: \mathbb{R}^d \to \mathbb{R}^d , \end{align}\] then the Jacobian of the softmax function is: \[\begin{align} \mathcal{M}(P_x) := \nabla \, \mathrm{softmax}\left(x\right) = \mathrm{diag}(P_x) - P_x P_x^\top, \end{align}\] where: \[\begin{align} P_x = \Biggl( \frac{\mathrm{e}^{x_i}}{\sum_j \mathrm{e}^{x_j}} \Biggr) . \end{align}\] [54] incorporate this Jacobian expression of the softmax function into self-attention, leading to a Lipschitz constant bound for a single head: \[\begin{align} \mathrm{Lip}\Bigl[\mathrm{Attention}\Bigr] \leq \|W_V\|_2 \Biggl( \|P\|_2 + 2 \|x\|_2^2 \; \|A\|_2 \; \max_i \|\mathcal{M}(P_{i, :})\|_2 \Biggr) \end{align}\] with: \[\begin{align} P = \mathrm{softmax}\Bigl( x A x^\top \Bigr) \in \mathbb{R}^{n \times n} \qquad \text{and} \qquad A = \frac{W_QW_K^\top}{\sqrt{d_K}} \in \mathbb{R}^{d \times d} \end{align}\] [54]. It is worth noting that the local Lipschitz constant contains a term: \[\begin{align} \mathrm{Lip}\Bigl[\mathrm{Attention}\Bigr] \leq O\Bigl( \max_i \|\mathcal{M}(P_{i, :})\|_2 \Bigr) , \end{align}\] which shows that softmax Jacobian norm determines an upper bound of the local Lipschitz constant.

4.5.4.2 Jacobian Softmax Norm Minimization.

To approximate the Jacobian norm more efficiently, [54] introduce a surrogate function \(g\) that captures the partial ordering of the singular values of softmax mappings. For \(x \in \mathbb{R}_{>0}^n\), let \(x_{(k)}\) be the \(k\)-th largest component of \(x\). Define: \[\begin{align} g_k(x) := x_{(k)}(1 - x_{(k)} + x_{(k+1)}) , \end{align}\] for \(k=1,2,\cdots,n-1\) [54]. Let \(A = \mathrm{diag}(x) - xx^\top\), then the following inequality holds true: \[\begin{align} x_{(1)} \geq g_1(x) \geq \sigma_1(A) \geq x_{(2)} \geq g_2(x) \geq \sigma_2(A) \geq x_{(n)} \geq g_n(x) \geq \sigma_n(A) \geq 0 \end{align}\] [54].

Building on the derived bound, [54] propose an efficient regularization loss expressed in terms of the function \(g\) with two forms: \[\begin{align} \mathcal{L}_{\mathrm{JaSMin}(k=0)} = \sum_{\ell=1}^{L} \sum_{h_i=1}^{h} \max_j\log \Bigl[ g_1\Bigl( P_{j,:}^{\ell, i} \Bigr) \Bigr] , \end{align}\] and: \[\begin{align} \mathcal{L}_{\mathrm{JaSMin}(k)} = \sum_{\ell=1}^{L} \sum_{h_i=1}^{h} \max_j\log \Biggl[ \frac{ g_1\Bigl( P_{j,:}^{\ell, i} \Bigr) }{ g_k\Bigl( P_{j,:}^{\ell, i} \Bigr) } \Biggr] , \end{align}\] where \(P_{j,:}^{\ell, i}\) represents the softmax map for the \(\ell\)-th attention module (\(1\leq \ell \leq L\)), \(i\)-th head (\(1 \leq i \leq h\)), and the \(j\)-th row of the map [54]. They also show that for \(\frac{g_1}{g_k}\) is bounded below \(\gamma\): \[\begin{align} \frac{ g_1\Bigl(P_{j,:}^{\ell,i}\Bigr) }{ g_k\Bigl(P_{j,:}^{\ell,i}\Bigr) } \leq \gamma , \end{align}\] for \(1 \leq \gamma \leq \frac{k}{4}\). Then the softmax Jacobian \(2\)-norm is bounded by: \[\begin{align} \| \mathcal{M}(P_{j,:}^{\ell,i}) \|_2 \leq O\Bigl(\frac{\gamma}{k}\Bigr) \end{align}\] [54].

4.6 Remarks↩︎

The regularization approaches surveyed above also involve a trade-off between enforceability, computational cost, expressive power, and guarantee strength. Weight-based methods such as weight clipping [25] (see Section 4.1) and spectral normalization [27], [57] (see Section 4.1) are simple to implement and scale well to large models, but they typically control only upper bounds induced by layer norms rather than the exact network Lipschitz constant (see Proposition 8 and Proposition 10). Gradient-based methods [26] (see Section 4.2) act more directly on the input–output sensitivity by penalizing gradient norms (see Lemma 1), yet they are often local in nature and can be computationally expensive due to repeated differentiation. Activation- and architecture-based approaches, such as GroupSort [75], invertible residual constructions [61], [79], and Lipschitz-continuous transformer variants [54], [55], [92] (see Sections 4.3 and 4.5), provide stronger structural control by design, but may restrict architectural flexibility or reduce expressive efficiency. Class-margin regularization [83] (see Section 4.4) connects Lipschitz control to certifiable robustness more directly through margin-based objectives, although its effectiveness still depends on the quality of the underlying Lipschitz bound (see Equation equation 9 ).

5 Certifiable Robustness↩︎

Certifiable robustness, a cornerstone of trustworthy deep learning, guarantees that a neural network’s predictions remain consistent under adversarial perturbations within a specified norm ball. Unlike empirical defenses like adversarial training, which lack formal guarantees, Lipschitz-based certification methods leverage the Lipschitz constant to bound sensitivity to input changes, ensuring provable robustness. While some literature, such as [49], have been introduced previously from the perspectives — such as theoretical foundations, estimation methods and regularization approaches, we may revisit some of them in this section from the perspective of certifiable robustness.

5.1 Formal Robustness Guarantees↩︎

To complement the results in Section 4.4, we now discuss formal robustness guarantees for the general case. In the setting of a \(k\)-class classifier \(f\) with Lipschitz constant \(K\), a central question in adversarial defense is to determine the largest perturbation norm that leaves the prediction unchanged.

5.2 Certifiable Robustness via Global Lipschitz Bound↩︎

For a network \(f\) with 1-Lipschitz activations and \(W^{(\ell)}\) is the \(\ell\)-th layer parameter matrix, the spectral-norm product bounds the Lipschitz constant of \(f\) globally: \[\begin{align} \mathrm{Lip}[f] \le \prod_{\ell=1}^L \|W^{(\ell)}\|_2,\notag \end{align}\] providing a certified global bound. This is computationally cheap but typically loose. SeqLip [43] refines this spectral product by incorporating singular-vector alignment between layers, also producing certified global bounds. LipSDP [70] interprets activation functions as gradients of convex potential functions, satisfying certain properties described by quadratic constraints. Therefore, Lipschitz constant estimation problem is treated as a semi-definite program (SDP) problem.

5.3 Certifiable Robustness via Local Lipschitz Bound↩︎

Local Lipschitz bounds are computed over a neighborhood. Fast-Lip [69] analyze the activation patterns of ReLU networks and derive a gradient norm based local Lipschitz bound. CLEVER [29] estimates local constants via extreme value theory on sampled gradient norms; this improves tightness but does not yield a formal certificate. MILP-based methods [39] encode the exact piecewise-linear constraints of ReLU networks within the ball, yielding exact local Lipschitz constants (and thus exact certified radii), but are limited to small networks.

5.4 Certifiable Robustness in Language Models↩︎

5.4.0.1 Proper Distance Metrics for Text.

Language models typically operate on text sequences as inputs, for which the conventional notion of a Lipschitz constant over \(\mathbb{R}\) with \(\ell_p\) norm does not directly apply. The Levenshtein distance [95] measures the number of character replacements, insertions or deletions needed in order to transform a sequence \(x\) to sequence \(y\). It is a proper metric that satisfies all axioms of a metric space:

1. Non-Negativity. \(d_{\mathrm{lev}}(x,y) \geq 0\);

2. Identity. \(d_{\mathrm{lev}}(x,y) = 0 \Longleftrightarrow x =y\);

3. Symmetry. \(d_{\mathrm{lev}}(x,y) = d_{\mathrm{lev}}(y,x)\);

4. Triangle Inequality. \(d_{\mathrm{lev}}(x,z) \leq d_{\mathrm{lev}}(x,y) + d_{\mathrm{lev}}(y,z)\);

where \(d_{\mathrm{lev}}(x,y)\) denotes the Levenshtein distance between two sequences \(x\) and \(y\). Thereby it provides a principled foundation for discussing Lipschitz continuity in the context of language models. For tokenized text — where each vocabulary item is represented as a one-hot vector — a variant of the Levenshtein distance, known as the ERP distance, has been proposed for comparing sequences of one-hot vectors [96].

[97] introduce LipsLev, a method for training convolutional text classifiers that are 1-Lipschitz in the sense of Levenshtein distance for Lipschitz constant based defense [97]. That is, for a classifier \(f: S \to \mathbb{R}^c\), two one-hot vector-valued sequences \(x \in \mathbb{R}^{n \times d} \subset S\) and \(y \in \mathbb{R}^{m \times d} \subset S\) where \(n,m\) are sequence lengths and \(d\) is one-hot vector dimension, the classification margin must hold: \[\begin{align} | f_i(x) - f_j(y)| \leq K_{\mathrm{lev}} \; d_{\mathrm{lev}}(x, y) \end{align}\] for a constant \(K_{\mathrm{lev}} \leq 1\), so that the prediction does not change under perturbation radius \(d_{\mathrm{lev}}(x, y) < R\) [97]. They then train a 1-Lipschitz classifier over text sequences — in the sense of the Levenshtein distance — by normalizing each layer’s outputs through division by its corresponding Lipschitz constant [97].

6 Conclusions↩︎

Lipschitz continuity serves as a cornerstone for enhancing the robustness, generalization, and optimization stability of deep neural networks, providing a rigorous mathematical framework to quantify and control sensitivity to input perturbations. This comprehensive survey synthesizes key insights across theoretical foundations, estimation techniques, regularization approaches, and certifiable robustness, highlighting their interconnections and practical implications. By unifying disparate research threads, we address a critical gap in the literature, offering a cohesive resource for researchers and practitioners. Challenges remain in balancing bound tightness with computational scalability, preserving model expressivity under constraints, and extending certifications to diverse norms and large-scale architectures like transformers, paving the way for future advancements in trustworthy deep learning systems.

Acknowledgments↩︎

This publication has emanated from research conducted with the financial support of Taighde Éireann – Research Ireland under grant number 18/CRT/6223. The authors gratefully acknowledge the partial in-kind support provided by the J.E. Cairnes School of Business & Economics at the University of Galway, Ireland. We especially thank Ping Song of Research Ireland – Centre for Research Training in AI & the University of Galway (Ireland) for the help. For the purpose of Open Access, the author has applied a CC BY public copyright license to any Author Accepted Manuscript version arising from this submission. We also thank reviewers for their constructive comments which can significantly improve our research quality.

7 Proofs: Lipschitz Constants of Activation Functions↩︎

7.1 Proof: Lipschitz Constant of Sigmoid↩︎

7.2 Proof: Lipschitz Constant of Tanh↩︎

7.3 Proof: Lipschitz Constant of Softplus↩︎

7.4 Proof: Lipschitz Constant of Swish↩︎

7.5 Proof: Lipschitz constant of GELU↩︎

7.6 Proof: Lipschitz constant of Softmax↩︎

References↩︎

---

1. Correspondence to: roisincrtai@gmail.com↩︎