On Lipschitz Explosion in Deep Neural Networks with Normalization: Consequences for Optimization and Adversarial Robustness

Introduction

Increasing the depth of neural networks is a key driver of modern performance, boosting image classification [1] and text generation in large language models [2]. Yet increased depth also makes optimization more challenging, motivating architectural stabilizers such as residual connections [1] and normalization layers [3], which are now standard components of modern architectures. In this work, we take a closer look at batch normalization in particular, a seemingly routine component of modern deep neural networks, and study its unexpected implications for optimization, generalization, and adversarial robustness.

A recurring quantity linking depth to both optimization and stability is a model’s (local) Lipschitz constant. Crucially, deep networks come with two different worst-case Lipschitz notions, depending on what is perturbed. Writing $F_\theta$ for the network map with parameters $\theta$ (and using Frobenius / Euclidean norms for concreteness), we consider

The parameter Lipschitz constant $\operatorname{Lip}_{\Theta}$ governs how sharply the network output (and hence the induced training objective) can change as the weights move, and it therefore enters worst-case iteration-complexity guarantees in nonsmooth optimization. The input Lipschitz constant $\operatorname{Lip}_{X}$ governs worst-case sensitivity to data perturbations, and it appears (often via upper bounds) in several generalization and robustness analyses based on margins and stability, including spectrally-normalized bounds and related Lipschitz controls [4, 5]. We use these doubly Lipschitz notions as a lens to study unexpected implications of normalization layers in optimization, generalization, and adversarial robustness.

Main result: both Lipschitz constants can be exponentially large. We show that for deep networks with normalization layers (with batch normalization as the canonical example), both $\operatorname{Lip}_{X}(F_{\theta})$ and $\operatorname{Lip}_{\Theta}(F; X_{\mathrm{in}})$ can grow exponentially with depth $m$, and consequently exponentially with the ambient parameter dimension $d$. This occurs even under strong per-layer norm control (e.g., $\lVert W_{\ell}\rVert_{2} \le1$, and in the input-Lipschitz construction even $W_{\ell}$ orthogonal), and already for linear activations. Put differently, we find that while normalization stabilizes many practical training pipelines [3, 6, 7, 8, 9, 10], it causes worst-case doubly Lipschitz explosion as depth increases.

Implications of parameter-Lipschitz explosion for nonsmooth optimization

For smooth objectives, first-order methods admit dimension-independent rates for reaching approximate stationarity [11]. Neural-network training objectives, however, are typically nonsmooth (e.g., due to ReLU activations, pooling layers or normalizations), so stationarity is captured by generalized-gradient notions. A classical tractable notion is Goldstein stationarity [12], and recent algorithms can find approximate Goldstein stationary points with oracle complexity polynomial in the Lipschitz constant $L$ of the objective [13, 14, 15]. In our setting, $L$ is not an external constant; instead, it is induced by the network and our lower bound shows it can be $\exp(\Omega(d))$ for deep normalized architectures. Consequently, Lipschitz-based “dimension-free” guarantees of the form $\operatorname{poly}(L,1/\varepsilon,1/\delta)$ can unfortunately translate into exponential dependence on parameter $d$ for this model class. Moreover, lower bounds in nonsmooth optimization show that some polynomial dependence on $L$ is unavoidable in the worst case [16], underscoring that one cannot expect a uniform theory for Goldstein stationarity that is both fully general and only polylogarithmic in $L$. Motivated by the optimization barrier created by exploding Lipschitz constants, we consider a relaxed directional notion (probabilistic Goldstein stationarity) and give a basic algorithm that finds such points with oracle complexity that does not depend polynomially on $L$. This begets a new quest to find more refined solution concepts (with probabilistic Goldstein stationarity as a basic starting point) that are theoretically tractable in $\operatorname{poly}(d,1/\varepsilon)$ time and hence in $\operatorname{poly}(d,1/\varepsilon,\log L)$ time.

Implications of input-Lipschitz explosion for generalization and adversarial robustness

The same phenomenon impacts learning-theoretic uses of Lipschitz control. Since $\operatorname{Lip}_{X}(F_{\theta})$ (or upper bounds on it) enters several margin-based generalization bounds and Lipschitz-based stability arguments [4, 5, 17], an exponential lower bound indicates that these worst-case guarantees can deteriorate rapidly with depth, even when each layer is individually norm-controlled. Therefore, we conclude that the Lipschitzness assumption may not be realistic for studying generalization in deep neural networks with normalization layers, since this assumption becomes vacuous exponentially fast with depth.

The input Lipschitzness $\operatorname{Lip}_{X}(F_{\theta})$ lens on batch normalization also sheds new light on recent empirical evidence that batch normalization can have an outsized effect on adversarial robustness [18, 19, 20]. By showing theoretically that the Lipschitzness of deep neural networks grows exponentially fast with depth as a byproduct of batch normalization, we show that normalization layers significantly increase the sensitivity of neural-network outputs to noise in their inputs. Therefore, we can partially explain this empirical observation.

At a high level, our hardness construction for exponential growth of $\operatorname{Lip}_{X}(F_{\theta})$ with depth exploits a rank-separation effect. The witness pair consists of a low-rank matrix and a nearby matrix that is full-rank but nearly rank-deficient. Although these two inputs are close in norm, normalization can amplify their separation across layers, leading to an exponential growth in the input Lipschitz constant. This mechanism suggests a corresponding empirical prediction for adversarial vulnerability under batch normalization: If a clean image, or a patch-based representation of it, has approximate low-rank structure (as suggested by classical low-rank patch models for natural images and by low-dimensional structure observed in datasets such as MNIST, Fashion-MNIST, CIFAR-10, and ImageNet [21, 22, 23, 24, 25]) then adding small i.i.d. noise can generically make the representation full-rank while keeping it close to the low-rank set. Our theory predicts that batch-normalized networks amplify this small perturbation through depth, resulting in a large change in the output. We experimentally validate this prediction (see Figure 1).

Our main contributions are as follows:

• We distinguish two natural Lipschitz notions for deep networks: parameter Lipschitzness and input Lipschitzness. Through this lens, we reveal unexpected implications of normalization layers for optimization, generalization, and adversarial robustness.

• For deep networks with normalization layers, we show that both input and parameter Lipschitz constants can grow exponentially with depth (and hence with parameter dimension) even under strong per-layer norm constraints and even with linear activations.

• We clarify how parameter-Lipschitz explosion propagates to worst-case nonsmooth optimization. Existing guarantees for Goldstein stationarity that are polynomial in the objective Lipschitz constant $L$ become exponential in the parameter dimension for deep normalized networks. This motivates a new search for tractable stationarity notions with complexity polynomial in dimension and accuracy, but only polylogarithmic in $L$. We initiate this search by introducing probabilistic Goldstein stationarity, a basic directional relaxation of Goldstein stationarity derived from its dual formulation. We give a simple function-value algorithm that finds such points with oracle complexity that avoids polynomial dependence on $L$.

• We connect input-Lipschitz explosion to generalization and adversarial robustness. In particular, our results show that worst-case Lipschitz-based generalization, margin, and stability arguments

can become vacuous exponentially fast with depth in normalized networks, even when each layer is norm-controlled.

• More broadly, our construction identifies a simple mechanism for adversarial sensitivity in batch-normalized networks. Normalization can amplify small input differences when two nearby inputs have different rank structure. Consequently, when the data are approximately low rank, perturbations that introduce new singular directions can be magnified across layers, while perturbations that remain within the original singular subspace may be much less disruptive. Our experiments support this prediction across MNIST, Fashion-MNIST, and CIFAR-10. Perturbations that change the input’s singular structure lead to substantially sharper drops in accuracy, confidence, and classification margins than perturbations at the same noise level that stay within the input’s original singular subspace. This suggests a geometry-dependent failure mode of normalized networks and helps explain why batch normalization can have an outsized effect on adversarial robustness.

Exponential Explosion of Two Lipschitz Constants with Dimension

In this section, we study the Lipschitz constant of deep neural networks. In particular, we show that the presence of normalization layers, with a focus on batch normalization, causes the Lipschitz constant of the network to grow exponentially with depth, and consequently with the dimension of the parameter space $d$. We observe this explosion in two notions of Lipschitz continuity: (i) the Lipschitz constant of the deep neural network with fixed weights, viewed as an input–output map with respect to the inputs, and (ii) the Lipschitz constant of the deep neural network with fixed input data, viewed as a function on the parameter space with respect to the weights.

Preliminaries and Notation

Notation. For a matrix $X$, $\lVert X\rVert_F$ is the Frobenius norm and $\lVert X\rVert_2$ the operator norm. For a square matrix $A$, $\lambda_{\min}(A)$ denotes its smallest eigenvalue. We write $X_{i:}$ for the $i$-th row of $X$ and $X_{:j}$ for the $j$-th column of $X$. We denote the set $\{1,2,\ldots,m\}$ by $[m]$.

Batch normalization (Column-wise normalization on $\mathbb{R}^{n\times k}$). Let $X\in\mathbb{R}^{n\times k}$. Define the diagonal matrix of column-squared norms

When every column of $X$ is nonzero (equivalently $D(X)\succ0$), the column-wise batch normalization map is $\operatorname{BN}(X):=XD(X)^{-1/2}$. Equivalently, $(\operatorname{BN}(X))_{ij}=X_{ij}/\lVert X_{:j}\rVert_2$. We emphasize that the batch normalization operator defined above omits the mean subtraction step for convenience, in line with prior formulations and theoretical studies [9, 26, 27]. Our results extend to mean-subtracted batch normalization by an isometric embedding into the mean-zero subspace; Appendix § 8.8 gives the construction. If $D(X)\succ0$, then every column of $\operatorname{BN}(X)$ has unit norm, hence

A deep network with normalization layers and linear activation. Fix an integer $n\ge1$, i.e., the batch size. Let $k_0,k_1,\ldots,k_m$ be positive integers (layer widths / feature dimensions). For each layer $\ell\in\{1,\ldots,m\}$ let $W_\ell\in\mathbb{R}^{k_{\ell-1}\times k_\ell}$ be a weight matrix. Given an input matrix $X^{(0)}\in\mathbb{R}^{n\times k_0}$, define the forward iterates

whenever the normalization is well-defined. The depth-$m$ network map is

Why linear activations? We focus on linear activations because they provide a conservative setting for a lower bound. By removing nonlinearities, we rule out amplification effects that could come from the activation function itself, so any growth we prove must come from the interaction between the weights and normalization. This choice is also consistent with the mean-field analysis of batch-normalized networks in Yang et al. [ [28], where linear activations exhibit the slowest gradient-norm growth among the activation functions considered, even though the growth can still be exponential in depth. Thus, showing exponential Lipschitz growth already in the linear case demonstrates that the phenomenon is not caused by nonlinear activations. We expect the same mechanism to persist for a broader class of standard nonlinearities, we show this result for networks with ReLU activations in Appendix § 8.9 and leave the rest to future work.

Input and parameter Lipschitzness. We use two Lipschitz notions, depending on whether the perturbation is applied to the input or to the weights. First, fix a sequence of weights $(W_1,\ldots,W_m)$ and an input domain $\mathcal{X}$ on which every normalization step is defined. View the network as a map from inputs to outputs. Its input Lipschitz constant is

Second, fix an input batch $X_{\mathrm{in}}\in\mathbb{R}^{n\times k_0}$ and view the network output as a function of the parameters. For a weight tuple $\theta=(W_1,\ldots,W_m)\in\Theta\subseteq\prod_{\ell\in[m]}\mathbb{R}^{k_{\ell-1}\times k_\ell}$, we use the Euclidean norm obtained by vectorizing all weights, $\lVert\theta\rVert_2:=\left(\sum_{\ell=1}^{m}\lVert W_\ell\rVert_F^2\right)^{1/2}$. The corresponding parameter-to-output map and its Lipschitz constant are

The first quantity measures sensitivity to input perturbations, while the second measures sensitivity to changes in the weights.

Remark 2.2. Normalization layer is invariant under positive rescaling since $\operatorname{BN}(\alpha X)=\operatorname{BN}(X)$ for every $\alpha>0$. Our exponential lower bounds are established on compact admissible domains with scales bounded away from zero uniformly in depth. Specifically, the input domains satisfy $1\le\lVert X\rVert_F\le\sqrt{k}$, while every weight matrix in the parameter witnesses has operator norm one.

Explosion of the Input Lipschitz Constant $\operatorname{Lip}_{\mathcal{X}}(F_m)$

We first show that the input Lipschitz constant $\operatorname{Lip}_{\mathcal{X}}(F_m)$ can grow exponentially with depth. This growth is not an artifact of a large or expansive linear map. In fact, it already occurs in the following highly controlled setting with orthogonal weights:

Here each weight matrix is orthogonal $W_\ell\in\mathbb{O}(k)$, so $\lVert W_\ell\rVert_2=1$ and therefore $\lVert XW_\ell\rVert_F\le\lVert X\rVert_F$. Thus, the linear part of every layer is non-expansive.

The setting is deliberately restrictive. The linear maps do not increase Frobenius norm, and by (1), every post-normalization activation has fixed norm, $\lVert X^{(\ell)}\rVert_F=\sqrt{k}$. In particular, the activations remain in a bounded region throughout the network. Nevertheless, we show below that the input–output map can still have an input Lipschitz constant that grows exponentially with the depth $m$. Moreover, this growth is witnessed by uncountably many input pairs.

Our proof construction uses orthogonal weights, so the linear maps themselves do not create any expansion. The instability comes from the interaction between normalization and near rank-deficiency.

We choose two inputs, $X_0$ and $Y_0$, that are exponentially close, with $Y_0$ rank-deficient and $X_0$ only slightly perturbed away from being rank-deficient. After $m$ carefully chosen normalization steps, the trajectory starting from $X_0$ becomes uniformly well-conditioned, while the trajectory starting from $Y_0$ remains rank-deficient. Thus the two outputs are separated by at least $\sqrt{\lambda_{\min}(X_m^{\mathsf{T}}X_m)}$ which is relatively high in a well-conditioned matrix, even though the original perturbation was exponentially small. This creates the exponential separation-to-perturbation ratio. We give the full proof of this result in Appendix § 8.6, and provide a proof sketch in Appendix § 8.§ 1.

The next corollary translates the depth-dependent lower bound into a lower bound in terms of the number of parameters. Thus, even when the parameter count is the measure of complexity, the input Lipschitz constant can still grow exponentially.

Explosion of the Parameter Lipschitz Constant $\operatorname{Lip}_{\Theta}(F_m)$

In Section 2.2, we showed that the input–output map can have an input Lipschitz constant $\operatorname{Lip}_{X}(F_m)$ that grows exponentially with depth. We now show that the same phenomenon occurs when the input is fixed and the network is viewed as a function of its parameters. In other words, normalization can also lead to exponential growth of the parameter Lipschitz constant $\operatorname{Lip}_{\Theta}(F_m)$.

The proof is a direct reduction from the input-Lipschitz construction in Theorem 2.3. We fix the input to be $X_{\mathrm{in}}=I_k$ and encode the two input witnesses from Theorem 2.3 into the first-layer weights. The two parameter choices differ only in this first layer, while all subsequent weights are kept the same as in the input construction. By the scale-invariance of batch normalization, the resulting forward pass reproduces the same two trajectories after the first layer. Thus the parameter perturbation is as small as the original input perturbation, while the output separation remains bounded below, giving the same exponential separation ratio. Full details are provided in Appendix § 8.7. As in Corollary 2.4, this depth-dependent lower bound can also be translated into an exponential lower bound in the parameter dimension $d$.

From the optimization viewpoint, the relevant quantity is the Lipschitz constant of the training objective that is fed by the network. To make this connection explicit, fix a realizable target parameter $\theta^\star\in\Theta$ and define

Here $F_m(\theta^\star)$ plays the role of the observed label generated by the ground-truth network. The objective function $f(\theta)$ is equivalent to the mean-squared loss evaluated on the outputs of the neural network.

The next corollary shows that the exponential growth of $\operatorname{Lip}_{\Theta}(F_m)$ transfers to the Lipschitz constant of this induced objective.

up to the constants in Corollary 2.7.

Thus, the same exponential sensitivity appears in the induced optimization objective and directly affects convergence guarantees that scale with its Lipschitz constant.

Experiments for Lipschitz Explosion

We now empirically validate the exponential input-Lipschitz growth predicted by our theory. The results are shown in Figure 2 on a logarithmic scale. We focus on the setting of Corollary 2.4 with $k = 2$. Since the proof of Theorem 2.3 establishes existence through a nonconstructive probabilistic argument over the weights $\{W_{\ell}\}_{\ell=1}^{m}$, we use a constructive heuristic to select the weights in our experiments. Specifically, at each layer $\ell\in[m]$, we choose $W_{\ell}$ greedily so as to maximize the isotropy functional from Definition B.1. In the two-dimensional setting, this amounts to a deterministic grid search over rotation angles. At each layer, we scan a grid of candidate rotations and select the one that maximizes the isotropy of the next activation.

Vector (PDF)

Figure 2a compares the empirical Lipschitz ratio obtained from the constructed witness pair $(X_0,Y_0)$ with the predicted exponential growth. The same plot also shows witnesses obtained by applying random rotations $Q$ to the input witness pair, yielding pairs of the form $(QX_0,QY_0)$. As predicted by the orthogonal equivariance of the construction, these rotated witnesses exhibit the same exponential growth behavior. Figure 2b compares the constructed witness with generic local input pairs. The random local pairs are generated from i.i.d. standard Gaussian perturbations, and we report both the median and the 95th percentile of their empirical Lipschitz ratios. Unlike the constructed witness, these generic local pairs do not exhibit exponential growth with depth. These experiments support two complementary conclusions. First, the exponential growth predicted by our theory is clearly observable for the constructed witnesses. Second, this behavior is not typical for arbitrary local perturbations. Thus, the witness construction captures a genuine worst-case instability of normalized networks rather than a numerical artifact.

Implications for Optimization Theory in Deep Neural Networks

We now turn to the optimization consequences of the lower bound $L := \operatorname{Lip}(f) \ge\exp(\Omega(d))$ from § 2.3. Dimension-free convergence rates are often viewed as especially desirable in deep learning, since modern networks can have a very large number of parameters. Our results show that, for normalized networks, this interpretation can be misleading. A guarantee may be formally independent of the parameter dimension $d$, but if it depends polynomially on the objective Lipschitz constant $L$, then Corollary 2.8 implies that the same guarantee can still scale exponentially with $d$. In this sense, Lipschitz-dependent “dimension-free” rates need not be truly dimension-free for deep normalized architectures.

This issue is particularly relevant in nonsmooth optimization. Neural-network training objectives are nonconvex, and they are often nonsmooth as well. Even if the outer loss is smooth, standard architectural components such as ReLU [29] and leaky-ReLU activations [30], max-pooling [31], absolute-value regularizers [32], and piecewise-linear combinatorial layers [33, 34] can introduce nonsmoothness into the objective. At the same time, these objectives are typically locally Lipschitz and differentiable almost everywhere, so generalized-gradient notions of stationarity provide the natural language for worst-case analysis. The relevant background on nonsmooth stationarity is recalled in Appendix § 9.§ 1; here, we focus on how the Lipschitz explosion above changes the interpretation of existing convergence guarantees.

Complexity of Finding Goldstein Stationary Points

Existing algorithms for finding $(\delta,\epsilon)$-approximate Goldstein stationary points have oracle complexity that depends polynomially on the objective Lipschitz constant $L$ [13, 14, 35, 15, 36]. Such a dependence is acceptable when $L$ is viewed as a fixed property of the optimization problem. In deep normalized networks, however, $L$ is not an external constant; it is determined by the architecture and can grow with depth. Our results in § 2 show that $L$ can be as large as $\exp(\Omega(d))$. Consequently, a guarantee of the form $O(\operatorname{poly}(L,1/\epsilon,1/\delta))$ can translate into exponential dependence on the parameter dimension $d$, even when the stated rate has no explicit dependence on $d$.

This raises a natural question: can one design algorithms for Goldstein stationarity whose dependence on $L$ is only logarithmic, for example with oracle complexity polynomial in $d$, $1/\epsilon$, $1/\delta$, and $\log L$?

We show that, in full generality, the answer is negative. One cannot hope for a randomized first-order method whose oracle complexity is polynomial in $d$, $1/\epsilon$, $1/\delta$, and only logarithmic in $L$, uniformly over all Lipschitz objectives. The obstruction is not specific to neural networks or to nonconvexity: it already appears for convex Lipschitz functions. We rely directly on Kornowski and Shamir [ [16], which gives a lower bound on the oracle complexity of finding $(\delta,\epsilon)$-Goldstein stationary points of 1-Lipschitz functions.

Toward Tractable Solution Concepts

These observations highlight a mismatch between existing nonsmooth optimization guarantees and the regimes that arise in deep normalized networks. As we observed, current guarantees for Goldstein stationarity do not yield truly polynomial-time bounds once the architecture-induced Lipschitz constant $L$ is taken into account. This suggests a natural quest: to identify weaker, but still meaningful, stationarity notions whose complexity is polynomial in $d$, $1/\epsilon$, and $1/\delta$, and polynomial only in $\log L$. As a first step in this direction, we introduce probabilistic Goldstein stationarity, a simple directional relaxation of Goldstein stationarity. We use it as a basic baseline rather than a final answer, and give a simple function-value algorithm that finds such points with number of queries independent of $L$.

Goldstein stationarity is a worst-case directional condition. To see this, fix $\delta>0$ and let $\partial_\delta f(x)$ denote the Goldstein subdifferential. For a direction $v\in\mathbb{S}^{d-1}$, consider the support function of the set of Goldstein subdifferentials $m_\delta(x;v):=\min_{g\in\partial_\delta f(x)}\langle g,v\rangle$. By Lemma C.3, Goldstein stationarity can be written in the equivalent form $\operatorname{dist}(0,\partial_\delta f(x))\leq\epsilon\Longleftrightarrow\max_{v\in\mathbb{S}^{d-1}}-m_\delta(x;v)\leq\epsilon$. Thus, a $(\delta,\epsilon)$-Goldstein stationary point is one for which no direction has a large negative margin. This dual viewpoint goes back to the work of Goldstein [ [12]. We now relax this worst-case requirement by asking only that large-margin directions occupy a small fraction of the sphere.

We view probabilistic Goldstein stationarity (PGS) as a first step toward understanding which stationarity notions are tractable for deep normalized networks, rather than as a definitive solution concept. Additionally, we give a simple algorithm for finding PGS points that uses only function evaluations and prove its convergence guarantee.

Implications for Generalization and Robustness in Deep Neural Networks

The implications of Lipschitz explosion results for generalization are relatively direct. As discussed in § 1, blindly assuming Lipschitzness $L$ for deep normalized neural networks can make the resulting generalization bounds vacuous with depth (and number of parameters). We now turn to the unexpected implications of our results for the robustness of neural networks.

We next connect Lipschitz explosion to decision flips, showing that it can directly undermine the robustness of deep normalized networks to input noise. The input-Lipschitz lower bound in Theorem 2.3 shows that an exponentially small input perturbation can create a large separation in the network output. One might suspect that this only makes global Lipschitz certificates overly pessimistic, without leading to an actual change in the predicted label. The next two informal consequences show that this is not the case. The same construction can be converted into an explicit binary robustness failure. We defer the formal statements and proofs to Appendix § 10. We first show that a constant-margin decision flip can occur within an exponentially small adversarial radius. This provides a theoretical mechanism that helps explain why batch normalization can have an outsized effect on adversarial robustness [18, 19, 20].

The readout is the hyperplane placed halfway between the two output features generated by the witness pair in the proof of Theorem 2.3. Thus, the result does not rely on a separate instability mechanism. It shows that the Lipschitz explosion itself can produce a genuine small-radius adversarial example once the learned representation is passed through a linear classifier. We now explain the mechanism behind the separation observed in Figure 1. The key point is that the adversarial perturbation is not only small in norm, but also rank-creating. This allows normalization to amplify the perturbation across layers, whereas perturbations that remain within the same low-rank branch do not experience the same amplification. This perspective is also consistent with the empirical geometry of standard image datasets, since MNIST, Fashion-MNIST, and CIFAR-10 images exhibit strong approximate low-rank structure in Figure 3 and Figure 5.

Vector (PDF)

Vector (PDF)

network. Thus perturbations that leave the low-rank branch and perturbations that stay within it have fundamentally different behavior under normalization.

We empirically validate the prediction of Proposition 4.2 in both realistic and controlled settings. In Figure 4, we use a ResNet56 model pretrained on CIFAR-10, while in Figure 8 we use toy CNNs with batch normalization trained on MNIST and Fashion-MNIST. In all experiments, we vary the noise variance and compare perturbations with matched pre-clipping energy but different geometric structure. We consider four perturbation families: i.i.d. pixel-wise noise, i.i.d. column-wise noise, i.i.d. row-wise noise, and singular-value noise, which perturbs only the singular values and therefore keeps the image in the same singular subspace.

Vector (PDF)

Vector (PDF)

The results are consistent with the mechanism suggested by Proposition 4.2. Perturbations that introduce new directions or alter the rank structure of the input lead to a much stronger degradation in accuracy, confidence, and margin. In contrast, same-subspace perturbations are substantially less effective, even at comparable noise energy. This supports the view that the vulnerability induced by normalization is not determined only by perturbation size, but also by whether the perturbation moves the input away from its underlying low-rank structure. We defer the detailed experimental setup and additional results to Appendix § 11.

Conclusion

We showed that normalization layers can induce exponential growth in both input and parameter Lipschitz constants, even under strong per-layer norm control and already with linear activations. This Lipschitz explosion changes the interpretation of several standard guarantees for deep normalized networks, making Lipschitz-based optimization, generalization, and robustness analyses potentially vacuous in worst-case regimes. Our results also identify a rank-separation mechanism behind adversarial sensitivity and show that this mechanism is reflected empirically. We hope these findings motivate new stability notions and analysis tools that better capture the behavior of modern normalized architectures.

Acknowledgement

The authors thank Suvrit Sra for introducing them to nonsmooth optimization and its literature.

A.S. and P.J. were partially supported by the ONR grant N00014-24-1-2470. Part of this research was performed while A.S. was visiting the Institute for Mathematical and Statistical Innovation (IMSI), which is supported by the National Science Foundation (Grant No. DMS-2425650).

Related Works

Deep learning theory. Our study is motivated by the fact that depth is both powerful and difficult to manage. Deeper networks have driven many empirical gains, but they also tend to make optimization more fragile. Many standard tools in modern architectures were introduced to make very deep models easier to train, including normalization layers [3, 7], residual skip connections [1], and modified activation functions [37]. In practice, these tools often stabilize training and make it possible to scale networks to much greater depths. We ask whether this practical stabilization can come with unanticipated costs for training, robustness, and statistical behavior. As a first step toward this broader question, in this work, we study batch normalization.

The Lipschitz constant of a neural network appears in several learning-theoretic analyses through Lipschitz- and margin-sensitive complexity measures. For example, Bartlett et al. [ [4] prove margin-based generalization bounds whose complexity term depends on a spectral complexity of the network, given by the product of the layer spectral norms together with an additional correction factor. Similarly, Bubeck and Sellke [ [5] study robustness through the Lipschitz constant, showing tradeoffs between interpolation, model size, and smoothness. These results make Lipschitz control a natural quantity to study when reasoning about generalization and robustness. Our contribution is complementary. We show that, in deep networks with normalization layers, the relevant input Lipschitz constant can itself grow exponentially with depth, even under strong per-layer norm control. Thus, Lipschitz-dependent generalization and robustness guarantees can become vacuous for deep normalized networks in the worst case.

Our analysis is closely related to mean-field studies of batch normalization. Yang et al. [ [28] show that, in the infinite-width limit, gradient norms in batch-normalized networks can grow exponentially with depth. More recently, Meterez et al. [ [27] showed that this explosion can be avoided for linear activations under suitable full-rank conditions on the input batch. Our results show that this conditioning requirement is essential. Even with linear activations and norm-controlled weights, inputs that are close to rank deficiency can lead to exponential growth of the worst-case input Lipschitz constant, and hence to large worst-case sensitivity. The mechanism is tied to the isotropizing effect of normalization layers. Batch normalization tends to push representations toward a more isotropic geometry [9], but near rank-deficiency this same effect can amplify small perturbations into large output separations.

Smooth optimization. Dimension-free convergence guarantees have a long history in smooth optimization. For smooth nonconvex objectives, gradient descent can find an approximate first-order stationary point at a rate that does not depend on the ambient dimension [11, 38]. A large body of work has tried to go beyond first-order stationarity and reach second-order stationary points, where the gradient is small and there are no directions of significant negative curvature [39, 40, 41]. Under smoothness assumptions, such guarantees are possible in polynomial time. However, dimension-free second-order guarantees cannot hold in full generality. Simchowitz et al. [ [42] show that some dependence on the dimension is unavoidable for general smooth functions. Motivated by this limitation, Daneshmand et al. [ [43] identify noise conditions under which stochastic gradient methods can still achieve dimension-free convergence to second-order stationary points, using structure that is natural in deep learning losses. These results rely on smoothness of the objective, an assumption that is often violated in modern neural networks because of ReLU activations, and max-pooling operations. For this reason, we focus here on nonsmooth optimization, where the objective need not be differentiable everywhere.

Nonsmooth optimization. The study of stationarity for nonsmooth Lipschitz functions goes back to the work of Goldstein [ [12], which introduced what are now known as Goldstein stationary points. Unlike in the smooth setting, exact first-order stationarity is not a finite-time target in general nonsmooth nonconvex optimization. This has led recent work to focus on relaxed notions such as $(\delta,\epsilon)$-Goldstein stationarity. In this direction, Zhang et al. [ [13] gave one of the first non-asymptotic algorithms for finding such points, with complexity scaling polynomially in the Lipschitz constant $L$, the radius parameter $1/\delta$, and the accuracy $1/\epsilon$. More recently, Davis et al. [ [14] developed a gradient-sampling approach under a standard oracle model, improving the dependence on $\epsilon$ and $L$ in certain regimes, at the cost of an additional dependence on the dimension. While these guarantees are polynomial for a fixed Lipschitz constant $L$, this distinction is important in deep learning, where $L$ is not an external constant. It is induced by the network architecture and can grow with depth. Our results show that, for deep networks with normalization layers, this growth can be exponential. Consequently, nonsmooth optimization guarantees that are polynomial in $L$ can become exponential in depth, even when they are formally dimension-free.

Proofs and Details for Section 2

We begin with a proof sketch of Theorem 2.3, and then give the full proof of Theorems 2.3 and 2.6 together with the required machinery and lemmas.

Proof Sketch of Theorem 2.3

We specialize to the square orthogonal regime (3) and the recursion (2). Since $W_\ell\in\mathrm{O}(k)$, the linear maps are non-expansive, so any growth must come from the normalization. Our argument exploits a tension between (i) the ability of $\operatorname{BN}(\cdot)$ to regularize anisotropy when combined with orthogonal mixing, and (ii) the instability of this regularization near rank-deficient inputs.

We follow Joudaki et al. [ [26] and track isotropy via the functional $I(\cdot)$ in Definition B.1. For normalized matrices $X = \operatorname{BN}(\cdot)$ we have $\operatorname{tr}(X^{\mathsf{T}}X) = k$ by (1), hence $I(X) = \det(X^{\mathsf{T}}X)^{1/k} \in(0, 1]$. A key input is the one-step lift inequality of Meterez et al. [ [27]: Theorem B.3 shows that for $X = \operatorname{BN}(\cdot)$ and Haar-uniform $W \in\mathrm{O}(k)$,

where $S(\cdot)$ is the spectral deviation from Definition B.2. Lemma B.4 provides a uniform gap $S(X) \ge s_{k,\eta} > 0$ whenever $I(X) \le1-\eta$. Combining these yields a deterministic multiplicative lift: Corollary B.5 guarantees that, as long as $I(X) \le1-\eta$, we can choose $W$ so that $I(\operatorname{BN}(XW)) \ge r_{k,\eta}I(X)$ with $r_{k,\eta} = (1-s_{k,\eta})^{-1} > 1$. Iterating this choice gives the greedy geometric increase Corollary B.6, so after $m$ layers we can enforce $I(X_m) \ge1-\eta$ provided the initial isotropy is at least $r_{k,\eta}^{-m}(1-\eta)$.

We realize such an initialization while keeping it exponentially close to rank deficiency. As in the appendix proof of Theorem 2.3, we set $G_\varepsilon= (1-\varepsilon)\mathbf{1}\mathbf{1}^{\mathsf{T}} + \varepsilon I$ and choose $X_0$ with $X_0^{\mathsf{T}}X_0 = G_\varepsilon$. Then $\operatorname{diag}(G_\varepsilon) = I$, so $X_0 = \operatorname{BN}(X_0)$, and $I(X_0) = \det(G_\varepsilon)^{1/k} \ge\varepsilon^{(k-1)/k}$. Choosing

ensures $r_{k,\eta}^{m}I(X_0) \ge1-\eta$, hence (by Corollary B.6) we obtain $I(X_m) \ge1-\eta$ along the resulting trajectory. Near-isotropy then forces uniform conditioning: Lemma B.7 yields $\lambda_{\min}(X_m^{\mathsf{T}}X_m) \ge c_k I(X_m)^k \ge c_k(1-\eta)^k =: \delta_{k,\eta}$.

To turn this into a Lipschitz lower bound, we compare $X_0$ to a rank-$(k-1)$ neighbor. Let $X_0 = U\Sigma V^{\mathsf{T}}$ be an SVD with singular values $\sigma_1 \ge\cdots\ge\sigma_k > 0$ and define the matrix $Y_0 := U\operatorname{diag}(\sigma_1,\ldots,\sigma_{k-1},0)V^{\mathsf{T}}$. Then $\operatorname{rank}(Y_0) < k$ and $\lVert X_0-Y_0\rVert_F = \sigma_k = \sqrt{\varepsilon}$. Let $Y_t := F_t(Y_0)$ under the same weights (one may, if desired, interpret BN via the standard $\varepsilon$-regularization from Remark 2.1 to avoid degenerate columns; this does not change the bound). Since each layer is an invertible right-multiplication followed by column-wise rescaling, rank cannot increase, so $\operatorname{rank}(Y_t) < k$ for all $t$. Applying Lemma B.8 at depth $m$ gives $\lVert X_m-Y_m\rVert_F \ge\lambda_{\min}(X_m^{\mathsf{T}}X_m) \ge\delta_{k,\eta}$, and therefore

with $\gamma_{k,\eta} = \frac{k}{2(k-1)}\log r_{k,\eta}$ and $C_{k,\eta}$ as in Theorem 2.3. Finally, the orthogonal equivariance $\operatorname{BN}((QX)W) = Q\operatorname{BN}(XW)$ for $Q \in\mathrm{O}(k)$ implies $F_m(QX) = QF_m(X)$, so the same lower bound is witnessed by the uncountable family of pairs $(QX_0,QY_0)$, as recorded at the end of the proof of Theorem 2.3.

Isotropy and a One-step Lift Inequality

We use the notion of isotropy from Joudaki et al. [ [26]¹, together with the one-step lift developed in Meterez et al. [ [27], to study the dynamics of linear neural networks with normalization across layers. Let $X \in\mathbb{R}^{k \times k}$ and assume $X^{\mathsf{T}}X \succ0$.

A Uniform Multiplicative Lift Below an Isotropy Cap

Fix $\eta\in(0,1)$ and define the compact set of eigenvalue vectors

Consider $S(\lambda):=\frac{1}{2k^2(k+2)}\sum_{i=1}^{k}(\lambda_i-1)^2$ from Definition B.2.

We now prove a uniform spectral gap away from $I(X)=1$.

Proof. $K_{\eta}$ is nonempty, closed, and bounded, hence compact, and $S$ is continuous, so the minimum is attained. Also $S(\lambda)=0$ iff $\lambda_1=\cdots=\lambda_k=1$, but that vector is not in $K_{\eta}$ since its geometric mean is $1>1-\eta$. Thus $s_{k,\eta}>0$. The final statement follows because $S(X)=S(\lambda(X^{\mathsf{T}}X))$ and, for $X=\operatorname{BN}(\cdot)$, $I(X)=\left(\prod_{i=1}^{k}\lambda_i\right)^{1/k}$.

In turn, we can convert this result into a deterministic multiplicative lift for the isotropy functional $I(X)$, as long as it is below $1-\eta$.

Proof. By Lemma B.4, $S(X)\geq s_{k,\eta}$. Theorem B.3 then gives $\mathbb{E}_{W}[I(\operatorname{BN}(XW))\mid X]\geq r_{k,\eta}I(X)$. Since $I(\operatorname{BN}(XW))\geq0$, there must exist a realization of $W$ achieving at least the expectation.

Consequently, we can show that the isotropy functional $I(X)$ increases geometrically until $X$ reaches near-isotropy. In the proof, we choose the weight matrices greedily according to Corollary B.5.

Proof. If $I(X_{t-1})\le1-\eta$, then by construction $I(X_t)\ge r_{k,\eta}I(X_{t-1})$. Once $I(X_{t-1})>1-\eta$ we do not require any further increase. An induction yields the claimed bound.

Near-isotropy Forces a Floor on $\lambda_{\min}$

A useful basic result is the inequality that connects the minimum eigenvalue in the spectrum of $X^\top X$ to the isotropy functional $I(X)$.

Proof. Set $\beta:=\lambda_k$ and note $\sum_{i=1}^{k-1}\lambda_i=k-\beta$. By AM–GM, for fixed sum $k-\beta$ the product of the first $k-1$ eigenvalues is maximized when they are equal, hence

Therefore

which rearranges to $\beta\ge((k-1)/k)^{k-1}I(X)^k$.

Distance to a Rank-deficient Trajectory

We next use a simple geometric fact that turns the rank mismatch between the two trajectories into a quantitative separation. In our construction, one trajectory becomes well-conditioned, while the other remains rank-deficient throughout the network. The following lemma shows that any rank-deficient matrix must stay at least a smallest-singular-value distance away from a well-conditioned full-rank matrix.

Proof. Since $\operatorname{rank}(Y)<k$ there exists a unit vector $v\in\ker(Y)$. Then

Finally, $\lVert Xv\rVert_2^2=v^\top X^\top Xv\ge\lambda_{\min}(X^\top X)$.

Exponential growth of $\operatorname{Lip}_X(F_m)$ with depth

Proof. Fix $m \in\mathbb{N}$ and let $r := r_{k,\eta} > 1$ from Corollary B.5. Define

Step 1: choose an initialization $X_0$ with controlled anisotropy. Let

Then $\operatorname{diag}(G_{\varepsilon})=I$ and the eigenvalues of $G_{\varepsilon}$ are $\lambda_1=k-(k-1)\varepsilon$ and $\lambda_2=\cdots=\lambda_k=\varepsilon$. Choose any $X_0\in\mathbb{R}^{p\times k}$ such that $X_0^{\mathsf{T}}X_0=G_{\varepsilon}$ (e.g. take $X_0=L^{\mathsf{T}}$ where $G_{\varepsilon}=LL^{\mathsf{T}}$ is a Cholesky factorization). Since $\operatorname{diag}(G_{\varepsilon})=I$, every column of $X_0$ has norm $1$, hence $X_0=\operatorname{BN}(X_0)$.

Moreover,

because $k-(k-1)\varepsilon\geq1$ for $\varepsilon\in(0,1]$.

Step 2: pick $\varepsilon$ so that $m$ greedy steps reach near-isotropy. Choose

Then $r^m I(X_0)\geq r^m\varepsilon^{(k-1)/k}=1-\eta$.

Now choose $W_1,\ldots,W_m$ greedily as in Corollary B.6 and set $X_t:=\operatorname{BN}(X_{t-1}W_t)$. By Corollary B.6, $I(X_m)\geq1-\eta$.

Step 3: eigenvalue floor after reaching near-isotropy. Lemma B.7 yields

Step 4: a rank-deficient companion $Y_0$ and distance growth. Let $X_0=U\Sigma V^{\mathsf{T}}$ be an SVD with singular values $\sigma_1\geq\cdots\geq\sigma_k>0$. Define the rank-$(k-1)$ matrix

Then $\operatorname{rank}(Y_0)<k$ and

Let $Y_t:=F_t(Y_0)$. Since each layer is multiplication by an invertible matrix ($W_t$) followed by multiplication by a diagonal matrix (coming from BN), the rank cannot increase, hence $\operatorname{rank}(Y_t)<k$ for all $t$.

Applying Lemma B.8 at depth $m$ gives

Therefore

Uncountably many witness pairs for the same weights. For any orthogonal $Q\in O(k)$, we have for all $X$ and all orthogonal $W$,

because $((QX)W)^{\mathsf{T}}((QX)W)=W^{\mathsf{T}}X^{\mathsf{T}}Q^{\mathsf{T}}QXW=W^{\mathsf{T}}X^{\mathsf{T}}XW$. Hence $F_m(QX)=QF_m(X)$ by induction. Consequently, the ratio $\lVert F_m(QX_0)-F_m(QY_0)\rVert_F/\lVert QX_0-QY_0\rVert_F$ equals the ratio for $(X_0,Y_0)$ for every $Q\in O(k)$. Since $O(k)$ is uncountable, this yields uncountably many witness pairs.

Exponential growth of $\operatorname{Lip}_{\Theta}(F_m)$ with depth

Proof. In the construction of Theorem 2.6, let $\theta^\star$ be one of the two parameter choices witnessing the lower bound on $\operatorname{Lip}_{\Theta}(F_m)$, and let $\theta$ be the other one. Since the target is realizable, $f(\theta^\star)=0$. Therefore,

The right-hand side is exactly the same separation ratio used in the proof of Theorem 2.6. Applying Corollary 2.7 gives the claimed exponential lower bound. $\square$

Extension to Normalization with Mean Subtraction

The results of Sections 2.2 and 2.3 extend to batch normalization with mean subtraction. Under our notation in § 2, the centered normalization map is defined, for $X \in\mathbb{R}^{n \times k}$, by

whenever every centered column is nonzero. This corresponds to the standard batch-variance normalization with a fixed affine scale of $1/\sqrt{n}$.

We now try to embed our hard instances into the centered subspace. Suppose $n \ge k+1$, and choose $U \in\mathbb{R}^{n \times k}$ satisfying

Such a matrix exists since the orthogonal complement of $\mathbf{1}$ in $\mathbb{R}^n$ has dimension $n-1 \ge k$. Left multiplication by $U$ is an isometry and maps every matrix into the subspace of column-centered matrices. In particular, for every $Z \in\mathbb{R}^{k \times k}$,

so the embedding preserves rank, Gram eigenvalues, column norms, and Frobenius distances.

Moreover, since $\mathbf{1}^{\mathsf T}U=0$, mean subtraction leaves $UZW$ unchanged. Consequently, for every $Z$ and $W$ for which normalization is defined,

Thus, centered normalization reproduces the original normalization dynamics exactly on the embedded subspace.

Let $\overline{F}_m$ denote the depth-$m$ network obtained by replacing each normalization layer of $F_m$ with $\operatorname{BN}_{\mathrm{ctr}}$, while keeping the same weight matrices. Applying (5) inductively yields

In particular, applying this identity to the witness pair $(X_0,Y_0)$ from Theorem 2.3 gives

Hence the exponential input Lipschitz lower bound carries over with exactly the same initial perturbation size and final output separation. The Gram-eigenvalue bounds underlying the construction are also preserved by the embedding.

The parameter Lipschitz lower bound follows from the same identity. Indeed, replacing the fixed input $I_k$ in Theorem 2.6 with $U$, and using the same two weight tuples, preserves both their parameter distance and the corresponding output separation. Since the weight matrices remain $k \times k$, the parameter count is unchanged.

Finally, the robustness construction of Theorem 4.1 extends directly to the centered network. Using the readout $UA_m$ and the same offset $b_m$, we have $\lVert UA_m\rVert_F = 1$ and

Thus, the embedded clean and adversarial inputs achieve exactly the same classification margins and perturbation distance as in Theorem 4.1. In particular, the exponential input and parameter sensitivity bounds, as well as the exponentially small adversarial radius, hold for mean-subtracted batch normalization.

Extension to ReLU Activations

Our construction in Appendix § 8.8 also extends to networks with centered batch normalization followed by ReLU activations. We use the standard affine shift of batch normalization, fixed to 2 in every channel, and define

where ReLU acts entrywise and $\mathbf{1}\mathbf{1}^{\mathsf{T}} \in\mathbb{R}^{n \times k}$ denotes the all-ones matrix. Due to normalization layers, each column of $\operatorname{BN}_{\mathrm{ctr}}(HW)$ has unit Euclidean norm, so every entry lies in $[-1,1]$. Consequently, the affine shift makes every entry strictly positive, and ReLU acts as the identity.

The key observation is that mean subtraction eliminates any matrix with identical rows, even after multiplication by a weight matrix. Combining this property with (5), we obtain

Thus, the affine shift introduced by one layer is removed by the centering operation in the next layer, while the centered component follows exactly the dynamics of the original network.

Writing $\widetilde{F}_m = \widetilde{f}_{W_m} \circ\cdots\circ\widetilde{f}_{W_1}$, induction yields

In particular, for the witness pair $(X_0,Y_0)$,

Since $U^{\mathsf{T}}U = I_k$, the input and output distances are preserved exactly. The exponential input Lipschitz lower bound therefore holds with the same witness ratio. Likewise, using the fixed input $U$ and the weight tuples from Theorem 2.6 preserves their parameter distance and output separation, yielding the same parameter Lipschitz lower bound. All affine shifts are fixed, so the number of free weight parameters is unchanged.

Finally, the robustness construction transfers under the same readout $UA_m$ and offset $b_m$. Indeed, since $\mathbf{1}^{\mathsf{T}}U = 0$,

Hence the clean and adversarial inputs keep the same score margins and perturbation distance as in Theorem 4.1.

We conclude that the exponential input and parameter Lipschitz lower bounds, together with the exponentially small adversarial radius, are also attained by networks with centered batch normalization and ReLU activations.

Embedding the Construction at Larger Widths

We show that the exponential lower bounds can be achieved with constants independent of width by embedding the width-two construction into networks of arbitrary width $k \ge2$.

Let $X_0,Y_0$ and $W_1,\ldots,W_m$ be the witnesses from Theorem 2.3 for $k=2$. For any $k\geq2$, consider the block-diagonal matrices

where the identity blocks are omitted when $k=2$. Since normalization acts columnwise and the additional columns have unit norm, the two trajectories at depth $t$ are precisely $\operatorname{diag}(X_t,I_{k-2})$ and $\operatorname{diag}(Y_t,I_{k-2})$. Their input distance and final output separation are therefore identical to those of the width-two construction. Consequently, the input Lipschitz lower bound

holds for every $k\geq2$, with constants independent of width.

Similarly, applying the same embedding to the parameter witnesses in Theorem 2.6, with fixed input $I_k$, preserves their parameter distance, output separation, and unit operator norms. Thus, the parameter Lipschitz lower bound also holds uniformly over width, with the constants of the width-two construction.

These bounds have width-independent exponential rates in depth. When all entries of the weight matrices are counted, however, the parameter dimension is $d=mk^2$, and the input lower bound becomes

Proofs and Details for Section 3

Background on Goldstein Stationarity and First-order Oracle

Notation and Background. Let $f:\mathbb{R}^d\to\mathbb{R}$ be (locally) $L$-Lipschitz with respect to $\|\cdot\|_2$:

Let $\partial f(x)$ denote the Clarke subdifferential at $x$, which coincides with $\{\nabla f(x)\}$ when $f$ is differentiable at $x$ and agrees with the usual convex subdifferential when $f$ is convex. The canonical solution concept studied in nonsmooth nonconvex optimization is that of a Clarke stationary point [44]. A point $x\in\mathbb{R}^d$ is $\epsilon$-Clarke stationary if $\operatorname{dist}(0,\partial f(x))\leq\epsilon$, where

Zhang et al. [ [13] showed that for general nonconvex Lipschitz functions, Clarke stationarity cannot be guaranteed in a finite number of oracle queries in the worst case, even under very strong first-order oracle access. This result shows that the consideration of more relaxed solution concepts, in favor of tractability, is inevitable.

Goldstein Stationarity. A classical and now-standard relaxation is due to Goldstein [ [12]. Fix $\delta>0$. The key idea is to convexify the subdifferential over a neighborhood:

We say that $x$ is $(\delta,\epsilon)$-Goldstein stationary if

This definition recovers the smooth optimality condition when $f$ is differentiable and $\delta$ is sufficiently small. The recent breakthrough of Zhang et al. [ [13] highlighted the tractability of this solution concept in nonsmooth settings, and subsequent works have further developed and studied it extensively [14, 35, 45, 15, 36, 46].

Oracle model. We consider the standard first-order oracle model used in the nonsmooth optimization literature. At each query point $x_t$, the oracle reveals the function value $f(x_t)$ together with subgradient information at $x_t$ (e.g., a Clarke subgradient $g_t\in\partial f(x_t)$ or even the full set $\partial f(x_t)$). The lower bounds stated in the paper holds even when the oracle returns the entire set $\partial f(x_t)$.

Proofs and Details for Section 3.1

We now immediately show that no randomized algorithm can avoid a polynomial dependence on $L$.

Proof. Set $\varepsilon' := \varepsilon/L$. Since $\varepsilon\in(0,L)$, we have $\varepsilon' \in(0,1)$. Applying Theorem C.1 with accuracy $\varepsilon'$ gives a convex 1-Lipschitz function $f:\mathbb{R}^{d} \to\mathbb{R}$, with $d = \widetilde{O}(1/\varepsilon'^{6}) = \widetilde{O}((L/\varepsilon)^{6})$, such that $f(x_{1})-\inf_{x} f(x) \le1$, and no randomized first-order algorithm can find a $(\delta,\varepsilon')$-Goldstein stationary point with constant probability unless

The condition on $\delta$ becomes

Now define $F := Lf$. Then $F$ is convex and $L$-Lipschitz, and

Moreover, by positive homogeneity of the Clarke and Goldstein subdifferentials,

and therefore

Thus, $x$ is $(\delta,\varepsilon)$-Goldstein stationary for $F$ if and only if it is $(\delta,\varepsilon/L)$-Goldstein stationary for $f$, i.e., $(\delta,\varepsilon')$-Goldstein stationary for $f$.

Consequently, any randomized first-order algorithm that finds a $(\delta,\varepsilon)$-Goldstein stationary point of $F$ with probability at least $2/3$ would also find a $(\delta,\varepsilon')$-Goldstein stationary point of $f$ with the same probability. By Theorem C.1, this requires

This proves the claim.

Remark C.2. The regime in Corollary 3.1 is compatible with our exploding-Lipschitz setting. After rescaling the lower bound from [16], the hard instance has dimension $\widetilde{O}((L/\varepsilon)^{6})$ and requires $\delta\le L/(12\varepsilon)$. In our deep normalized networks, $L$ can be $\Omega(\exp(d))$, so these conditions are mild for fixed or polynomially scaled accuracy parameters. Thus the lower-bound regime overlaps with the regime arising in our construction, and substituting $L = \Omega(\exp(d))$ into $\Omega(L^{2}/\varepsilon^{2})$ yields exponential dependence on the network dimension.

Algorithm 1 Probabilistic Greedy Goldstein Descent (PGGD)

Require: $x_0 \in\mathbb{R}^d$, $\delta> 0$, $\tau> 0$, $\rho\in(0,1)$, $\gamma\in(0,1)$, $f_{\mathrm{lb}} \le\inf f$.

1. $\Delta\gets f(x_0)-f_{\mathrm{lb}}$

2. $T_{\max} \gets\left\lceil\frac{\Delta}{\tau\delta} \right\rceil$

3. $m \gets\left\lceil\frac{1}{\rho}\log\frac{T_{\max}+1}{\gamma} \right\rceil$

4. $x \gets x_0$

5. for $t=0$ to $T_{\max}$ do

6. &nbsp;&nbsp;&nbsp;&nbsp;Sample i.i.d. $v_1,\ldots,v_m \sim\operatorname{Unif}(\mathbb{S}^{d-1})$

7. &nbsp;&nbsp;&nbsp;&nbsp;Query $f(x-\delta v_i)$ for all $i$ and let $i^\star\in\arg\min_i f(x-\delta v_i)$

8. &nbsp;&nbsp;&nbsp;&nbsp;if $f(x-\delta v_{i^\star}) \le f(x)-\tau\delta$ then

9. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$x \gets x-\delta v_{i^\star}$

10. &nbsp;&nbsp;&nbsp;&nbsp;else

11. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;return $x$

12. &nbsp;&nbsp;&nbsp;&nbsp;end if

13. end for

14. return $x$

Why exploding $L$ undermines “dimension-free” rates. The discussion in § 2 shows that the Lipschitz constant cannot be treated as a harmless constant in deep normalized networks. Existing guarantees for finding $(\delta,\epsilon)$-Goldstein stationary points scale polynomially in the objective Lipschitz constant $L$, and Corollary 3.1 shows that some polynomial dependence on $L$ is unavoidable in the worst case. But in our setting, $L$ is induced by the network itself. By Theorem 2.6, it can grow exponentially with depth, and hence with the parameter dimension. Therefore, a guarantee of the form $\operatorname{poly}(L,1/\delta,1/\epsilon)$ may be formally independent of dimension, but it can still translate into exponential dependence on the dimension for deep normalized networks.

Proofs and Details for Section 3.2

For fixed $x$ and $\delta>0$, recall

The following elementary identity is the dual form of Goldstein stationarity.

Proof. Let $Q=\partial_\delta f(x)$ and let $g^\star$ be the projection of $0$ onto $Q$. For any unit vector $v$,

Thus $\max_v m_\delta(x;v)\le\operatorname{dist}(0,Q)$. If $g^\star\ne0$, take $v^\star=g^\star/\|g^\star\|_2$. The projection optimality condition gives $\langle g-g^\star,g^\star\rangle\ge0$ for all $g\in Q$, hence

Therefore $m_\delta(x;v^\star)=\|g^\star\|_2=\operatorname{dist}(0,Q)$. If $g^\star=0$, then $\operatorname{dist}(0,Q)=0$, and the stated equivalence for $\epsilon\ge0$ follows immediately.

Proof. Fix $x$. By Davis et al. [ [14], there is a set $\mathcal{V}_x\subset\mathbb{S}^{d-1}$ of full surface measure such that, for every $v\in\mathcal{V}_x$, the function $t\mapsto f(x-tv)$ is differentiable for a.e. $t\in[0,\delta]$ and satisfies

Now suppose $m_{\delta}(x;v)\ge\tau$. For a.e. $t\in[0,\delta]$, whenever $f$ is differentiable at $x-tv$, we have

since $x-tv\in B_{\delta}(x)$ for all $t<\delta$ and the endpoint is irrelevant for the integral. Therefore,

Integrating gives

which is equivalent to the desired inequality.

Proof. Each accepted step decreases $f$ by at least $\tau\delta$. Since $f(x_{0})-f_{\mathrm{lb}}=\Delta$, the algorithm can accept at most $T_{\max}\le\Delta/(\tau\delta)$ steps, and hence performs at most $T_{\max}+1$ stopping checks.

Call a point $x$ bad if

At any bad point, the probability that all $m$ sampled directions fail to give sufficient decrease is at most

A union bound over the stopping checks shows that, with probability at least $1-\gamma$, the algorithm never stops at a bad point. Thus the returned point $x_{\mathrm{out}}$ satisfies

By Lemma C.4, for almost every $V$,

Therefore,

so $x_{\mathrm{out}}$ is $(\delta,\tau,\rho)$-PGS. The query bound follows from at most $T_{\max}+1$ iterations and $m+1$ function evaluations per iteration.

Proofs and Details for Section 4

This appendix makes the two informal robustness statements from § 4 precise. We reuse the rank-separated witness pair constructed in the proof of Theorem 2.3.

Proof. Use the construction in the proof of Theorem 2.3. With $r=r_{k,\eta}>1$ from Corollary B.5, that proof constructs inputs $Y_0$ and $X_0$ such that

It also constructs orthogonal weights for which, writing $Y_m=F_m(Y_0)$ and $X_m=F_m(X_0)$,

where $c_k=((k-1)/k)^{k-1}$ is the constant from Lemma B.7. Set

Then $\lVert X_{\mathrm{adv}}-X_{\mathrm{clean}}\rVert_F\leq a_{k,\eta}e^{-\gamma_{k,\eta}m}$.

It remains to choose the readout. Let

Then $\lVert A_m\rVert_F=1$, and direct calculation gives

Thus both inequalities hold with

Finally, the perturbation $E=X_{\mathrm{adv}}-X_{\mathrm{clean}}$ is feasible for the displayed infimum because $s_m(X_{\mathrm{clean}}+E)=s_m(X_{\mathrm{adv}})\leq0$. This proves the adversarial-radius bound.

Proof. The rank identities for $X_{\mathrm{clean}}$ and $X_{\mathrm{adv}}$ were proved in Theorem 4.1. Now let $Z$ be rank-deficient and suppose the forward pass is well-defined. Each layer maps

Here $W_t$ is invertible and, because the normalization is well-defined, $D(Z_{t-1}W_t)^{-1/2}$ is an invertible diagonal matrix. Hence the rank cannot increase through a layer, so $F_m(Z)$ is rank-deficient.

On the full-rank branch of the construction, the proof of Theorem 2.3 gives

Applying Lemma B.8 with $X=F_m(X_{\mathrm{adv}})$ and $Y=F_m(Z)$ yields

which proves the claim.

Experiments

Low-rank structure in image datasets. To support the rank-separation mechanism studied in the main text, we examine the spectral structure of standard image datasets. As shown in Figure 3 and Figure 5, images from MNIST, Fashion-MNIST, and CIFAR-10 exhibit strong approximate low-rank structure. Their spectra decay rapidly, and most of the spectral energy is concentrated in a small number of singular components. Thus, while the images are not necessarily exactly rank-deficient, they typically lie close to a low-rank set. This observation is consistent with classical low-rank models of image data [21, 22, 23, 24, 25] and motivates our focus on perturbations that move an input away from its low-rank structure.

Perturbation model for MNIST and Fashion-MNIST. In the MNIST and Fashion-MNIST experiments, we view each image as a matrix $X \in[0,1]^{28\times28}$ and compare several perturbation families at matched pre-clipping energy. Given a perturbation $\Delta$, the perturbed image is

where $\Pi_{[0,1]}$ denotes entrywise clipping to the valid pixel range. The noise level is controlled by a variance parameter $\sigma^2$, and before clipping every sampled perturbation is rescaled so that

Thus all noise families have the same root-mean-square pixel energy before clipping, making the curves directly comparable across perturbation types. We consider four perturbation structures. IID pixel noise samples a full-dimensional Gaussian direction and is full rank with probability one, so it directly tests rank-creating perturbations. Singular-value (or same-subspace) noise perturbs only the singular values of the image and therefore stays within the row and column subspaces already present in $X$. Column-offset and row-offset noise broadcast one Gaussian vector across rows or columns, respectively, producing structured rank-one perturbations. Together, these perturbation families allow us to separate the effect of perturbation magnitude from the effect of perturbation geometry.

Rank effects of matched-energy perturbations. We next examine how different perturbation geometries affect the intrinsic rank structure of the images when their pre-clipping energy is matched. As shown in Figure 6 and Figure 7, i.i.d. pixel noise consistently increases effective rank, stable rank, and the number of components needed to capture most of the spectral energy. In contrast, same-subspace perturbations have a much weaker effect, while row- and column-structured perturbations often preserve or even reduce these rank-based quantities. These observations support the mechanism in the main text that perturbations that create new spectral directions are more disruptive to the low-rank structure of images than perturbations that remain aligned with it.

Vector (PDF)

Vector (PDF)

Robustness under matched-energy perturbations. We next evaluate whether the geometric differences between perturbation families translate into different prediction behavior. Figure 8 shows that, even when perturbations are matched in pre-clipping energy, their effect on BatchNorm networks can differ substantially. Across MNIST, Fashion-MNIST, and CIFAR-10, subspace-preserving perturbations tend to preserve accuracy, confidence, and margin for larger noise levels, and they produce lower prediction-flip rates. In contrast, i.i.d. pixel noise and structured row or column perturbations often cause earlier confidence collapse, margin degradation, and prediction instability. This is consistent with the rank-separation mechanism from the main text: perturbations that move images away from their dominant low-rank structure are amplified more strongly than perturbations that remain aligned with that structure.

Batch normalization versus matched controls. We next compare models with batch normalization to matched architectures without batch normalization under the same perturbation families. As shown in Figure 9, BatchNorm models often maintain high accuracy at very small noise levels but can undergo a much sharper accuracy drop once the perturbation variance crosses a moderate threshold. This effect is especially clear on MNIST under i.i.d. pixel noise, column offsets, and row offsets. In contrast, the matched NoBN models tend to degrade more gradually. Same-subspace perturbations remain comparatively less damaging, which is consistent with the rank-separation mechanism discussed in the main text.

Vector (PDF)