Stochastic Rounding Increases Small Singular Values

Linkai Ma1
Department of Computer Science
Purdue University
ma856@purdue.edu
Tingzhou Yu
Department of Mathematics
University of Alberta
tingzho1@ualberta.ca
Petros Drineas
Department of Computer Science
Purdue University
pdrineas@purdue.edu


Abstract

Over the past half-dozen years, stochastic rounding (SR) has regained significant attention as a quantization scheme for low-precision floating-point arithmetic, with applications spanning numerical analysis and modern machine learning systems. Recent work has shown that SR acts as an implicit regularizer by increasing the smallest singular value of extremely tall-and-thin (or, symmetrically, short-and-fat) matrices. In this work, we substantially sharpen and extend this understanding in two directions. First, we show that the regularization effect of SR is not restricted to extreme aspect ratio regimes: it persists for matrices with constant aspect ratio. Second, we demonstrate that SR does not merely regularize the smallest singular value, but instead lifts entire clusters of singular values at the tail of the spectrum. Together, these results provide a more general characterization of stochastic rounding as a spectral regularizer, revealing that its effects extend beyond extremal aspect ratios and act on a broader portion of the singular value spectrum.

1 Introduction↩︎

In recent years, the rapid scaling of machine learning models, such as deep neural networks and Large Language Models, has driven an increasing reliance on low-precision arithmetic to meet the computational and memory demands of modern AI workloads [1], [2]. Quantization techniques [3][5], which map high-precision values to low-bit representations, have become important for both training and inference. Among these techniques, stochastic rounding (SR) has recently re-emerged as a powerful alternative to traditional deterministic rounding schemes [6][8], owing to its unbiasedness and favorable error propagation properties.

SR was introduced in the Numerical Analysis community in the 1950s in [9], [10] and was mostly forgotten until recently. In the early 2020s [11], [12] highlighted its role in mitigating numerical error accumulation, followed by a growing body of work suggesting that stochastic rounding has non-trivial structural effects on the matrices it is applied to. In particular, recent theoretical results [13] demonstrate that SR can act as an implicit regularizer: when applied entrywise to a matrix, it increases the smallest singular value and promotes full-rank structure with high probability. This phenomenon has important, yet under-explored, implications for downstream tasks, including regression, optimization, and deep learning, where spectral properties of matrices play a central role in numerical stability and generalization.

However, despite this progress, our current theoretical understanding of SR-induced regularization remains incomplete in at least two key aspects. First, existing results are largely confined to extreme aspect ratio regimes, such as very tall-and-thin matrices (or, symmetrically, very short-and-fat matrices), where randomness can accumulate across a large number of rows or columns. While these settings are analytically convenient, they do not capture the near-square matrices that frequently arise in modern machine learning pipelines, including self-attention weights, fully-connected and convolutional weight matrices, etc. Second, prior analyses focus almost exclusively on the smallest singular value, treating regularization as a one-dimensional phenomenon. In practice, however, many applications, particularly in deep learning, depend not only on the extremal singular values but also on the structure of the spectral tail, where clusters of small singular values encode fine-grained information. Indeed, recent studies highlight that the transition of this spectral tail into a heavy-tailed distribution serves as a critical signature of implicit self-regularization, often dictating the network’s generalization capacity [14][17].

These limitations raise a research question: To what extent does stochastic rounding reshape the entire tail of the singular value spectrum, and does its regularization effect persist beyond extreme aspect ratio regimes? In this work, we provide a substantially sharper and more general understanding of stochastic rounding as a spectral regularizer. We show that the regularization effect of SR is not restricted to pathological matrix geometries, but instead extends to matrices with constant aspect ratio. Moreover, we move beyond the smallest singular value and demonstrate that SR systematically regularizes clusters of singular values at the tail of the spectrum, revealing a richer and more structured form of implicit regularization than previously recognized.

1.1 The Uniformly Dithered Quantizer and Stochastic Rounding↩︎

Prior to formally stating our contributions, we define our notion of SR via the so-called Uniformly Dithered Quantizer. Indeed, SR admits several concrete instantiations. Throughout this work we focus on the uniformly dithered quantizer as a canonical and widely used choice [18]. Our analysis extends to any SR scheme producing independent, unbiased, bounded entrywise noise. Fixing this instantiation lets us state quantitative bounds in terms of a single resolution parameter \(\rho\).

We adopt the same definition as in Section 2.1 of [18] and define the quantizer for a scalar \(x\). Given a dynamic range \(R\) (with \(|x| \le R\)) and a bit-budget \(B\), the quantizer uses \(M = 2^B\) distinct quantization points \(q_1, q_2, \ldots, q_M\), evenly spaced between \(-R\) and \(R\) with resolution \[\label{eqn:pd1} \rho = \frac{2R}{M-1}.\tag{1}\] The quantization operation \(Q_{\rho}(x)\) is stochastic: for a value \(x \in [q_k, q_{k+1}]\), the output is determined probabilistically by \[\label{eqn:Q95rho} Q_{\rho}(x) = \begin{cases} q_{k+1} & \text{with probability } p = \frac{(x-q_k)}{\rho}, \\ q_k & \text{with probability } 1-p. \end{cases}\tag{2}\] The quantizer enjoys several properties: it is unbiased; its quantization error is bounded by the resolution \(\rho\); the variance of the quantization error is bounded by \(\frac{\rho^2}{4}\); and the quantization error is sub-gaussian (see Appendix 6.2 for details).

1.2 Our contributions↩︎

We summarize our two main contributions below. Let \(\mathbf{A}\in \mathbb{R}^{n \times d}\), with \(n > d\), be a matrix that is stochastically rounded (or quantized) by applying the Uniformly Dithered Quantizer of the previous section with resolution \(\rho\) to form the matrix \(\tilde{\mathbf{A}}\in \mathbb{R}^{n \times d}\) as follows: \[\tilde{\mathbf{A}}_{ij} = Q_{\rho}(\mathbf{A}_{ij}),\;\;\text{for all}\;i,j.\]

Our first contribution (Section 3) is to resolve the main open problem of [13] by showing that the regularization effect of SR applies to matrices with constant aspect ratio. [13] established that SR increases the smallest singular value primarily in very tall-and-thin or short-and-fat regimes. We prove that this phenomenon persists even when the number of rows and columns are of the same order, substantially broadening the applicability of SR-based regularization. Roughly speaking, we prove that the smallest singular value of \(\tilde{\mathbf{A}}\), denoted as \(\sigma_{\min}(\tilde{\mathbf{A}})\), is, with high probability, at least2 \[\sigma_{\min}(\tilde{\mathbf{A}}) \gtrsim \rho \sqrt{n \nu},\] for \(d \leq c_0\cdot n\), for some constant \(c_0\)[13] proved that the same bound holds,3 with comparably high probability, only when \(d \ll n^{\frac{1}{4}}\). Precisely, Corollary 4.2 of [13] requires \(d = o\left((\frac{n}{\log n})^{\frac{1}{4}}\right)\). At \(n = 10^6\) this caps \(d\) by \(16\), enforcing an aspect ratio of at least \(6\times 10^4\). Additionally, their lower order terms decay very slowly, and [13] themselves noted that \(n > 10^{50}\) is required just to drop the lower order term factor below \(\frac{1}{2}\).

Our new bound demonstrates that stochastic rounding is a general-purpose spectral regularizer, rather than one that relies on extreme dimensional imbalances. See Theorem 2 and Corollary 1 for precise statements of our results and further comparisons with the bounds of [13].

Our second contribution (Section 4) moves beyond the focus on the smallest singular value and shows that stochastic rounding regularizes entire clusters of small singular values. Specifically, we establish that SR lifts not just the minimum singular value, but a non-trivial portion of the spectral tail. This reveals that SR induces a spectral tail form of regularization, affecting the geometry of the low-energy subspace of the matrix rather than a single direction. This contribution provides a more accurate description of how stochastic rounding interacts with high-dimensional data and model representations. More precisely, suppose \(\mathbf{A}\) has a large spectral gap at index \(k\), i.e., \(\sigma_k (\mathbf{A}) \gg \sigma_{k+1}(\mathbf{A})\). We analyze the expected excess tail energy in the spectrum of \(\tilde{\mathbf{A}}\), namely \[\label{eqn:Tdef} \mathcal{T}:= \mathbb{E}\left[\sum_{i = k+1}^{d} \sigma_i(\tilde{\mathbf{A}})^2 \right] - \sum_{i = k+1}^{d} \sigma_i(\mathbf{A})^2.\tag{3}\] Here \(\sigma_i(\cdot)\) denotes the \(i\)-th singular value of the respective matrix. Theorem 5 and Corollary 2 show that, excluding lower order terms, the expected excess tail energy is lower bounded by the smallest \((d-k)\) column variances of the noise matrix \(\tilde{\mathbf{A}}-\mathbf{A}\). Our analysis uses resolvents and contour integration to isolate the relevant set of singular values of the matrix, as well as other tools from random matrix theory and high-dimensional probability. An alternative interpretation of our result is that SR injects noise that is diffused across directions, preventing concentration in low-dimensional subspaces and thereby lifting not just a single singular value, but an entire portion of the spectrum.

We conclude by highlighting connections between our results and the AI/ML literature. (See also Appendix 6.1 for a discussion of relevant prior work on SR.) Recent work has explained that the spectral structure of weight and embedding matrices plays a central role in optimization and generalization in modern machine learning systems, particularly through the emergence of heavy-tailed spectra and implicit self-regularization phenomena [14][16]. Our results suggest that stochastic rounding can actively reshape this structure, providing a new mechanism for influencing spectral geometry. While quantization noise is typically viewed as harmful, we show that stochastic rounding introduces structured noise that can improve spectral properties of matrices, in contrast to conventional deterministic quantization schemes [3], [4]. Our findings move towards providing a theoretical foundation for the empirical success of stochastic rounding in low-precision training [8], [19], [20].

2 Background↩︎

2.1 Notation↩︎

Let \([k]\) denote the set of positive integers \({1, 2, \cdots, k}\). Let \(\mathrm{i}\) denote the imaginary unit. Scalars are denoted by lowercase letters (e.g., \(x\), \(y\)), vectors by bold lowercase letters (e.g., \(\mathbf{x}\), \(\mathbf{y}\)), and matrices by bold uppercase letters (e.g., \(\mathbf{A}\), \(\mathbf{B}\)). We use \(\mathbf{I}_k\) to denote the \(k \times k\) identity matrix (or simply \(\mathbf{I}\) when the dimension is clear from context). For a matrix \(\mathbf{A}\in \mathbb{R}^{n \times d}\) with \(n > d\), we denote its column space (range) by \(\mathcal{R}(\mathbf{A})\). We use \(\mathbf{P}_{\mathbf{A}}\) to denote the orthogonal projector onto the column space of \(\mathbf{A}\). We write \(\|\mathbf{A}\|_2\) for the spectral norm and \(\|\mathbf{A}\|_F\) for the Frobenius norm. Let \(\mathbf{A}\in \mathbb{R}^{n \times d}\) with \(n>d\) have rank \(r\). We denote its thin Singular Value Decomposition (SVD) by \(\mathbf{A}= \mathbf{U}\mathbf{\Sigma}\mathbf{V}^\top,\) where \(\mathbf{U}\in \mathbb{R}^{n \times r}\), \(\mathbf{\Sigma}\in \mathbb{R}^{r \times r}\), and \(\mathbf{V}\in \mathbb{R}^{d \times r}\). We use \(\mathbf{u}_i\) (respectively, \(\mathbf{v}_i\)) to denote the \(i\)-th left (respectively, right) singular vector of \(\mathbf{A}\), and \(\sigma_i\) to denote its \(i\)-th singular value. Similarly, for the perturbed matrix \(\tilde{\mathbf{A}}= \mathbf{A}+ \mathbf{E}\) (with \(\mathbf{E}\in \mathbb{R}^{n \times d}\)), we denote its thin SVD by \(\tilde{\mathbf{A}}= \tilde{\mathbf{U}}\tilde{\mathbf{\Sigma}} \tilde{\mathbf{V}}^\top\). Again, \(\tilde{\mathbf{u}}_i\) (respectively, \(\tilde{\mathbf{v}}_i\)) denotes the \(i\)-th left (respectively, right) singular vector of \(\tilde{\mathbf{A}}\), and \(\tilde{\sigma}_i\) denotes its \(i\)-th singular value. We will focus on the asymptotic scaling of terms rather than their exact magnitudes. Therefore, we do not keep track of constants closely; absolute constants such as \(C, c_1, c_2,\) etc., may be absorbed or change value from line to line without explicit mention.

2.2 Noise Model and Perturbation Structure↩︎

Applying \(Q_{\rho}\) 2 entrywise to \(\mathbf{A}\) yields the stochastically rounded matrix \[\tilde{\mathbf{A}}= Q_{\rho} (\mathbf{A}) = \mathbf{A}+ \mathbf{E}, \label{eq:SR-model}\tag{4}\] where the entries \(\mathbf{E}_{ij}\) satisfy several properties: They are independent across \((i,j)\); they have zero mean; they are bounded in absolute value by \(\rho\); their variance is at most \(\frac{\rho^2}{4}\); and they are sub-gaussian. See properties 1720 in Appendix 6.2 for details.

2.2.0.1 Column variance matrix.

Because the entries of \(\mathbf{E}\) are centered and independent, the expected Gram of \(\mathbf{E}\) is diagonal: \(\mathbf{D}:= \mathbb{E}\left[\mathbf{E}^\top \mathbf{E}\right] = \operatorname{Diag}(\nu_1, \dots, \nu_d)\), where \(\nu_j := \sum_{i=1}^{n} \operatorname{Var}(\mathbf{E}_{ij})\). The quantity \(\nu_j\) is the column variance of column \(j\). In particular, \(\nu_j \le n\rho^2/4\) by 19 . We denote the decreasing rearrangement of \(\nu_1,\dots,\nu_d\) as \(\nu_1^\downarrow \ge \cdots \ge \nu_d^\downarrow\).

2.2.0.2 Normalized variance parameters.

Our analysis of \(\sigma_d(\tilde{\mathbf{A}})\) will invoke the following normalized variance quantities, rescaled by the resolution \(\rho\): \[\label{def:nu} \nu := \frac{1}{n\rho^2}\min_{1\le j\le d}\{\nu_j\},\quad \bar\nu := \frac{1}{n\rho^2}\max_{1\le j\le d}\{\nu_j\},\quad \mu := \frac{1}{\rho^2}\max_{1\le i\le n}\sum_{j=1}^{d} \operatorname{Var}(\mathbf{E}_{ij}),\quad \kappa := \max\left\{\mu,\; n\bar\nu\right\}.\tag{5}\]

2.2.0.3 Gram matrices and the perturbation.

We denote the Gram matrices of \(\mathbf{A}\) and \(\tilde{\mathbf{A}}\) by \(\mathbf{G}:= \mathbf{A}^\top \mathbf{A}\) and \(\tilde{\mathbf{G}}:= \tilde{\mathbf{A}}^\top \tilde{\mathbf{A}}\), and the corresponding perturbation by \[\boldsymbol{\Delta}:= \tilde{\mathbf{G}}- \mathbf{G}= \mathbf{A}^\top \mathbf{E}+ \mathbf{E}^\top \mathbf{A}+ \mathbf{E}^\top \mathbf{E}.\] Since the entries of \(\mathbf{E}\) are centered, \(\mathbb{E}\left[\mathbf{A}^\top \mathbf{E}\right] = \mathbb{E}\left[\mathbf{E}^\top \mathbf{A}\right] = 0\), and hence \[\mathbb{E}\left[\boldsymbol{\Delta}\right] = \mathbb{E}\left[\mathbf{E}^\top \mathbf{E}\right] = \mathbf{D}= \operatorname{Diag}(\nu_1, \dots, \nu_d).\] We denote the \(i\)-th largest eigenvalue of \(\mathbf{G}\) as \(\lambda_i = \sigma_i^2(\mathbf{A})\) and of \(\tilde{\mathbf{G}}\) as \(\tilde{\lambda}_i = \sigma_i^2(\tilde{\mathbf{A}})\).

2.3 A Lower Bound on \(\sigma_d(\mathbf{E})\)↩︎

A key ingredient for our analysis of \(\sigma_d(\tilde{\mathbf{A}})\) is a sharp4 lower bound on the smallest singular value of a noise matrix with independent, centered, bounded, nonhomogeneous entries, due to [21]. We restate it here, specialized to our setting.

Theorem 1. Let \(\mathbf{E}\) be as in equation 4 and normalized variance parameters \(\nu, \mu, \kappa\) be as before. Assume \(\kappa \ge \log n\). Then for every \(a > 0\) there exists a constant \(C_a > 0\), depending only on \(a\), such that \[\mathbb{P}\left(\sigma_d(\mathbf{E}) \ge \rho\left(\sqrt{n\nu} - \sqrt{\mu} - C_a \kappa^{1/3}(\log n)^{2/3}\right)\right) \ge 1 - n^{-a}.\]

Set \(\mathbf{Y}:= \mathbf{E}/\rho\). Then \(\mathbf{Y}\in \mathbb{R}^{n \times d}\) has independent centered entries with \(\left|\mathbf{Y}_{ij}\right| \le 1\). Applying [21] to \(\mathbf{Y}\) yields the claim. See Appendix 6.4 for details.


2.4 Resolvent and Contour Integration↩︎

A central tool in our analysis of the small-cluster regularization is the classical contour-integral representation of spectral projectors, which was recently utilized to derive sharp matrix perturbation bounds [22][27]. We record the essential pieces here; the cluster-bound argument later applies them to the Gram matrices \(\mathbf{G}\) and \(\tilde{\mathbf{G}}\) of Section 2.2. For simplicity, assume that \(\mathbf{A}\) and \(\tilde{\mathbf{A}}\) both have full column rank.

2.4.0.1 Resolvent.

For \(z \in \mathbb{C}\) outside the spectrum of \(\mathbf{G}\), the resolvent of \(\mathbf{G}\) is \(\mathbf{R}_\mathbf{G}(z) := (z\mathbf{I}- \mathbf{G})^{-1}.\) Writing the eigendecomposition of the Gram matrix as \(\mathbf{G}= \mathbf{V}\mathbf{\Lambda}\mathbf{V}^\top\) with \(\mathbf{\Lambda}= \operatorname{Diag}(\lambda_1,\dots,\lambda_d)\) and \(\mathbf{V}= [\mathbf{v}_1,\dots,\mathbf{v}_d]\), the resolvent admits the spectral decomposition \[\label{eq:resolvent95decomp} \mathbf{R}_\mathbf{G}(z) = \sum_{i=1}^{d}\frac{\mathbf{v}_i\mathbf{v}_i^\top}{z-\lambda_i}, \qquad \|\mathbf{R}_\mathbf{G}(z)\|_2 = \frac{1}{\mathrm{dist}(z,\,\operatorname{spec}(\mathbf{G}))}.\tag{6}\]

2.4.0.2 Neumann expansion for the perturbed resolvent.

For \(\tilde{\mathbf{G}}= \mathbf{G}+ \boldsymbol{\Delta}\), the factorization \(z\mathbf{I}- \tilde{\mathbf{G}}= (z\mathbf{I}- \mathbf{G})(\mathbf{I}- \mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta})\) gives \(\mathbf{R}_{\tilde{\mathbf{G}}}(z) = (\mathbf{I}- \mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta})^{-1}\mathbf{R}_\mathbf{G}(z)\). Whenever \(\|\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\|_2 < 1\), this unrolls into the Neumann series \[\label{eq:neumann} \mathbf{R}_{\tilde{\mathbf{G}}}(z) = \sum_{\ell=0}^{\infty}\left(\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\right)^\ell \mathbf{R}_\mathbf{G}(z).\tag{7}\]

2.4.0.3 Spectral projectors via contour integrals.

Suppose \(\mathbf{G}\) has a spectral gap at index \(k\), i.e., \(\lambda_k > \lambda_{k+1}\), and write \(g := \lambda_k - \lambda_{k+1} > 0\). Let \(\Gamma \subset \mathbb{C}\) be a positively oriented closed contour that encloses \(\lambda_1,\dots,\lambda_k\) and excludes \(\lambda_{k+1},\dots,\lambda_d\). For any function \(f\) holomorphic inside \(\Gamma\), Cauchy’s integral formula applied term-wise to the spectral decomposition of \(\mathbf{R}_\mathbf{G}(z)\) yields \[\label{eq:contour-projector} \frac{1}{2\pi \mathrm{i}}\oint_\Gamma f(z)\,\mathbf{R}_\mathbf{G}(z)\,dz = \sum_{\ell=1}^{k} f(\lambda_\ell)\,\mathbf{v}_\ell\mathbf{v}_\ell^\top.\tag{8}\] Choosing \(f(z) = 1\) recovers the orthogonal projector onto the top-\(k\) eigenspace of \(\mathbf{G}\), while \(f(z) = z\) and taking traces yields \[\label{eq:trace95contour} \sum_{i=1}^{k}\sigma_i(\mathbf{A})^2 \;=\; \sum_{i=1}^{k}\lambda_i \;=\; \mathop{\mathrm{tr}}\left(\sum_{i=1}^k \lambda_i \mathbf{v}_i \mathbf{v}_i^\top\right) \;=\; \frac{1}{2\pi \mathrm{i}}\oint_\Gamma z\,\mathop{\mathrm{tr}}\left(\mathbf{R}_\mathbf{G}(z)\right)\,dz,\tag{9}\] which expresses the top-\(k\) energy as a contour integral of the resolvent. The analogous formula holds for \(\tilde{\mathbf{G}}\) whenever \(\Gamma\) still separates its top-\(k\) eigenvalues from the rest, so combining it with 7 represents the energy shift on the top-\(k\) block as a convergent series in \(\boldsymbol{\Delta}\).

3 Regularization of the Smallest Singular Value↩︎

In this section, we present and sketch the proof of our first main contribution: A high-probability lower bound on \(\sigma_d(\tilde{\mathbf{A}})\) for stochastically rounded matrices that only requires a constant aspect ratio.

Similar to [13], we start from the inequality \[\label{eq:sigma-decomp} \sigma_d(\tilde{\mathbf{A}}) \ge\; \sigma_d(\mathbf{E}) - \|\mathbf{P}_\mathbf{A}\mathbf{E}\|_2,\tag{10}\] where \(\mathbf{P}_\mathbf{A}\) is the orthogonal projector onto \(\mathcal{R}(\mathbf{A})\). We then proceed to bound the two terms on the right separately. We significantly improve the result of 14 by sharpening the lower bound on \(\sigma_d(\mathbf{E})\) and the upper bound on \(\|\mathbf{P}_\mathbf{A}\mathbf{E}\|_2\). Specifically, the Brailovskaya–van Handel bound (Theorem 1) already controls \(\sigma_d(\mathbf{E})\) from below; it remains to tightly control the projection \(\|\mathbf{P}_\mathbf{A}\mathbf{E}\|_2\) from above. Throughout this section, we set \(r\) to be the rank of \(\mathbf{A}\) and state the projection bound in terms of \(r\) rather than \(d\), so that the result remains sharp when \(\mathbf{A}\) is rank-deficient.

The following lemma gives a net-based bound on the \(\|\mathbf{P}_\mathbf{A}\mathbf{E}\|_2\). We sketch its proof below. For details, see Appendix 6.5.

Lemma 1. Let \(\mathbf{A}, \mathbf{E}\) be as in equation 4 and let \(r = \mathop{\mathrm{rank}}(\mathbf{A})\). There exist absolute constants \(C_0, c_0 > 0\) such that, for every \(u \ge 0\), \[\label{eq:PAE-bound} \mathbb{P}\left(\|\mathbf{P}_\mathbf{A}\mathbf{E}\|_2 > C_0\,\rho\left(\sqrt{d+r} + u\right)\right) \le 2 e^{-c_0 u^2}.\qquad{(1)}\]

Let \(\mathbf{U}\in \mathbb{R}^{n \times r}\) have orthonormal columns spanning \(\mathcal{R}(\mathbf{A})\), so that \(\mathbf{P}_\mathbf{A}= \mathbf{U}\mathbf{U}^\top\) and \(\|\mathbf{P}_\mathbf{A}\mathbf{E}\|_2 = \|\mathbf{U}^\top \mathbf{E}\|_2\). Set \(\mathbf{Y}:= \mathbf{U}^\top \mathbf{E}\in \mathbb{R}^{r \times d}\). For fixed \(\mathbf{x}\in \mathbb{S}^{d-1}\) and \(\mathbf{y}\in \mathbb{S}^{r-1}\), the scalar \(\mathbf{y}^\top \mathbf{Y}\mathbf{x}= (\mathbf{U}\mathbf{y})^\top \mathbf{E}\mathbf{x}\) is a weighted sum of independent centered bounded entries; Hoeffding’s lemma yields the sub-gaussian moment-generating bound \(\mathbb{E}\left[\exp(\lambda \mathbf{y}^\top \mathbf{Y}\mathbf{x})\right] \le \exp(\lambda^2 \rho^2/2)\), and the Chernoff bound gives \(\mathbb{P}\left(\left|\mathbf{y}^\top \mathbf{Y}\mathbf{x}\right| > t\right) \le 2 e^{-t^2/(2\rho^2)}\). Choosing \(1/4\)-nets of \(\mathbb{S}^{d-1}\) and \(\mathbb{S}^{r-1}\) of cardinality at most \(9^d\) and \(9^r\) [28], the two-net lemma [28] combined with a union bound over the at most \(9^{d+r}\) pairs and the choice \(t = C\rho(\sqrt{d+r}+u)\) (with \(C\) large enough) yields ?? .


Next, we derive our main theorem, lower-bounding \(\sigma_d(\tilde{\mathbf{A}})\) by combining 10 with Theorem 1 and Lemma 1 via a union bound.

Theorem 2. Under the noise model of Section 2.2, suppose \(\mathbf{A}\in \mathbb{R}^{n \times d}\) has rank \(r\) and \(\kappa \ge \log n\). Then for every \(a > 0\) there exist a constant \(C_a > 0\) depending only on \(a\), an absolute constant \(C_1 > 0\), and an absolute constant \(c_1 > 0\) such that, for every \(u \ge 0\), \[\label{eq:combined-bound} \mathbb{P}\left(\sigma_d(\tilde{\mathbf{A}}) \ge \rho\left(\sqrt{n\nu} - \sqrt{\mu} - C_a\,\kappa^{1/3}(\log n)^{2/3} - C_1\left(\sqrt{d+r}+u\right)\right)\right) \ge 1 - n^{-a} - 2 e^{-c_1 u^2}.\qquad{(2)}\]

By Theorem 1, on an event \(\mathcal{E}_1\) with \(\mathbb{P}\left(\mathcal{E}_1\right) \ge 1 - n^{-a}\), \[\sigma_d(\mathbf{E}) \ge \rho\left(\sqrt{n\nu} - \sqrt{\mu} - C_a\,\kappa^{1/3}(\log n)^{2/3}\right).\] By Lemma 1, on an event \(\mathcal{E}_2\) with \(\mathbb{P}\left(\mathcal{E}_2\right) \ge 1 - 2 e^{-c_1 u^2}\), \[\|\mathbf{P}_\mathbf{A}\mathbf{E}\|_2 \le C_1\,\rho\left(\sqrt{d+r} + u\right),\] where \(C_1 := C_0\) and \(c_1 := c_0\) are the constants from Lemma 1. On \(\mathcal{E}_1 \cap \mathcal{E}_2\)10 gives \[\sigma_d(\tilde{\mathbf{A}}) \ge \sigma_d(\mathbf{E}) - \|\mathbf{P}_\mathbf{A}\mathbf{E}\|_2 \ge \rho\left(\sqrt{n\nu} - \sqrt{\mu} - C_a \kappa^{1/3}(\log n)^{2/3} - C_1\left(\sqrt{d+r}+u\right)\right).\] A union bound yields \(\mathbb{P}\left(\mathcal{E}_1 \cap \mathcal{E}_2\right) \ge 1 - n^{-a} - 2 e^{-c_1 u^2}\), establishing ?? .


Since \(r \le d\), the bound ?? immediately implies a slightly cruder but rank-free version with \(\sqrt{d+r}\) replaced by \(\sqrt{2d}\). Setting \(u = c_1^{-1/2}\sqrt{a\log n}\) in the resulting inequality and imposing a mild lower bound on \(\nu\) together with a constant aspect-ratio constraint \(d \le \eta_0 n\) yields the following headline result (see also Appendix 6.6).

Corollary 1. Let \(\mathbf{A}\in \mathbb{R}^{n \times d}\). Fix \(a > 0\) and \(\nu_0 > 0\) and let \(\eta_0 = \eta_0(a,\nu_0) > 0\), \(c_0 = c_0(\nu_0) > 0\), and \(n_0 = n_0(a,\nu_0) \in \mathbb{N}\) be constants. Suppose \(\nu \ge \nu_0\), \(d \le \eta_0 n\), and \(n \ge n_0\). Then \[\label{eq:rectangular-bound} \mathbb{P}\left(\sigma_d(\tilde{\mathbf{A}}) \ge c_0\,\rho\sqrt{n}\right) \ge 1 - 3 n^{-a}.\qquad{(3)}\]

In contrast to the \(d = o\left((n/\log n)^{1/4}\right)\) requirement of [13], Corollary 1 accommodates any constant aspect ratio \(d/n \le \eta_0\), thus covering slightly rectangular matrices that dominate numerical linear algebra workloads. We do emphasize that the regime where \(d\approx n\) is still not covered by our bounds.

4 Regularization of a Small Singular Value Cluster↩︎

We now turn to our second main contribution: SR regularizes not just the smallest singular value, but a small cluster of singular values at the tail of the spectrum. Throughout this section we work in the setup of Section 2.2 and assume that \(\mathbf{A}\) has a large spectral gap at index \(k\), \(\sigma_k(\mathbf{A}) \gg \sigma_{k+1}(\mathbf{A}),\) or, equivalently in terms of the Gram-matrix eigenvalues, \(\lambda_k \gg \lambda_{k+1}\). We denote the spectral gap by \(g := \lambda_k - \lambda_{k+1}\). Throughout this section, we assume \(\mathbf{A}\) and \(\tilde{\mathbf{A}}\) both have full column rank. The quantity of interest is the expected excess tail energy \(\mathcal{T}\), as defined in 3 .

Our main results in this section are Theorem 5 and Corollary 2. We present them at the end of the section, after outlining their proofs, as the additional notation and intuition developed in the proof sketches will help the reader better interpret and appreciate the results.

The proof of Theorem 5 (from which Corollary 2 follows) consists of four main steps.

  1. An expected energy identity for the entire spectrum (Theorem 3) reduces a lower bound on \(\mathcal{T}\) to an upper bound on the expected top-\(k\) singular-value energy shift (see 11 ).

  2. Cauchy’s integral formula along a contour \(\Gamma\) separating \(\lambda_1,\ldots,\lambda_k\) from the rest of the spectrum, combined with the Neumann expansion of the perturbed resolvent, writes the expected top-\(k\) shift as a series of contour integrals \(\sum_{\ell\ge 1} I_\ell\) (Lemma 2).

  3. The leading terms \(I_1\) and \(I_2\) are computed exactly via residue calculus and then upper bounded via concentration inequalities (Lemmas 3 and 4); together with Step 1, \(I_1\) produces the leading term \(\sum_{j > k}\nu_j^\downarrow\) in the bound.

  4. The remainder term \(\sum_{\ell\ge 3} I_\ell\) is bounded via triangle inequalities for contour integrals, where we use the high-probability bound on \(\|\boldsymbol{\Delta}\|_2\) from Lemma 6 to bound the integrand.

As a first observation, the expected total squared singular-value mass increases after rounding. This was implicitly stated in Theorem 9 of [29].5 We state it below and prove it in Appendix 6.7.

Theorem 3. Let \(\mathbf{A}\in \mathbb{R}^{n\times d}\) and \(\tilde{\mathbf{A}}= \mathbf{A}+\mathbf{E}\), where \(\mathbf{E}\) has independent and centered entries. Then, \[\mathbb{E}\left[\sum_{i=1}^{d}\sigma_i(\tilde{\mathbf{A}})^2\right] = \sum_{i=1}^{d}\sigma_i(\mathbf{A})^2 + \mathbb{E}\left[\|\mathbf{E}\|_F^2\right].\]

By Theorem 3 and splitting the full singular-value energy into the top and tail parts, \[\label{eq:tail-top-identity} \mathcal{T}= \mathbb{E}\left[\|\mathbf{E}\|_F^2\right] - \left(\mathbb{E}\left[\sum_{i=1}^{k}\sigma_i(\tilde{\mathbf{A}})^2\right] - \sum_{i=1}^{k}\sigma_i(\mathbf{A})^2\right).\tag{11}\] So any upper bound on the change in the top-\(k\) singular-value energy immediately yields a lower bound on \(\mathcal{T}\).

To control the top-\(k\) block we apply the contour-integral framework of Section 2.4 along a contour tailored to the spectral gap at index \(k\). Recall that the spectral gap is defined as \(g := \lambda_k - \lambda_{k+1} > 0.\) Let \(\Gamma \subset \mathbb{C}\) be the positively oriented closed contour formed by joining two half-circles of radius \(g/2\) centered at \(\lambda_1\) and \(\lambda_k\) by two horizontal segments of length \(\lambda_1 - \lambda_k\) at heights \(\pm g/2\) (Figure 1). By construction \(\Gamma\) encloses exactly \(\lambda_1, \ldots, \lambda_k\) and excludes \(\lambda_{k+1}, \ldots, \lambda_d\).

Figure 1: An illustration of the contour \Gamma.

Every point of \(\Gamma\) is at distance at least \(g/2\) from \(\operatorname{spec}(\mathbf{G})\): on either half-circle \(\Gamma\) has radius \(g/2\), while on the horizontal segments \(\left|\Im z\right| = g/2\) and the eigenvalues are real. Hence by 6 , \[\label{eq:resolvent-norm-Gamma} \forall z \in \Gamma, \qquad \mathrm{dist}(z,\,\operatorname{spec}(\mathbf{G})) \ge \frac{g}{2}, \qquad \|\mathbf{R}_\mathbf{G}(z)\|_2 \le \frac{2}{g}.\tag{12}\]

Combining the contour-integral representation 9 of the top-\(k\) energy with the Neumann expansion 7 of the perturbed resolvent yields the following series representation of the expected top-\(k\) energy shift. We state it as Lemma 2 and prove it in Appendix 6.8.

Lemma 2 (Top-\(k\) energy expansion). Suppose \(\|\boldsymbol{\Delta}\|_2 < g/2\). Then \[\label{eq:topk-series} \mathbb{E}\left[\sum_{i=1}^{k}\sigma_i(\tilde{\mathbf{A}})^2\right] - \sum_{i=1}^{k}\sigma_i(\mathbf{A})^2 \;=\; \sum_{\ell=1}^{\infty} I_\ell, \qquad I_\ell \;:=\; \frac{1}{2\pi \mathrm{i}}\oint_\Gamma z\,\mathbb{E}\left[\mathop{\mathrm{tr}}\left(\left(\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\right)^\ell \mathbf{R}_\mathbf{G}(z)\right)\right]\,dz.\qquad{(4)}\]

The next subsection integrates and bounds the leading terms \(I_1, I_2\). For the higher-order tail \(\sum_{\ell\ge 3} I_\ell\), we provide an upper bound via a contour-integral argument.

We collect five lemmas: a closed-form computation of \(I_1\), a residue plus concentration bound on \(I_2\), a high probability upper bound of \(\|\mathbf{E}\|_2\), a high probability upper bound of \(\|\boldsymbol{\Delta}\|_2\), and a Neumann-tail bound on \(\sum_{\ell\ge3} I_\ell\).

Lemma 3 (First-order term). Let \(\mathbf{V}_k := [\mathbf{v}_1, \ldots, \mathbf{v}_k] \in \mathbb{R}^{d \times k}\) collect the top-\(k\) right singular vectors of \(\mathbf{A}\). Then \[I_1 = \mathop{\mathrm{tr}}\left(\mathbf{V}_k \mathbf{V}_k^\top\mathbf{D}\right) \le \sum_{j=1}^{k} \nu_j^\downarrow.\]

Because \(\mathbb{E}\left[\boldsymbol{\Delta}\right] = \mathbf{D}\), residue calculus on \(z/(z-\lambda_i)^2\) around \(\Gamma\) collapses the contour integral to \(\mathop{\mathrm{tr}}(\mathbf{V}_k \mathbf{V}_k^\top\mathbf{D})\). The upper bound follows from Ky Fan’s maximum principle applied to the diagonal matrix \(\mathbf{D}\). For details, see Appendix 6.9.


The proofs of the next two lemmas are in Appendices 6.10 and 6.11, respectively.

Lemma 4 (Second-order term). There is an absolute constant \(C > 0\) such that \[I_2 \le \frac{C}{g}\left[\rho^2 \sum_{i=1}^{k}\sum_{j=k+1}^{d}(\lambda_i + \lambda_j) + \sum_{j=1}^{k}\left(\nu_j^\downarrow\right)^2 + \rho^4 n\, k(d-k)\right].\]

Cyclic property of trace and the spectral decomposition of \(\mathbf{R}_\mathbf{G}(z)\) reduce \(I_2\) to \(\sum_{i=1}^k \sum_{j=k+1}^d \mathbb{E}\left[(\mathbf{v}_i^\top\boldsymbol{\Delta}\mathbf{v}_j)^2\right]/(\lambda_i - \lambda_j)\) via case-by-case residue calculus. For each \((i,j)\), the linear part of \(\boldsymbol{\Delta}= \mathbf{A}^\top\mathbf{E}+ \mathbf{E}^\top\mathbf{A}+ \mathbf{E}^\top\mathbf{E}\) contributes \(C\rho^2(\lambda_i + \lambda_j)\) via Hoeffding, and the quadratic part \(\mathbf{E}^\top\mathbf{E}\) contributes \((\mathbf{v}_i^\top\mathbf{D}\mathbf{v}_j)^2 + C\rho^4 n\) via Hanson–Wright. Summing over \((i,j)\), applying \(\lambda_i - \lambda_j \ge g\), and bounding \(\sum_{i \le k < j}(\mathbf{v}_i^\top\mathbf{D}\mathbf{v}_j)^2\) by \(\sum_{j=1}^{k}(\nu_j^\downarrow)^2\) via Ky Fan’s maximum principle gives the claim. For details, see Appendix 6.10.


The next lemma bounds the operator norm of the error matrix \(\mathbf{E}\), which models the quantization error \(\tilde{\mathbf{A}}-\mathbf{A}\). See Appendix 6.11 for its detailed proof.

Lemma 5 (Operator norm of \(\mathbf{E}\)). Under the noise model of Section 2.2, there exist positive constants \(C,c\) such that, for every \(t\ge 0\), \[\mathbb{P}\left\{ \|\mathbf{E}\|_2\le C\rho(\sqrt{n+d}+t) \right\} \ge 1- 2e^{-ct^2}.\] Moreover, there exists a positive constant \(C'\) such that \[\mathbb{P}\left\{ \|\mathbf{E}\|_2\le C'\rho\sqrt n \right\} \ge 1-2e^{-(n+d)}.\]

A standard \(\varepsilon\)-net argument: choose \(1/4\)-nets \(\mathcal{N}\subset \mathbb{S}^{n-1}\) and \(\mathcal{D}\subset \mathbb{S}^{d-1}\) with cardinalities \(|\mathcal{N}|\le 9^n\) and \(|\mathcal{D}|\le 9^d\), apply Hoeffding’s inequality to the bounded sum \(\mathbf{x}^\top\mathbf{E}^\top\mathbf{y}\) for each fixed pair \(\mathbf{x}\in \mathcal{N}\) and \(\mathbf{y}\in \mathcal{D}\), and union bound over the \(9^{n+d}\) pairs. For details, see Appendix 6.11.


Lemma 6 (Operator norm of \(\boldsymbol{\Delta}\)). Under the noise model of Section 2.2, with \(n \ge d\) and \(\mathbf{A}\) of full column rank, there exist positive constants \(C, c_1, c_2\), such that: \[\mathbb{P}\left(\|\boldsymbol{\Delta}\|_2 \le C (\sigma_1(\mathbf{A}) \rho \sqrt{d} + \rho^2 n)\right) \ge 1 - c_1 e^{-c_2 d}.\]

Decompose \(\|\boldsymbol{\Delta}\|_2\) using triangle inequality. Then apply Lemma 1 and Lemma 5. A union bound yields the final result. For details, see Appendix 6.12.


With Lemma 6 in hand we can define the contraction parameters that drive the Neumann tail bound below: let \(C\) be the constant from Lemma 6 and set \[\label{eq:alpha-beta} \alpha(\rho) \;:=\; \frac{2}{g}\,C\left(\sigma_1(\mathbf{A})\,\rho\sqrt{d} + \rho^2 n\right), \qquad \beta(\rho) \;:=\; \frac{2}{g}\left(2\,\sigma_1(\mathbf{A})\,\rho\sqrt{nd} + \rho^2 nd\right).\tag{13}\] On the high-probability event of Lemma 6, \(\|\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\|_2 \le \alpha(\rho)\) uniformly in \(z\in\Gamma\); on its complement, the deterministic Frobenius bound \(\|\mathbf{E}\|_2\le\|\mathbf{E}\|_F\le\rho\sqrt{nd}\) gives \(\|\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\|_2\le\beta(\rho)\).

Remark 4. We could alternatively use the deterministic bound \(\|\boldsymbol{\Delta}\|_2 \le 2 \sigma_1(\mathbf{A}) \|\mathbf{E}\|_2 + \|\mathbf{E}\|_2^2\) and then invoke Lemma 5. This yields a looser bound with a higher success probability. See Appendix 6.19 for details.

Lemma 7 (Higher-order remainder). Assume \(\alpha(\rho) < 1\) and \(\beta(\rho) < 1\) with \(\alpha,\beta\) as in 13 . Then there is an absolute constant \(C > 0\) such that \[\left|\sum_{\ell=3}^{\infty} I_\ell\right| \le \frac{Cd}{g}\left(2(\lambda_1 - \lambda_k) + \pi g\right)\left(\lambda_1 + \frac{g}{2}\right)\left(\frac{\alpha(\rho)^3}{1 - \alpha(\rho)} + c_1 e^{-c_2 d}\,\frac{\beta(\rho)^3}{1 - \beta(\rho)}\right),\] where \(c_1, c_2\) are the constants from Lemma 6.

Let \(\mathcal{G}\) be the event of Lemma 6, on which \(\|\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\|_2\le \alpha(\rho)\). On \(\mathcal{G}^c\), we have \(\|\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\|_2\le \beta(\rho)\). Using triangle inequality for each integral \(I_\ell\) and summing the resulting geometric series completes the proof. For details, see Appendix 6.13.


We are now ready for our main result.

Theorem 5. Assume \(\alpha(\rho) < 1\) and \(\beta(\rho) < 1\) with \(\alpha, \beta\) as in 13 . Then \[\mathcal{T}\ge \sum_{j=k+1}^{d}\nu_j^\downarrow - \frac{C}{g}\left[\rho^2 \sum_{i=1}^{k}\sum_{j=k+1}^{d}(\lambda_i + \lambda_j) + \sum_{j=1}^{k}\left(\nu_j^\downarrow\right)^2 + \rho^4 n\, k(d-k)\right] - \mathcal{R}(\rho),\] where the higher-order remainder \(\mathcal{R}(\rho)\) is the right-hand side of Lemma 7.

Combine the reduction 11 with the series ?? , observing \(\mathbb{E}\left[\|\mathbf{E}\|_F^2\right] = \mathop{\mathrm{tr}}(\mathbf{D}) = \sum_{j=1}^{d}\nu_j\) and \(\sum_{j=1}^{d}\nu_j - \sum_{j=1}^{k}\nu_j^\downarrow = \sum_{j=k+1}^{d}\nu_j^\downarrow\). Then apply Lemmas 3, 4, and 7. For details, see Appendix 6.14.


Specializing the theorem under mild scaling assumptions yields a clean asymptotic statement. (See Appendix 6.15 for details.)

Corollary 2. There exist absolute constants \(c_1, c_2, c_3, c_g > 0\) such that, if \(\rho\sqrt{nd} \le c_1\sigma_k\), \(\sigma_1 \le c_2\sigma_k\), \(d \ge c_3\sqrt n, g= \lambda_k - \lambda_{k+1} > c_g \lambda_k\), then there is a constant \(C > 0\) such that \[\mathcal{T}\ge \sum_{j=k+1}^{d}\nu_j^\downarrow - C\left(k(d-k)\,\rho^2 + \frac{k\,n^2 \rho^4}{\sigma_k^2} + \frac{d^{5/2}\rho^3}{\sigma_k}\right).\]

Remark 6. When \(n,d,\frac{n}{d}\) are all sufficiently large, the leading term \(\sum_{j=k+1}^{d}\nu_j^\downarrow\) dominates. Equivalently, the expected excess tail energy is captured by the smallest \((d-k)\) column variance of \(\mathbf{E}\). In the special case \(k = d-1\), the leading term of this bound matches that of [13]. See Appendix 6.16.

Specializing further to the square case \(d=n\) removes one assumption and simplifies the bound. (See Appendix 6.17 for details.)

Corollary 3. There exist absolute constants \(c_1, c_2, c_g > 0\) such that, if \(d = n, \rho n \le c_1\sigma_k, \sigma_1 \le c_2\sigma_k\) and \(g>c_g \lambda_k\), then there is a constant \(C > 0\) such that \[\mathcal{T}\ge \sum_{j=k+1}^{n}\nu_j^\downarrow - C\left(k(n-k)\,\rho^2 + \frac{n^{5/2}\rho^3}{\sigma_k}\right).\]

Remark 7. When \(n,\frac{n}{k}\) are all sufficiently large, the leading term \(\sum_{j=k+1}^{d}\nu_j^\downarrow\) dominates, equivalently, the expected excess tail energy is captured by the smallest \((d-k)\) column variances of \(\mathbf{E}\). See Appendix 6.18 for a detailed statement and precise conditions.

5 Conclusions↩︎

In this work, we developed a sharper understanding of stochastic rounding (SR) as a spectral regularizer. Our results leave certain regimes uncovered. Corollary 1 requires a constant aspect ratio \(d \le \eta_0 n\) and therefore does not cover the near-square regime \(d \approx n\). Corollary 2 relies on restrictive assumptions, namely a strong spectral gap, a bounded ratio \(\frac{\sigma_1}{\sigma_k}\), and a small perturbation level \(\rho\sqrt{nd} \le c_1 \sigma_k\). While restrictive, these assumptions enable what is, to our knowledge, the first rigorous analysis showing that stochastic rounding regularizes not only the smallest singular value, but entire clusters at the tail of the spectrum. Understanding whether these assumptions can be relaxed, whether they are information-theoretically necessary, and how frequently they arise in practical machine learning workloads are important questions for future work. The corollary further imposes \(d \ge c_3 \sqrt{n}\), thus excluding extremely tall-and-thin matrices. However, we believe this restriction is an artifact of our analysis, since the looser bound in Corollary 4 dispenses with it. Finally, the nature of our work is theoretical and further empirical validations, beyond [6], [30], of the effects of SR in downstream deep learning applications would be useful.

We conclude by noting that our analysis suggests a unifying perspective in which SR injects noise that diffuses across directions, preventing concentration in low-dimensional subspaces and strengthening weak spectral modes. An important direction for future theoretical work is to understand how these spectral effects demonstrate themselves when SR is used to convert matrices to low precision in optimization problems, ranging from simple least-squares regression to linear and semidefinite programming. A much more ambitious and complex goal would be to understand the impact of SR to generalization error in modern machine learning systems, particularly in low-precision training.

6 Appendix↩︎

6.1 Prior Work↩︎

Stochastic rounding (SR) is a probabilistic approach to rounding that has proven effective in large-scale computations and low-precision arithmetic. Despite its illustrious beginnings in the 1950s [9], [10], [31], SR has been largely overlooked by the numerical analysis community. Crucially, the random nature of SR promotes error cancellation rather than deterministic error accumulation. In many core numerical applications, ranging from floating point summation [32], [33], inner product computation [11], [34] to polynomial evaluation via Horner’s algorithm [34], the accumulation of rounding errors under SR can be modeled as a sum of zero-mean random variables. Consequently, these errors concentrate around zero substantially faster than the worst-case bounds characteristic of deterministic rounding. Furthermore, SR inherently prevents stagnation problems typical of traditional deterministic rounding modes [12], [35], where numerous tiny updates to a comparatively large quantity are completely obliterated. For an extensive overview of the historical development, hardware implementations, and probabilistic error analysis of stochastic rounding, we refer the interested reader to the survey [36]. The SIAM News article “Stochastic Rounding 2.0”6 [37] offers a broader view, connecting the use of low-precision arithmetic with emerging complexity analyses and proposing that quantization-aware techniques may shape the future theoretical landscape of numerical algorithms for ML/AI.

The aforementioned theoretical advantages of SR extend into the optimization of modern machine learning architectures, particularly in the large-scale training of Large Language Models (LLMs) under stringent numerical precision constraints. As the parameter counts of foundational models scale into the billions, the adoption of resource-efficient, low-precision arithmetic formats, such as BF16, FP8, and emerging 4-bit structures like MXFP4, has become critical [38][41]. In these very low-precision regimes, traditional deterministic rounding inherently introduces biased quantization noise, which distorts gradient calculations, hinders algorithmic convergence, and frequently causes training stagnation [7]. Stochastic rounding, conversely, functions as an unbiased quantizer and yields accurate gradient estimates that effectively preserve critical update signals [20]. Recent theoretical and empirical analyses demonstrate that SR fundamentally smooths the optimization landscapes of quantized objectives [19] and acts as an implicit regularizer, yielding convergence bounds that seem to mitigate quantization errors [8], [42].

We now discuss in more detail two prior works that are particularly relevant to this paper: A probabilistic analysis of stochastic rounding as an implicit regularizer of the smallest singular value, and a deterministic perturbation analysis for clusters of small singular values.

6.1.0.1 Stochastic rounding as an implicit regularizer.

[13] established the first probabilistic lower bound on \(\sigma_d(\tilde{\mathbf{A}})\) after SR. Specialized to the uniformly dithered quantizer so that their entrywise bound \(\mathcal{R}\) coincides with our resolution \(\rho\), their main result (Theorem 4.1) states that if the entries of \(\mathbf{E}\) are independent, centered, satisfy \(\left|\mathbf{E}_{ij}\right|\le\rho\), and have minimal normalized column variance \(\nu\) as defined in Section 2.2, and if \(n \ge 836\), then with probability at least \(1 - n^{-c} - 2d^2/n^2\), \[\label{eq:dexter-bound} \sigma_d(\tilde{\mathbf{A}}) \ge \rho\sqrt{n}\left(\sqrt{\nu} - \varepsilon_{n,d}\right), \qquad \varepsilon_{n,d} := \sqrt{\tfrac{d}{n}} + 2d^2\sqrt{\tfrac{\log n}{n}} + \tfrac{C(\log n)^{2/3}}{n^{1/30}}\left(\tfrac{d}{n}\right)^{1/54}.\tag{14}\]

Corollary 4.2 of [13] requires \(d = o\left((n/\log n)^{1/4}\right)\). At \(n = 10^6\) this caps \(d\) by approximately \(16\), enforcing an aspect ratio of at least \(6\times 10^4\). The term \(\frac{(\log n)^{2/3}}{n^{1/30}}\) also decays slowly, and [13] themselves noted that \(n > 10^{50}\) is required just to drop this factor below \(1/2\).

6.1.0.2 Deterministic bounds for small singular value clusters.

Motivated by the empirical observation that downcasting a matrix to lower arithmetic precision tends to lift its smallest singular values, [29] obtained a deterministic lower bound for a cluster of small singular values under perturbation. Specialized to the present setting (Assumption 3.1 of [29]), let \(\mathbf{A}\in \mathbb{R}^{n \times d}\) with \(n \ge d\) have \(\mathop{\mathrm{rank}}(\mathbf{A}) \ge d - r\) for some \(r \ge 1\). Write the full SVD of \(\mathbf{A}\) as \(\mathbf{A}= \mathbf{U}\mathbf{\Sigma}\mathbf{V}^\top\) with \(\mathbf{U}\in \mathbb{R}^{n\times n}\), \(\mathbf{V}\in \mathbb{R}^{d\times d}\) orthogonal and \(\mathbf{\Sigma}\in \mathbb{R}^{n\times d}\) diagonal, and partition commensurately \[\mathbf{\Sigma}= \begin{bmatrix} \mathbf{\Sigma}_1 & \mathbf{0} \\ \mathbf{0} & \mathbf{\Sigma}_2 \\ \mathbf{0} & \mathbf{0} \end{bmatrix}, \qquad \mathbf{E}= \mathbf{U}\begin{bmatrix} \mathbf{E}_{11} & \mathbf{E}_{12} \\ \mathbf{E}_{21} & \mathbf{E}_{22} \\ \mathbf{E}_{31} & \mathbf{E}_{32} \end{bmatrix} \mathbf{V}^\top,\] where \(\mathbf{\Sigma}_1 \in \mathbb{R}^{(d-r)\times(d-r)}\) is nonsingular diagonal and \(\mathbf{\Sigma}_2 \in \mathbb{R}^{r\times r}\) is diagonal. If \(1/\|\mathbf{\Sigma}_1^{-1}\|_2 > 4\|\mathbf{E}\|_2\) and \(\|\mathbf{\Sigma}_2\|_2 < \|\mathbf{E}\|_2\), then for all \(1 \le j \le r\), \[\label{eq:boutsikas-bound} \sigma_{d-r+j}(\mathbf{A}+ \mathbf{E})^2 \ge \lambda_j\left(\mathbf{E}_{32}^\top\mathbf{E}_{32} + (\mathbf{\Sigma}_2 + \mathbf{E}_{22})^\top(\mathbf{\Sigma}_2 + \mathbf{E}_{22}) - \mathbf{R}_3\right) - r_4,\tag{15}\] where \(\mathbf{R}_3\) and \(r_4\) collect third-order and fourth-and-higher order terms in \(\mathbf{E}\).

Motivated by this result, our second contribution examines the same setting, namely a cluster of small singular values separated from the rest by a spectral gap. However, we approach the problem from the perspective of stochastic rounding. Instead of providing a deterministic lower bound for each singular value in the small cluster, we analyze the expected excess tail energy (as defined in 3 ) directly. We show that the leading term of the error corresponds to the sum of the \(r\) smallest column variances of the noise matrix \(\mathbf{E}\), namely7 \(\sum_{j=d-r+1}^d \nu_j^\downarrow\), using the notation of Section 2.2. This provides a much more interpretable characterization than the bound given in 15 .

6.2 The Uniformly Dithered Quantizer↩︎

We adopt the same definition as in Section 2.1 of [18]. We first define the quantizer for a scalar \(x\). Given a dynamic range \(R\) (with \(|x| \le R\)) and a bit-budget \(B\), the quantizer uses \(M = 2^B\) distinct quantization points \(q_1, q_2, \dots, q_M\), evenly spaced between \(-R\) and \(R\) with resolution \(\rho = \frac{2R}{M-1}\).

The quantization operation \(Q_{\rho}(x)\) is stochastic: for a value \(x \in [q_k, q_{k+1}]\), the output is determined probabilistically by \[Q_{\rho}(x) = \begin{cases} q_{k+1} & \text{with probability } p, \\ q_k & \text{with probability } 1-p, \end{cases}\] where \(p = \frac{x-q_k}{\rho}\).

Before stating the properties of the quantizer, we recall the sub-gaussian norm of a real-valued random variable \(X\) [28]: \[\label{def:subn} \|X\|_{\psi_2} := \inf\left\{t > 0 : \mathbb{E}\left[\exp(X^2/t^2)\right] \le 2\right\}.\tag{16}\]

The quantizer has the following four statistical properties:

  • Unbiased: the quantizer is unbiased, \[\mathbb{E}\left[Q_{\rho}(x)\right] = x. \label{prop:quant-unbiased}\tag{17}\]

  • Bounded Error: the quantization error is bounded by the resolution, \[\left|Q_{\rho}(x) - x\right| \le \rho. \label{prop:quant-bounded}\tag{18}\]

  • Bounded Variance: the variance of the quantization error is bounded, \[\mathbb{E}\left[(Q_{\rho}(x) - x)^2\right] \le \frac{\rho^2}{4}. \label{prop:quant-var}\tag{19}\]

  • Sub-Gaussian: the quantization error is sub-gaussian, \[\|Q_{\rho}(x) - x\|_{\psi_2} \le C\rho, \label{prop:quant-subg}\tag{20}\] for an absolute constant \(C = 1/\sqrt{\log 2}\).

Properties 17 and 19 are established in Lemma C.1 of [18]; property 18 is immediate from the definition; property 20 is proven in Appendix 6.3.


6.3 Proof of the Sub-Gaussian Property of the Quantization Error↩︎

Here we prove property 20 : the quantization error \(Q_{R,B}(x) - x\) is sub-gaussian with \(\|Q_{R,B}(x) - x\|_{\psi_2} \le C\rho\).

Recall that the sub-gaussian norm is defined in 16 . If \(\left|X\right|\le\rho\), let \(t=\frac{\rho}{\sqrt{\log 2}}\), we get \[\frac{X^2}{t^2}\le \log 2,\] and therefore \[\mathbb{E}\left[\exp(X^2/t^2)\right]\le 2.\] Hence \[\|X\|_{\psi_2}\le \frac{\rho}{\sqrt{\log 2}}\le C\rho.\] In particular, applying this to \(X = Q_{R,B}(x) - x\), which satisfies \(\left|X\right| \le \rho\) by 18 , we obtain \[\|Q_{R,B}(x) - x\|_{\psi_2}\le C\rho.\] Applying this to each entry of \(\mathbf{E}\) gives \(\|\mathbf{E}_{ij}\|_{\psi_2} \le C\rho\) for all \(i, j\).


6.4 Proof of Theorem 1↩︎

Recall that \(\rho\) is the quantization resolution. Set \(\mathbf{Y}:= \mathbf{E}/\rho\). Then \(\mathbf{Y}\in \mathbb{R}^{n \times d}\) has independent centered entries with \(\left|\mathbf{Y}_{ij}\right| \le 1\). For each entry \(\mathbf{A}_{ij}\), let \(q_{ij}\) be the lower grid point such that \(\mathbf{A}_{ij}\in [q_{ij},q_{ij}+\rho]\), and set \[p_{ij}:=\frac{\mathbf{A}_{ij}-q_{ij}}{\rho}\in[0,1].\] Under the uniformly dithered quantizer, \[\tilde{\mathbf{A}}_{ij}=q_{ij}+\rho B_{ij}, \qquad B_{ij}\sim \mathrm{Bernoulli}(p_{ij}),\] with the variables \(B_{ij}\) independent over \((i,j)\). Hence \[\mathbf{E}_{ij}=\tilde{\mathbf{A}}_{ij}-\mathbf{A}_{ij} =\rho(B_{ij}-p_{ij}),\] and therefore, with \(\mathbf{Y}:=\mathbf{E}/\rho\), \[\mathbf{Y}_{ij}=B_{ij}-p_{ij}.\] In particular, the entries of \(\mathbf{Y}\) are independent, centered, with variance \(p_{ij}(1-p_{ij})\). Thus \(\mathbf{Y}\) is a nonhomogeneous centered Bernoulli matrix. Applying [21] to \(\mathbf{Y}\) with variance parameters \[\min_{1\le j\le d}\sum_{i=1}^n \operatorname{Var}(\mathbf{Y}_{ij}) = n\nu, \qquad \max_{1\le i\le n}\sum_{j=1}^d \operatorname{Var}(\mathbf{Y}_{ij}) = \mu,\] and \[\max\left\{ \max_{1\le i\le n}\sum_{j=1}^d \operatorname{Var}(\mathbf{Y}_{ij}), \max_{1\le j\le d}\sum_{i=1}^n \operatorname{Var}(\mathbf{Y}_{ij}) \right\} = \kappa\] yields the claim.


6.5 Proof of Lemma 1↩︎

The proof uses two standard tools from high-dimensional probability.

Lemma 8 (Two-net lemma, [28] Lemma 4.4.2). Let \(\mathbf{Y}\in \mathbb{R}^{d_1 \times d_2}\) and \(\varepsilon \in [0, 1/2)\). For any \(\varepsilon\)-nets \(\mathcal{N}_1 \subset \mathbb{S}^{d_1 - 1}\) and \(\mathcal{N}_2 \subset \mathbb{S}^{d_2 - 1}\), \[\|\mathbf{Y}\|_2 \le \frac{1}{1 - 2\varepsilon}\max_{\mathbf{x}\in \mathcal{N}_2,\, \mathbf{y}\in \mathcal{N}_1} \left|\mathbf{y}^\top \mathbf{Y}\mathbf{x}\right|.\]

Lemma 9 (Hoeffding’s lemma, [43] Lemma 2.2). Let \(X\) be a real-valued random variable with \(\mathbb{E}\left[X\right] = 0\) and \(a \le X \le b\) almost surely. Then for all \(\lambda \in \mathbb{R}\), \[\mathbb{E}\left[\exp(\lambda X)\right] \le \exp\left(\frac{\lambda^2 (b-a)^2}{8}\right).\]

[Proof of Lemma 1] Let \(r := \mathop{\mathrm{rank}}(\mathbf{A})\) and let \(\mathbf{U}\in \mathbb{R}^{n \times r}\) have orthonormal columns spanning \(\mathcal{R}(\mathbf{A})\), so that \(\mathbf{P}_\mathbf{A}= \mathbf{U}\mathbf{U}^\top\) and \(\mathbf{U}^\top \mathbf{U}= \mathbf{I}_r\). Hence \[\|\mathbf{P}_\mathbf{A}\mathbf{E}\|_2 = \|\mathbf{U}\mathbf{U}^\top \mathbf{E}\|_2 = \|\mathbf{U}^\top \mathbf{E}\|_2.\] Set \[\mathbf{Y}:= \mathbf{U}^\top \mathbf{E}\in \mathbb{R}^{r \times d}.\] We now bound \(\|\mathbf{Y}\|_2\). Fix \(\mathbf{x}\in \mathbb{S}^{d-1}\) and \(\mathbf{y}\in \mathbb{S}^{r-1}\), and write \(\mathbf{u}:= \mathbf{U}\mathbf{y}\in \mathbb{R}^n\), so that \(\|\mathbf{u}\|_2 = 1\). Denote the \(j\)-th column of \(\mathbf{E}\) by \(\mathbf{E}_{\cdot j}\) and define \[\xi_j := \mathbf{u}^\top \mathbf{E}_{\cdot j} = \sum_{i=1}^{n} \mathbf{u}_i \, \mathbf{E}_{ij}.\] Because the entries of \(\mathbf{E}\) are independent and centered, the variables \(\xi_1, \dots, \xi_d\) are independent and centered.

We claim that each \(\xi_j\) for \(1 \le j \le d\) is sub-gaussian. Indeed, since \(\mathbf{u}_i \mathbf{E}_{ij}\) is centered and lies in the interval \([-\rho \left|\mathbf{u}_i\right|, \rho \left|\mathbf{u}_i\right|]\), Lemma 9 gives, for every \(\lambda \in \mathbb{R}\), \[\mathbb{E}\left[\exp(\lambda \, \mathbf{u}_i \mathbf{E}_{ij})\right] \le \exp\left(\frac{\lambda^2 \rho^2 \mathbf{u}_i^2}{2}\right).\] By independence over \(i\) and \(\sum_{i=1}^{n} \mathbf{u}_i^2 = 1\), \[\mathbb{E}\left[\exp(\lambda \xi_j)\right] = \prod_{i=1}^{n} \mathbb{E}\left[\exp(\lambda \, \mathbf{u}_i \mathbf{E}_{ij})\right] \le \exp\left(\frac{\lambda^2 \rho^2}{2} \sum_{i=1}^{n} \mathbf{u}_i^2\right) = \exp\left(\frac{\lambda^2 \rho^2}{2}\right).\] Thus each \(\xi_j\) is centered and sub-gaussian with proxy \(\rho^2\).

Now \[\mathbf{y}^\top \mathbf{Y}\mathbf{x}= \mathbf{y}^\top \mathbf{U}^\top \mathbf{E}\mathbf{x}= \mathbf{u}^\top \mathbf{E}\mathbf{x}= \sum_{j=1}^{d} \mathbf{x}_j \xi_j.\] Applying the same moment-generating-function computation, now to the independent variables \(\xi_j\) with \(\sum_{j=1}^{d} \mathbf{x}_j^2 = 1\), we obtain \[\mathbb{E}\left[\exp(\lambda \mathbf{y}^\top \mathbf{Y}\mathbf{x})\right] = \prod_{j=1}^{d} \mathbb{E}\left[\exp(\lambda \mathbf{x}_j \xi_j)\right] \le \exp\left(\frac{\lambda^2 \rho^2}{2} \sum_{j=1}^{d} \mathbf{x}_j^2\right) = \exp\left(\frac{\lambda^2 \rho^2}{2}\right).\]

By Markov’s inequality, \[\begin{align} \mathbb{P}\left(\mathbf{y}^\top \mathbf{Y}\mathbf{x}> t\right) &= \mathbb{P}\left(\exp \left(\lambda \mathbf{y}^\top \mathbf{Y}\mathbf{x}\right) > \exp \left(\lambda t\right)\right) \le \frac{\mathbb{E}\left[\exp(\lambda \mathbf{y}^\top \mathbf{Y}\mathbf{x})\right]}{ \exp \left(\lambda t\right)}\le \exp \left( \frac{\lambda^2 \rho^2}{2} - \lambda t\right). \end{align}\] Optimizing in \(\lambda\) yields: \[\mathbb{P}\left(\mathbf{y}^\top \mathbf{Y}\mathbf{x}> t\right) \le \exp\left(-\frac{t^2}{2 \rho^2}\right).\] Bounding the other side of the tail symmetrically and taking a union bound, we get, for every fixed \(\mathbf{x}\in \mathbb{S}^{d-1}\), \(\mathbf{y}\in \mathbb{S}^{r-1}\) and every \(t \ge 0\), \[\label{eq:chernoff-PAE} \mathbb{P}\left(\left|\mathbf{y}^\top \mathbf{Y}\mathbf{x}\right| > t\right) \le 2 \exp\left(-\frac{t^2}{2 \rho^2}\right).\tag{21}\]

Choose \(1/4\)-nets \(\mathcal{N}_d \subset \mathbb{S}^{d-1}\) and \(\mathcal{N}_r \subset \mathbb{S}^{r-1}\) with \(\left|\mathcal{N}_d\right| \le 9^d\) and \(\left|\mathcal{N}_r\right| \le 9^r\) [28]. By Lemma 8 with \(\varepsilon = 1/4\) and a union bound applied to 21 , \[\mathbb{P}\left(\|\mathbf{Y}\|_2 > 2t\right) \le \mathbb{P}\left(\max_{\mathbf{x}\in \mathcal{N}_d, \, \mathbf{y}\in \mathcal{N}_r} \left|\mathbf{y}^\top \mathbf{Y}\mathbf{x}\right| > t\right) \le 2 \cdot 9^{d+r} \exp\left(-\frac{t^2}{2 \rho^2}\right).\] Now choose \[t := C \rho \left(\sqrt{d+r} + u\right)\] with \(C > 0\) large enough. Since \((a+b)^2 \ge a^2 + b^2\) for \(a, b \ge 0\), \[\frac{t^2}{2 \rho^2} \ge C_1 (d+r) + C_2 u^2\] for some positive constants \(C_1, C_2\). Choosing \(C_1 > \log 9\) yields \[2 \cdot 9^{d+r} \exp\left(-\frac{t^2}{2 \rho^2}\right) \le 2 e^{-c_0 u^2}\] for some absolute \(c_0 > 0\). This proves \[\mathbb{P}\left(\|\mathbf{P}_\mathbf{A}\mathbf{E}\|_2 = \|\mathbf{Y}\|_2 > 2 C \rho \left(\sqrt{d+r} + u\right)\right) \le 2 e^{-c_0 u^2}.\]


6.6 Proof of Corollary 1↩︎

We first verify that the condition \(\kappa\ge \log n\) in Theorem 2. Since \(\bar\nu\ge\nu\ge\nu_0\), for sufficiently large \(n\) \[\kappa=\max\{\mu,n\bar\nu\}\ge n\bar\nu\ge n\nu_0\ge \log n.\] Since \(r = \mathop{\mathrm{rank}}(\mathbf{A}) \le d\), we have \(\sqrt{d+r} \le \sqrt{2d}\), and Theorem 2 yields the rank-free relaxation \[\mathbb{P}\left(\sigma_d(\tilde{\mathbf{A}}) \ge \rho\left(\sqrt{n\nu} - \sqrt{\mu} - C_a\,\kappa^{1/3}(\log n)^{2/3} - C_1\sqrt{2}\left(\sqrt{d}+u\right)\right)\right) \ge 1 - n^{-a} - 2 e^{-c_1 u^2}.\] Set \(u := c_1^{-1/2}\sqrt{a\log n}\) in this inequality, so that \(2 e^{-c_1 u^2} = 2 n^{-a}\) and the total failure probability is at most \(3 n^{-a}\). Absorbing \(\sqrt{2}\) into the constant, on an event of probability at least \(1 - 3n^{-a}\), \[\label{eq:rect-step1} \sigma_d(\tilde{\mathbf{A}}) \ge \rho\left(\sqrt{n\nu} - \sqrt{\mu} - C_a\,\kappa^{1/3}(\log n)^{2/3} - \widetilde{C}_1\sqrt{d} - C_1'\sqrt{\log n}\right),\tag{22}\] where \(\widetilde{C}_1 := C_1 \sqrt{2}\) and \(C_1' := C_1 \sqrt{2} c_1^{-1/2}\sqrt{a}\).

We now verify that each subtracted term on the right-hand side is asymptotically dominated by the leading \(\sqrt{n\nu}\) term under the assumptions of the corollary.

6.6.0.1 Lower bound on the leading term.

From \(\nu \ge \nu_0\), \(\sqrt{n\nu} \ge \sqrt{n\nu_0}\).

6.6.0.2 Bound on \(\sqrt{\mu}\).

Since \(\operatorname{Var}(\mathbf{E}_{ij}) \le \rho^2/4\) by 19 , the definition \(\mu := \rho^{-2}\max_i \sum_j \operatorname{Var}(\mathbf{E}_{ij})\) gives \(\mu \le d/4\). Combined with \(d \le \eta_0 n\), this yields \(\sqrt{\mu} \le \tfrac{1}{2}\sqrt{\eta_0 n}\).

6.6.0.3 Bound on the \(\kappa^{1/3}(\log n)^{2/3}\) term.

By the same variance bound, \(\bar\nu \le 1/4\), so \(n\bar\nu \le n/4\). Similarly \(\mu \le d/4 \le n/4\). Hence \(\kappa := \max\left\{\mu, n\bar\nu\right\} \le n/4\), and \(\kappa^{1/3}(\log n)^{2/3} \le (n/4)^{1/3}(\log n)^{2/3} = o(\sqrt{n})\) as \(n \to \infty\).

6.6.0.4 Bound on \(\sqrt{d}\).

Under \(d \le \eta_0 n\), \(\sqrt{d} \le \sqrt{\eta_0 n}\).

6.6.0.5 Bound on \(\sqrt{\log n}\).

Trivially \(\sqrt{\log n} = o(\sqrt{n})\).

Combining these bounds, the right-hand side of 22 is at least \[\rho\left(\sqrt{n\nu_0} - \tfrac{1}{2}\sqrt{\eta_0 n} - \widetilde{C}_1\sqrt{\eta_0 n} - o(\sqrt{n})\right) \;=\; \rho\sqrt{n}\left(\sqrt{\nu_0} - \left(\tfrac{1}{2} + \widetilde{C}_1\right)\sqrt{\eta_0} - o(1)\right).\] Choose \(\eta_0 = \eta_0(a, \nu_0) > 0\) small enough that \(\left(\tfrac{1}{2} + \widetilde{C}_1\right)\sqrt{\eta_0} \le \tfrac{1}{2}\sqrt{\nu_0}\), and choose \(n_0 = n_0(a, \nu_0)\) large enough that the \(o(1)\) term is at most \(\tfrac{1}{4}\sqrt{\nu_0}\) for all \(n \ge n_0\). Then on the same event, \(\sigma_d(\tilde{\mathbf{A}}) \ge c_0\,\rho\sqrt{n}\) with \(c_0 := \tfrac{1}{4}\sqrt{\nu_0}\), establishing ?? .


6.7 Proof of Theorem 3↩︎

For any \(\mathbf{M}\in \mathbb{R}^{n\times d}\), \(\sum_{i=1}^{d}\sigma_i(\mathbf{M})^2 = \|\mathbf{M}\|_F^2 = \mathop{\mathrm{tr}}(\mathbf{M}^\top \mathbf{M})\). Therefore, \[\begin{align} \mathbb{E}\left[\sum_{i=1}^{d}\sigma_i(\tilde{\mathbf{A}})^2\right] &= \mathbb{E}\left[\|\tilde{\mathbf{A}}\|_F^2\right] = \mathbb{E}\left[\mathop{\mathrm{tr}}\left((\mathbf{A}+\mathbf{E})^\top(\mathbf{A}+\mathbf{E})\right)\right]\\ &= \mathop{\mathrm{tr}}(\mathbf{A}^\top \mathbf{A}) + \mathbb{E}\left[\mathop{\mathrm{tr}}(\mathbf{A}^\top \mathbf{E})\right] + \mathbb{E}\left[\mathop{\mathrm{tr}}(\mathbf{E}^\top \mathbf{A})\right] + \mathbb{E}\left[\mathop{\mathrm{tr}}(\mathbf{E}^\top \mathbf{E})\right]. \end{align}\] Since the entries of \(\mathbf{E}\) are centered, \(\mathbb{E}\left[\mathop{\mathrm{tr}}(\mathbf{A}^\top\mathbf{E})\right] = \mathbb{E}\left[\mathop{\mathrm{tr}}(\mathbf{E}^\top\mathbf{A})\right] = 0\). Hence \[\mathbb{E}\left[\sum_{i=1}^{d}\sigma_i(\tilde{\mathbf{A}})^2\right] = \mathop{\mathrm{tr}}(\mathbf{A}^\top\mathbf{A}) + \mathbb{E}\left[\mathop{\mathrm{tr}}(\mathbf{E}^\top\mathbf{E})\right] = \sum_{i=1}^{d}\sigma_i(\mathbf{A})^2 + \mathbb{E}\left[\|\mathbf{E}\|_F^2\right].\]


6.8 Proof of Lemma 2↩︎

6.8.0.1 Step 1: Top-\(k\) energy of \(\mathbf{A}\) as a contour integral.

Applied to \(\mathbf{G}\), the trace contour-integral identity 9 gives \[\label{eq:topk-A} \sum_{i=1}^{k}\sigma_i(\mathbf{A})^2 \;=\; \frac{1}{2\pi \mathrm{i}}\oint_\Gamma z\,\mathop{\mathrm{tr}}\left(\mathbf{R}_\mathbf{G}(z)\right)\,dz.\tag{23}\]

6.8.0.2 Step 2: \(\Gamma\) also separates the top-\(k\) eigenvalues of \(\tilde{\mathbf{G}}\).

By Weyl’s inequality, \(\left|\lambda_i(\tilde{\mathbf{G}}) - \lambda_i(\mathbf{G})\right| \le \|\boldsymbol{\Delta}\|_2 < g/2\) for every \(i\). Since each point of \(\Gamma\) is at distance at least \(g/2\) from \(\operatorname{spec}(\mathbf{G})\), the contour \(\Gamma\) continues to enclose exactly \(\lambda_1(\tilde{\mathbf{G}}),\dots,\lambda_k(\tilde{\mathbf{G}})\) and to exclude \(\lambda_{k+1}(\tilde{\mathbf{G}}),\dots,\lambda_d(\tilde{\mathbf{G}})\). Applying 9 to \(\tilde{\mathbf{G}}\) along the same \(\Gamma\), \[\label{eq:topk-Atil} \sum_{i=1}^{k}\sigma_i(\tilde{\mathbf{A}})^2 \;=\; \frac{1}{2\pi \mathrm{i}}\oint_\Gamma z\,\mathop{\mathrm{tr}}\left(\mathbf{R}_{\tilde{\mathbf{G}}}(z)\right)\,dz.\tag{24}\]

6.8.0.3 Step 3: Subtraction.

Subtracting 23 from 24 and using linearity of the trace, \[\label{eq:topk-shift} \sum_{i=1}^{k}\sigma_i(\tilde{\mathbf{A}})^2 - \sum_{i=1}^{k}\sigma_i(\mathbf{A})^2 \;=\; \frac{1}{2\pi \mathrm{i}}\oint_\Gamma z\,\mathop{\mathrm{tr}}\left(\mathbf{R}_{\tilde{\mathbf{G}}}(z) - \mathbf{R}_\mathbf{G}(z)\right)\,dz.\tag{25}\]

6.8.0.4 Step 4: Neumann expansion.

Combining the hypothesis \(\|\boldsymbol{\Delta}\|_2<g/2\) with the resolvent-norm bound 12 , \[\label{eq:RGD-bound} \|\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\|_2 \;\le\; \|\mathbf{R}_\mathbf{G}(z)\|_2\,\|\boldsymbol{\Delta}\|_2 \;\le\; \tfrac{2}{g}\,\|\boldsymbol{\Delta}\|_2 \;<\; 1, \qquad z\in\Gamma.\tag{26}\] Hence the Neumann series 7 converges and \[\label{eq:resolvent-diff} \mathbf{R}_{\tilde{\mathbf{G}}}(z) - \mathbf{R}_\mathbf{G}(z) \;=\; \sum_{\ell=1}^{\infty}\left(\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\right)^\ell \mathbf{R}_\mathbf{G}(z).\tag{27}\]

6.8.0.5 Step 5: Interchange integral and series

Substituting the Neumann series 7 into 25 yields: \[\begin{align} \sum_{i=1}^{k}\sigma_i(\tilde{\mathbf{A}})^2 - \sum_{i=1}^{k}\sigma_i(\mathbf{A})^2 \;&=\; \frac{1}{2\pi \mathrm{i}}\oint_\Gamma z\,\mathop{\mathrm{tr}}\left(\sum_{\ell=1}^{\infty}\left(\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\right)^\ell \mathbf{R}_\mathbf{G}(z)\right)\,dz \\ &= \; \frac{1}{2\pi \mathrm{i}}\oint_\Gamma \sum_{\ell=1}^{\infty} \,z\,\mathop{\mathrm{tr}}\left(\left(\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\right)^\ell \mathbf{R}_\mathbf{G}(z)\right)\,dz. \end{align}\]

Denote \(f_m(z) := \sum_{\ell = 1}^m z\,\mathop{\mathrm{tr}}\left(\left(\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\right)^\ell \mathbf{R}_\mathbf{G}(z)\right).\) Using \(\left|\mathop{\mathrm{tr}}(\mathbf{M})\right|\le d\|\mathbf{M}\|_2\), we get: \[\begin{align} \left|f_m(z)\right| \le \left|z\right| d \sum_{\ell = 1}^m \|\left(\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\right)^\ell \mathbf{R}_\mathbf{G}(z)\|_2. \end{align}\]

By 12 and triangle inequality, \[\left|f_m(z)\right| \le \left|z\right|\, d \, \frac{2}{g} \sum_{\ell = 1}^m \|\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\|_2 ^\ell,\]

Since \(\|\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\|_2 < 1\), the geometric series is finite: \[\sum_{\ell = 1}^m \|\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\|_2 ^\ell \leq \sum_{\ell = 1}^\infty \|\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\|_2 ^\ell = \frac{\|\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\|_2}{1-\|\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\|_2}.\]

Hence, \[\left|f(z)\right| \le h(z) := \left|z\right|\, d \, \frac{2}{g} \frac{\|\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\|_2}{1-\|\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\|_2}.\]

By the boundedness of \(\Gamma\), \(h(z)\) is clearly integrable. The Lebesgue’s dominated convergence theorem then gives: \[\sum_{i=1}^{k}\sigma_i(\tilde{\mathbf{A}})^2 - \sum_{i=1}^{k}\sigma_i(\mathbf{A})^2 \;= \; \frac{1}{2\pi \mathrm{i}}\sum_{\ell=1}^{\infty}\oint_\Gamma \,z\,\mathop{\mathrm{tr}}\left(\left(\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\right)^\ell \mathbf{R}_\mathbf{G}(z)\right)\,dz.\]

Finally, taking expectation yields: \[\mathbb{E}\left[\sum_{i=1}^{k}\sigma_i(\tilde{\mathbf{A}})^2\right] - \sum_{i=1}^{k}\sigma_i(\mathbf{A})^2 \;=\; \sum_{\ell=1}^{\infty} I_\ell, \qquad I_\ell = \frac{1}{2\pi \mathrm{i}}\oint_\Gamma z\,\mathbb{E}\left[\mathop{\mathrm{tr}}\left(\left(\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\right)^\ell \mathbf{R}_\mathbf{G}(z)\right)\right]\,dz.\]


6.9 Proof of Lemma 3↩︎

By linearity and \(\mathbb{E}\left[\boldsymbol{\Delta}\right]=\mathbf{D}\) (Section 2.2), \[I_1 = \frac{1}{2\pi \mathrm{i}}\oint_\Gamma z\,\mathop{\mathrm{tr}}\left(\mathbf{R}_\mathbf{G}(z)\mathbf{D}\mathbf{R}_\mathbf{G}(z)\right)\,dz = \mathop{\mathrm{tr}}\left(\frac{1}{2\pi \mathrm{i}}\oint_\Gamma z\,\mathbf{R}_\mathbf{G}(z)^2\,dz \;\mathbf{D}\right).\] Writing \(\mathbf{V}= [\mathbf{v}_1,\dots,\mathbf{v}_d]\), \[\mathbf{R}_\mathbf{G}(z)^2 = \mathbf{V}\mathop{\mathrm{diag}}\left((z-\lambda_1)^{-2},\dots,(z-\lambda_d)^{-2}\right)\mathbf{V}^\top.\] Since \[\operatorname*{Res}_{z=\lambda_i}\frac{z}{(z-\lambda_i)^2} = \left.\tfrac{d}{dz}z\right|_{z=\lambda_i}=1,\] Cauchy’s integral formula gives \[\frac{1}{2\pi \mathrm{i}}\oint_\Gamma \frac{z}{(z-\lambda_i)^2}dz = \mathbb{1}_{i\le k}.\] Therefore \[\frac{1}{2\pi \mathrm{i}}\oint_\Gamma z\,\mathbf{R}_\mathbf{G}(z)^2 dz = \mathbf{V}_k\mathbf{V}_k^\top,\] with \[\mathbf{V}_k=[\mathbf{v}_1,\dots,\mathbf{v}_k],\] and \[I_1 = \mathop{\mathrm{tr}}\left(\mathbf{V}_k\mathbf{V}_k^\top\mathbf{D}\right) = \mathop{\mathrm{tr}}\left(\mathbf{V}_k^\top\mathbf{D}\mathbf{V}_k\right).\] Ky Fan’s maximum principle [44] gives \[\mathop{\mathrm{tr}}\left(\mathbf{V}_k^\top\mathbf{D}\mathbf{V}_k\right) \le \max_{\mathbf{Q}\in\mathbb{R}^{d\times k}\,:\,\mathbf{Q}^\top\mathbf{Q}=\mathbf{I}_k}\mathop{\mathrm{tr}}\left(\mathbf{Q}^\top\mathbf{D}\mathbf{Q}\right) = \sum_{j=1}^{k}\nu_j^\downarrow,\] where \(\nu_1^\downarrow \ge \cdots \ge \nu_d^\downarrow\) is the decreasing rearrangement of \(\nu_1,\dots,\nu_d\).


6.10 Proof of Lemma 4↩︎

6.10.0.1 Step 1: Residue calculus.

Using cyclicity and linearity of trace plus the spectral decomposition 6 \(\mathbf{R}_\mathbf{G}(z) = \sum_{i=1}^d\frac{\mathbf{v}_i\mathbf{v}_i^\top}{z-\lambda_i}\), \[\begin{align} &\mathop{\mathrm{tr}}\left(\mathbf{R}_\mathbf{G}(z) \boldsymbol{\Delta}\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\mathbf{R}_\mathbf{G}(z)\right)\\ &= \mathop{\mathrm{tr}}\left(\mathbf{R}_\mathbf{G}(z)^2 \boldsymbol{\Delta}\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\right)\\ &= \mathop{\mathrm{tr}}\left(\sum_{i,j=1}^{d} \frac{\mathbf{v}_i \mathbf{v}_i^\top\boldsymbol{\Delta}\mathbf{v}_j \mathbf{v}_j^\top\boldsymbol{\Delta}}{(z-\lambda_i)^2(z-\lambda_j)}\right)\\ &=\sum_{i,j=1}^{d} \mathop{\mathrm{tr}}\left( \frac{\mathbf{v}_i \mathbf{v}_i^\top\boldsymbol{\Delta}\mathbf{v}_j \mathbf{v}_j^\top\boldsymbol{\Delta}}{(z-\lambda_i)^2(z-\lambda_j)}\right)\\ &= \sum_{i,j=1}^{d}\frac{(\mathbf{v}_i^\top\boldsymbol{\Delta}\mathbf{v}_j)(\mathbf{v}_j^\top\boldsymbol{\Delta}\mathbf{v}_i)}{(z-\lambda_i)^2(z-\lambda_j)}. \end{align}\] Since \(\boldsymbol{\Delta}\) is symmetric, \((\mathbf{v}_i^\top\boldsymbol{\Delta}\mathbf{v}_j)(\mathbf{v}_j^\top\boldsymbol{\Delta}\mathbf{v}_i) = (\mathbf{v}_i^\top\boldsymbol{\Delta}\mathbf{v}_j)^2\), so \[\label{eq:I2-master} I_2 = \sum_{i,j=1}^{d}\mathbb{E}\left[(\mathbf{v}_i^\top\boldsymbol{\Delta}\mathbf{v}_j)^2\right]\frac{1}{2\pi \mathrm{i}}\oint_\Gamma\frac{z}{(z-\lambda_i)^2(z-\lambda_j)}\,dz.\tag{28}\]

We now compute the contour integral case by case.

6.10.0.2 Case 1: \(i,j\le k.\)

Clearly both poles lie inside \(\Gamma\).

If \(i\neq j\), the residues are \[\operatorname*{Res}_{z=\lambda_i}\frac{z}{(z-\lambda_i)^2(z-\lambda_j)} = \left.\frac{d}{dz}\frac{z}{z-\lambda_j}\right|_{z=\lambda_i} = -\frac{\lambda_j}{(\lambda_i-\lambda_j)^2},\] and \[\operatorname*{Res}_{z=\lambda_j}\frac{z}{(z-\lambda_i)^2(z-\lambda_j)} = \left.\frac{z}{(z-\lambda_i)^2}\right|_{z=\lambda_j} = \frac{\lambda_j}{(\lambda_i-\lambda_j)^2},\] so the integral is zero.

If \(i=j\le k\), then there is a triple pole at \(z=\lambda_i\), and \[\operatorname*{Res}_{z=\lambda_i}\frac{z}{(z-\lambda_i)^3} = \frac{1}{2!}\left.\frac{d^2}{dz^2} z\right|_{z=\lambda_i} = 0,\] so the integral is again zero.

6.10.0.3 Case 2: \(i,j>k.\)

Both poles lie outside \(\Gamma\), so the integral is zero.

6.10.0.4 Case 3: \(i\le k<j.\)

Only the double pole at \(z=\lambda_i\) lies inside \(\Gamma\), and \[\frac{1}{2\pi \mathrm{i}}\oint_\Gamma \frac{z}{(z-\lambda_i)^2(z-\lambda_j)}\,dz = \operatorname*{Res}_{z=\lambda_i}\frac{z}{(z-\lambda_i)^2(z-\lambda_j)} = \left.\frac{d}{dz}\frac{z}{z-\lambda_j}\right|_{z=\lambda_i} = -\frac{\lambda_j}{(\lambda_i-\lambda_j)^2}.\]

6.10.0.5 Case 4: \(j\le k<i.\)

Only the simple pole at \(z=\lambda_j\) lies inside \(\Gamma\), and \[\frac{1}{2\pi \mathrm{i}}\oint_\Gamma \frac{z}{(z-\lambda_i)^2(z-\lambda_j)}\,dz = \operatorname*{Res}_{z=\lambda_j}\frac{z}{(z-\lambda_i)^2(z-\lambda_j)} = \left.\frac{z}{(z-\lambda_i)^2}\right|_{z=\lambda_j} = \frac{\lambda_j}{(\lambda_i-\lambda_j)^2}.\]

Therefore \[\begin{align} I_2 &= \sum_{i=1}^k \sum_{j=k+1}^d \mathbb{E}\left[(\mathbf{v}_i^\top\boldsymbol{\Delta}\mathbf{v}_j)^2\right] \left(-\frac{\lambda_j}{(\lambda_i-\lambda_j)^2}\right) \\ &\quad+ \sum_{j=1}^k \sum_{i=k+1}^d \mathbb{E}\left[(\mathbf{v}_i^\top\boldsymbol{\Delta}\mathbf{v}_j)^2\right] \frac{\lambda_j}{(\lambda_i-\lambda_j)^2}. \end{align}\] Relabeling indices in the second sum and using \((\mathbf{v}_i^\top\boldsymbol{\Delta}\mathbf{v}_j)^2=(\mathbf{v}_j^\top\boldsymbol{\Delta}\mathbf{v}_i)^2\), we obtain \[\label{eq:I2-clean} I_2 = \sum_{i=1}^k \sum_{j=k+1}^d \frac{\mathbb{E}\left[(\mathbf{v}_i^\top\boldsymbol{\Delta}\mathbf{v}_j)^2\right]}{\lambda_i-\lambda_j}.\tag{29}\]

6.10.0.6 Step 2: Bounding \(\mathbb{E}\left[(\mathbf{v}_i^\top\boldsymbol{\Delta}\mathbf{v}_j)^2\right]\).

For \(i\ne j\), \(\boldsymbol{\Delta}= \mathbf{A}^\top\mathbf{E}+ \mathbf{E}^\top\mathbf{A}+ \mathbf{E}^\top\mathbf{E}\) and \(\mathbf{A}\mathbf{v}_\ell = \sigma_\ell\mathbf{u}_\ell\) give \[\mathbf{v}_i^\top\boldsymbol{\Delta}\mathbf{v}_j = \sigma_i\,\mathbf{u}_i^\top\mathbf{E}\mathbf{v}_j + \sigma_j\,\mathbf{u}_j^\top\mathbf{E}\mathbf{v}_i + \mathbf{v}_i^\top\mathbf{E}^\top\mathbf{E}\mathbf{v}_j.\] Applying the scalar inequality \((a+b+c)^2 \le 3(a^2+b^2+c^2)\) yields:8 \[\label{eq:vidv-split} \mathbb{E}\left[(\mathbf{v}_i^\top\boldsymbol{\Delta}\mathbf{v}_j)^2\right] \le 3\sigma_i^2\,\mathbb{E}\left[(\mathbf{u}_i^\top\mathbf{E}\mathbf{v}_j)^2\right] + 3\sigma_j^2\,\mathbb{E}\left[(\mathbf{u}_j^\top\mathbf{E}\mathbf{v}_i)^2\right] + 3\,\mathbb{E}\left[(\mathbf{v}_i^\top\mathbf{E}^\top\mathbf{E}\mathbf{v}_j)^2\right].\tag{30}\]

6.10.0.7 Linear terms.

For unit vectors \(\mathbf{u}\in\mathbb{R}^n\), \(\mathbf{v}\in\mathbb{R}^d\), notice that: \[\mathbf{u}^\top\mathbf{E}\mathbf{v}= \sum_{a=1}^n \sum_{b=1}^d \mathbf{u}_a\mathbf{v}_b\mathbf{E}_{ab},\] which is a sum of independent centered sub-gaussian random variables.

By 20 , \[\|\mathbf{u}_a\mathbf{v}_b\mathbf{E}_{ab}\|_{\psi_2}\le C\rho\left|\mathbf{u}_a\mathbf{v}_b\right|.\] The general Hoeffding inequality [28] gives, for some absolute constant \(c\), \[\mathbb{P}\left(\left|\mathbf{u}^\top\mathbf{E}\mathbf{v}\right|\ge t\right) \le 2\exp\left(-\frac{ct^2}{\rho^2}\right),\] since \(\sum_{a,b}\mathbf{u}_a^2\mathbf{v}_b^2 = 1\). Therefore, for some absolute constant \(C,\) \[\mathbb{E}\left[(\mathbf{u}^\top\mathbf{E}\mathbf{v})^2\right] = \int_0^\infty 2t\,\mathbb{P}\left(\left|\mathbf{u}^\top\mathbf{E}\mathbf{v}\right|\ge t\right)\,dt \le 4\int_0^\infty t e^{-ct^2/\rho^2}\,dt = \frac{2\rho^2}{c} \le C\rho^2.\]

6.10.0.8 Quadratic term.

Write \(\mathbf{x}=\mathop{\mathrm{vec}}(\mathbf{E})\in\mathbb{R}^{nd}\). Then \[\mathbf{v}_i^\top\mathbf{E}^\top\mathbf{E}\mathbf{v}_j = \mathbf{x}^\top\mathbf{B}^{ij}\mathbf{x}, \qquad \mathbf{B}^{ij} := \tfrac{1}{2}\bigl((\mathbf{v}_i\mathbf{v}_j^\top)\otimes \mathbf{I}_n + (\mathbf{v}_j\mathbf{v}_i^\top)\otimes \mathbf{I}_n\bigr).\] Indeed, by the cyclic property of trace:

\[\mathbf{v}_i^\top\mathbf{E}^\top\mathbf{E}\mathbf{v}_j = \mathop{\mathrm{tr}}(\mathbf{E}^\top\mathbf{E}\mathbf{v}_j\mathbf{v}_i^\top).\]

Recall  [45], \[\begin{align} \mathop{\mathrm{vec}}(\mathbf{A}\mathbf{X}\mathbf{B}) &= (\mathbf{B}^\top\otimes \mathbf{A})\,\mathop{\mathrm{vec}}(\mathbf{X}), \tag{31}\\ \mathop{\mathrm{tr}}(\mathbf{A}^\top\mathbf{B}) &= \mathop{\mathrm{vec}}(\mathbf{A})^\top\mathop{\mathrm{vec}}(\mathbf{B}), \tag{32} \end{align}\]

First, use 32 with \(\mathbf{A}= \mathbf{E}, \mathbf{B}= \mathbf{E}\mathbf{v}_j \mathbf{v}_i^\top\): \[\mathop{\mathrm{tr}}(\mathbf{E}^\top\mathbf{E}\mathbf{v}_j\mathbf{v}_i^\top) = \mathop{\mathrm{vec}}(\mathbf{E})^\top\mathop{\mathrm{vec}}( \mathbf{E}\mathbf{v}_j \mathbf{v}_i^\top).\]

Second, use 31 with \(\mathbf{A}= \mathbf{I}_n, \mathbf{X}= \mathbf{E}, \mathbf{B}= \mathbf{v}_j \mathbf{v}_i^\top\): \[\mathop{\mathrm{vec}}( \mathbf{E}\mathbf{v}_j \mathbf{v}_i^\top) = \bigl((\mathbf{v}_i\mathbf{v}_j^\top)\otimes\mathbf{I}_n\bigr)\mathop{\mathrm{vec}}(\mathbf{E}).\]

Substituting the above two results, we get: \[\mathbf{v}_i^\top\mathbf{E}^\top\mathbf{E}\mathbf{v}_j =\mathop{\mathrm{vec}}(\mathbf{E})^\top\bigl((\mathbf{v}_i\mathbf{v}_j^\top)\otimes\mathbf{I}_n\bigr)\mathop{\mathrm{vec}}(\mathbf{E}).\]

Finally, since a quadratic form depends only on the symmetric part of its matrix: \[\begin{align} \mathbf{v}_i^\top\mathbf{E}^\top\mathbf{E}\mathbf{v}_j &= \frac{1}{2}\mathop{\mathrm{vec}}(\mathbf{E})^\top\left((\mathbf{v}_i\mathbf{v}_j^\top)\otimes\mathbf{I}_n + \left((\mathbf{v}_i\mathbf{v}_j^\top)\otimes\mathbf{I}_n\right)^\top\right)\mathop{\mathrm{vec}}(\mathbf{E})\\ &= \frac{1}{2}\mathop{\mathrm{vec}}(\mathbf{E})^\top\left((\mathbf{v}_i\mathbf{v}_j^\top)\otimes\mathbf{I}_n + (\mathbf{v}_j\mathbf{v}_i^\top)\otimes\mathbf{I}_n\right)\mathop{\mathrm{vec}}(\mathbf{E})\\ &= \mathop{\mathrm{vec}}(\mathbf{E})^\top\mathbf{B}^{ij} \mathop{\mathrm{vec}}(\mathbf{E}). \end{align}\]

Next, notice that since \(\mathbf{v}_i\) and \(\mathbf{v}_j\) are unit vectors, both \((\mathbf{v}_i\mathbf{v}_j^\top)\otimes\mathbf{I}_n\) and \((\mathbf{v}_j\mathbf{v}_i^\top)\otimes\mathbf{I}_n\) have \(n\) \(1\)-valued singular values and \(n(d-1)\) \(0\)-valued singular values. Therefore, \[\begin{align} &\|(\mathbf{v}_i\mathbf{v}_j^\top)\otimes\mathbf{I}_n\|_2 = \|(\mathbf{v}_j\mathbf{v}_i^\top)\otimes\mathbf{I}_n\|_2 =1 \\ &\|(\mathbf{v}_i\mathbf{v}_j^\top)\otimes\mathbf{I}_n\|_F = \|(\mathbf{v}_j\mathbf{v}_i^\top)\otimes\mathbf{I}_n\|_F = \sqrt{n}. \end{align}\]

Triangle inequality yields: \[\|\mathbf{B}^{ij}\|_2 \le 1, \quad\|\mathbf{B}^{ij}\|_F \le \sqrt{n}.\]

The Hanson–Wright inequality [28] therefore gives, for any \(t\ge 0\), there exists a positive constant \(c\) such that: \[Y := \mathbf{v}_i^\top\mathbf{E}^\top\mathbf{E}\mathbf{v}_j - \mathbb{E}\left[\mathbf{v}_i^\top\mathbf{E}^\top\mathbf{E}\mathbf{v}_j\right], \quad \mathbb{P}\left(\left|Y\right|\ge t\right) \le 2\exp\left[-c\min\left(\frac{t^2}{\rho^4 n},\,\frac{t}{\rho^2}\right)\right].\]

We can then calculate: \[\mathbb{E}[Y^2] = \int_0^\infty 2t\,\mathbb{P}(|Y|\ge t)\,dt.\]

The two arguments of the minimum are equal at \(t=\rho^2 n\), so \[\mathbb{E}\left[Y^2\right] \le 4\int_0^{\rho^2 n} t e^{-ct^2/(\rho^4 n)}\,dt + 4\int_{\rho^2 n}^{\infty} t e^{-ct/\rho^2}\,dt.\]

For the first term, with the substitution \[u=\frac{ct^2}{\rho^4 n}, \qquad du=\frac{2ct}{\rho^4 n}\,dt,\] we get \[\int_0^{\rho^2 n} t e^{-ct^2/(\rho^4 n)}dt = \frac{\rho^4 n}{2c}\int_0^{cn} e^{-u}\,du \le \frac{\rho^4 n}{2c}.\] For the second term, with \[s=\frac{ct}{\rho^2}, \qquad ds = \frac{c}{\rho^2}\,dt,\] we obtain \[\int_{\rho^2 n}^{\infty} t e^{-ct/\rho^2}dt = \frac{\rho^4}{c^2}\int_{cn}^{\infty} s e^{-s}\,ds = \frac{\rho^4}{c^2}(cn+1)e^{-cn}.\] Since \(cn\ge 0\) and \(x+1\le e^x\) for all \(x\ge 0\), we have \[(cn+1)e^{-cn}\le 1,\] and therefore \[\frac{\rho^4}{c^2}(cn+1)e^{-cn} \le C\rho^4.\]

Adding the two integrals, we get: \[\operatorname{Var}[\mathbf{v}_i^\top\mathbf{E}^\top\mathbf{E}\mathbf{v}_j] = \mathbb{E}\left[\left(\mathbf{v}_i^\top\mathbf{E}^\top\mathbf{E}\mathbf{v}_j - \mathbb{E}\left[\mathbf{v}_i^\top\mathbf{E}^\top\mathbf{E}\mathbf{v}_j\right]\right)^2\right] \le C \rho^4 n.\]

Since \(\mathbb{E}\left[\mathbf{E}^\top\mathbf{E}\right] = \mathbf{D}\), \(\mathbb{E}\left[\mathbf{v}_i^\top\mathbf{E}^\top\mathbf{E}\mathbf{v}_j\right] = \mathbf{v}_i^\top\mathbf{D}\mathbf{v}_j\).

Finally, the decomposition of the variance gives:

\[\begin{align} \mathbb{E}\left[(\mathbf{v}_i^\top\mathbf{E}^\top\mathbf{E}\mathbf{v}_j)^2\right] &= \mathbb{E}\left[\mathbf{v}_i^\top\mathbf{E}^\top\mathbf{E}\mathbf{v}_j\right]^2 + \operatorname{Var}[\mathbf{v}_i^\top\mathbf{E}^\top\mathbf{E}\mathbf{v}_j]\\ &\le (\mathbf{v}_i^\top\mathbf{D}\mathbf{v}_j)^2 + C\rho^4 n. \end{align}\]

Combining the linear and quadratic parts in 30 , \[\label{eq:vidv-final} \mathbb{E}\left[(\mathbf{v}_i^\top\boldsymbol{\Delta}\mathbf{v}_j)^2\right] \le C\left[\rho^2(\sigma_i^2+\sigma_j^2) + (\mathbf{v}_i^\top\mathbf{D}\mathbf{v}_j)^2 + \rho^4 n\right].\tag{33}\]

6.10.0.9 Step 3: Plugging into \(I_2\).

Substituting 33 into 29 and using \(\lambda_i=\sigma_i^2\), \[I_2 \le C\sum_{i=1}^{k}\sum_{j=k+1}^{d}\frac{\rho^2(\lambda_i+\lambda_j) + (\mathbf{v}_i^\top\mathbf{D}\mathbf{v}_j)^2 + \rho^4 n}{\lambda_i-\lambda_j}.\] For the middle term, with \(\mathbf{P}:= \mathbf{V}_k\mathbf{V}_k^\top\) and \(\mathbf{V}_{k,\perp}=[\mathbf{v}_{k+1},\dots,\mathbf{v}_d]\), \(\mathbf{V}_{k,\perp}\mathbf{V}_{k,\perp}^\top= \mathbf{I}- \mathbf{P}\), we can express the double sum in matrix form:

\[\begin{align} &\sum_{i=1}^{k}\sum_{j=k+1}^{d}(\mathbf{v}_i^\top\mathbf{D}\mathbf{v}_j)^2 = \|\mathbf{V}_k^\top\mathbf{D}\mathbf{V}_{k,\perp}\|_F^2 \\ &= \mathop{\mathrm{tr}}\left( \mathbf{V}_k^\top\mathbf{D}\mathbf{V}_{k,\perp} \mathbf{V}_{k,\perp}^\top\mathbf{D}\mathbf{V}_k \right) = \mathop{\mathrm{tr}}\left(\mathbf{V}_k^\top\mathbf{D}(\mathbf{I}-\mathbf{P}) \mathbf{D}\mathbf{V}_k\right)\\ &= \mathop{\mathrm{tr}}\left(\mathbf{P}\mathbf{D}(\mathbf{I}-\mathbf{P}) \mathbf{D}\right)\\ &= \mathop{\mathrm{tr}}(\mathbf{P}\mathbf{D}^2) - \mathop{\mathrm{tr}}(\mathbf{P}\mathbf{D}\mathbf{P}\mathbf{D}) \le \mathop{\mathrm{tr}}(\mathbf{P}\mathbf{D}^2), \end{align}\]

since \(\mathop{\mathrm{tr}}(\mathbf{P}\mathbf{D}\mathbf{P}\mathbf{D}) = \mathop{\mathrm{tr}}(\mathbf{D}^{1/2} \mathbf{P}\mathbf{D}^{1/2} \mathbf{D}^{1/2} \mathbf{P}\mathbf{D}^{1/2}) = \|\mathbf{D}^{1/2}\mathbf{P}\mathbf{D}^{1/2}\|_F^2 \ge 0\).

Since \(\mathbf{D}^2 = \mathop{\mathrm{diag}}(\nu_1^2,\dots,\nu_d^2)\), Ky Fan’s maximum principle yields \(\mathop{\mathrm{tr}}(\mathbf{P}\mathbf{D}^2) \le \sum_{j=1}^{k}(\nu_j^\downarrow)^2\).

Combining with the spectral-gap bound \(\lambda_i-\lambda_j \ge g\) for \(i\le k<j\), \[I_2 \le \frac{C}{g}\left[\rho^2\sum_{i=1}^{k}\sum_{j=k+1}^{d}(\lambda_i+\lambda_j) + \sum_{j=1}^{k}(\nu_j^\downarrow)^2 + \rho^4 n\,k(d-k)\right].\]


6.11 Proof of Lemma 5↩︎

We will use Hoeffding’s inequality for sums of bounded independent random variables.

Theorem 8 (Hoeffding’s inequality, [28] Theorem 2.2.6). Let \(X_1,\dots,X_N\) be independent random variables with \(m_i \le X_i \le M_i\) for each \(i\). Then for every \(t > 0\), \[\mathbb{P}\left(\left|\sum_{i=1}^{N}\left(X_i - \mathbb{E}\left[X_i\right]\right)\right| \ge t\right) \le 2\exp\left(-\frac{2 t^2}{\sum_{i=1}^{N}(M_i - m_i)^2}\right).\]

[Proof of Lemma 5] The proof is based on an \(\varepsilon\)-net argument; see Section 4.2 of [28] for a review.

6.11.0.1 Step 1: Approximation.

Fix \(\varepsilon = 1/4\). By [28], there exist \(\varepsilon\)-nets \(\mathcal{N}\subset \mathbb{S}^{n-1}\) and \(\mathcal{D}\subset \mathbb{S}^{d-1}\) with bounded cardinality: \[\left|\mathcal{N}\right| \le 9^n, \qquad \left|\mathcal{D}\right| \le 9^d.\] By Lemma 8, \(\|\mathbf{E}\|_2\) can be bounded using the nets as \[\label{eq:E-net-relaxation} \|\mathbf{E}\|_2 \le 2 \max_{\mathbf{x}\in\mathcal{N},\;\mathbf{y}\in\mathcal{D}} |\mathbf{x}^\top \mathbf{E}\mathbf{y}|.\tag{34}\]

6.11.0.2 Step 2: Concentration.

Fix \(\mathbf{x}\in \mathcal{N}\) and \(\mathbf{y}\in \mathcal{D}\). Then \[\mathbf{x}^\top \mathbf{E}\mathbf{y} = \sum_{i=1}^{n}\sum_{j=1}^{d} \mathbf{x}_i\mathbf{y}_j\mathbf{E}_{ij},\] which is a sum of independent bounded random variables. Using the uniformly dithered quantizer of Section 6.2, assuming \(\mathbf{A}_{ij}\in [q_{ij}, q_{ij}+\rho]\), \[\mathbf{E}_{ij} = \begin{cases} \rho(1 - p_{ij}) & \text{with probability } p_{ij},\\ -\rho\,p_{ij} & \text{with probability } 1 - p_{ij}, \end{cases}\] where \(p_{ij}:=(\mathbf{A}_{ij}-q_{ij})/\rho\). Therefore \(\mathbf{E}_{ij}\) lies in the interval \([-\rho p_{ij},\,\rho(1-p_{ij})]\) of width \(\rho\). Applying Theorem 8, \[\begin{align} \mathbb{P}\left(\left|\mathbf{x}^\top \mathbf{E}\mathbf{y}\right| \ge t\right) &\le 2\exp\left(-\frac{2 t^2}{\sum_{i,j}\rho^2 \mathbf{x}_i^2 \mathbf{y}_j^2}\right)\\ &= 2\exp\left(-\frac{2 t^2}{\rho^2 \left(\sum_{i=1}^{n}\mathbf{x}_i^2\right)\left(\sum_{j=1}^{d}\mathbf{y}_j^2\right)}\right) = 2\exp\left(-\frac{2 t^2}{\rho^2}\right). \end{align}\]

6.11.0.3 Step 3: Union bound.

By a union bound on the event \(\left\{\max_{\mathbf{x}\in\mathcal{N},\,\mathbf{y}\in\mathcal{D}}\left|\mathbf{x}^\top \mathbf{E}\mathbf{y}\right| \ge t\right\} = \bigcup_{\mathbf{x}\in\mathcal{N},\,\mathbf{y}\in\mathcal{D}}\left\{\left|\mathbf{x}^\top \mathbf{E}\mathbf{y}\right| \ge t\right\}\), \[\begin{align} \mathbb{P}\left(\max_{\mathbf{x}\in\mathcal{N},\,\mathbf{y}\in\mathcal{D}}\left|\mathbf{x}^\top \mathbf{E}\mathbf{y}\right| \ge t\right) &\le \sum_{\mathbf{x}\in\mathcal{N},\,\mathbf{y}\in\mathcal{D}}\mathbb{P}\left(\left|\mathbf{x}^\top \mathbf{E}\mathbf{y}\right| \ge t\right)\\ &\le 9^{n+d} \cdot 2\exp\left(-\frac{2 t^2}{\rho^2}\right). \end{align}\] Now choose \[t=C_0\rho(\sqrt{n+d}+s),\] where \(C_0>0\) is a sufficiently large constant. Since \[(\sqrt{n+d}+s)^2\ge n+d+s^2,\] we obtain \[2\cdot 9^{n+d} \exp\left(-\frac{2t^2}{\rho^2}\right) \le 2e^{-cs^2}.\] for a positive constant \(c\). Hence \[\mathbb{P}\left\{ \|\mathbf{E}\|_2>C_0\rho(\sqrt{n+d}+s) \right\} \le 2e^{-cs^2}.\] Finally, taking \(s=C'\sqrt{n+d}\) with \(C'\) sufficiently large gives \[\mathbb{P}\left\{ \|\mathbf{E}\|_2\le C'\rho\sqrt{n+d} \right\} \ge 1-2e^{-(n+d)}.\] As \(n\ge d\), we have for some positive constant \(C\) \[\mathbb{P}\left\{ \|\mathbf{E}\|_2\le C\rho\sqrt n \right\} \ge 1-2e^{-(n+d)}.\]


6.12 Proof of Lemma 6↩︎

6.12.0.1 Step 1: Triangle inequality.

Recall \(\boldsymbol{\Delta}= \mathbf{A}^\top\mathbf{E}+ \mathbf{E}^\top\mathbf{A}+ \mathbf{E}^\top\mathbf{E}\). Since \(\|\mathbf{M}^\top\|_2 = \|\mathbf{M}\|_2\) for any \(\mathbf{M}\), the triangle inequality gives \[\label{eq:Delta-triangle} \|\boldsymbol{\Delta}\|_2 \le \|\mathbf{A}^\top\mathbf{E}\|_2 + \|\mathbf{E}^\top\mathbf{A}\|_2 + \|\mathbf{E}^\top\mathbf{E}\|_2 = 2\,\|\mathbf{A}^\top\mathbf{E}\|_2 + \|\mathbf{E}\|_2^2.\tag{35}\]

6.12.0.2 Step 2: Bound on \(\|\mathbf{A}^\top\mathbf{E}\|_2\) via Lemma 1.

Using the thin SVD \(\mathbf{A}= \mathbf{U}\mathbf{\Sigma}\mathbf{V}^\top,\) \[\|\mathbf{A}^\top\mathbf{E}\|_2 = \|\mathbf{V}\mathbf{\Sigma}\mathbf{U}^\top\mathbf{E}\|_2 \le \|\mathbf{V}\|_2\,\|\mathbf{\Sigma}\|_2\,\|\mathbf{U}^\top\mathbf{E}\|_2 = \sigma_1(\mathbf{A})\,\|\mathbf{U}^\top\mathbf{E}\|_2,\] Since \(\mathbf{P}_\mathbf{A}= \mathbf{U}\mathbf{U}^\top\), \[\|\mathbf{U}^\top\mathbf{E}\|_2 = \|\mathbf{U}\mathbf{U}^\top\mathbf{E}\|_2 = \|\mathbf{P}_\mathbf{A}\mathbf{E}\|_2.\] Since \(\mathbf{A}\) has full column rank, \(r := \mathop{\mathrm{rank}}(\mathbf{A}) = d\). Applying Lemma 1 with \(u = \sqrt{d}\) and \(r = d\) gives \[\mathbb{P}\left(\|\mathbf{P}_\mathbf{A}\mathbf{E}\|_2 > C_0\,\rho\left(\sqrt{2d} + \sqrt{d}\right)\right) \le 2 e^{-c_0 d},\] so on an event \(\mathcal{E}_1\) with \(\mathbb{P}\left(\mathcal{E}_1\right) \ge 1 - 2 e^{-c_0 d}\), \[\label{eq:AE-bound} \|\mathbf{A}^\top\mathbf{E}\|_2 \le C_1\,\sigma_1(\mathbf{A})\,\rho\sqrt{d}, \qquad C_1 := C_0(1 + \sqrt{2}).\tag{36}\]

6.12.0.3 Step 3: Bound on \(\|\mathbf{E}\|_2^2\) via Lemma 5.

By Lemma 5, on an event \(\mathcal{E}_2\) with \(\mathbb{P}\left(\mathcal{E}_2\right) \ge 1 - 2 e^{-(n+d)}\), \[\label{eq:E-square-bound} \|\mathbf{E}\|_2^2 \le C_2\,\rho^2 n,\tag{37}\] where \(C_2\) is the square of the constant from Lemma 5.

6.12.0.4 Step 4: Union bound.

On \(\mathcal{E}_1 \cap \mathcal{E}_2\), combining 3536 , and 37 , \[\|\boldsymbol{\Delta}\|_2 \le 2 C_1\,\sigma_1(\mathbf{A})\,\rho\sqrt{d} + C_2\,\rho^2 n \le C\left(\sigma_1(\mathbf{A})\,\rho\sqrt{d} + \rho^2 n\right),\] with \(C := \max\{2 C_1, C_2\}\). By a union bound, \[\mathbb{P}\left(\mathcal{E}_1 \cap \mathcal{E}_2\right) \ge 1 - 2 e^{-c_0 d} - 2 e^{-(n+d)} \ge 1 - 4 e^{-c_0 d},\] where the last step uses \(n+d \ge d\), so \(e^{-(n+d)} \le e^{-d}\), and absorbs constants. Setting \(c_1 := 4\) and \(c_2 := \min\{c_0, 1\}\) gives \[\mathbb{P}\left(\|\boldsymbol{\Delta}\|_2 \le C\left(\sigma_1(\mathbf{A})\,\rho\sqrt{d} + \rho^2 n\right)\right) \ge 1 - c_1 e^{-c_2 d}.\]


6.13 Proof of Lemma 7↩︎

For each \(z\in\Gamma\), set \[X(z) := \|\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\|_2,\] and let \(\mathcal{G}\) denote the high-probability event of Lemma 6, on which \(\|\boldsymbol{\Delta}\|_2\le C\left(\sigma_1\rho\sqrt d+\rho^2 n\right)\), with \(\mathbb{P}\left(\mathcal{G}^c\right)\le c_1 e^{-c_2 d}\).

6.13.0.1 Step 1: Pointwise contour bound on \(\left|I_\ell\right|\).

Fix \(\ell\ge 3\). Recall from ?? that \[I_\ell = \frac{1}{2\pi \mathrm{i}}\oint_\Gamma z\,\mathbb{E}\left[\mathop{\mathrm{tr}}\left((\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta})^\ell\mathbf{R}_\mathbf{G}(z)\right)\right]\,dz,\] so by the triangle inequality on the contour integral and Jensen’s inequality on the expectation, \[\label{eq:Iell-step1} \left|I_\ell\right| \le \frac{1}{2\pi}\oint_\Gamma \left|z\right|\,\mathbb{E}\left[\left|\mathop{\mathrm{tr}}\left((\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta})^\ell\mathbf{R}_\mathbf{G}(z)\right)\right|\right]\,\left|dz\right|.\tag{38}\] For any \(d\times d\) matrix \(\mathbf{M}\) one has \(\left|\mathop{\mathrm{tr}}(\mathbf{M})\right|\le d\|\mathbf{M}\|_2\), and submultiplicativity of the operator norm gives \[\|(\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta})^\ell\mathbf{R}_\mathbf{G}(z)\|_2 \le \|\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\|_2^\ell\,\|\mathbf{R}_\mathbf{G}(z)\|_2 = X(z)^\ell\,\|\mathbf{R}_\mathbf{G}(z)\|_2.\] Substituting into 38 , \[\label{eq:Iell-step1-final} \left|I_\ell\right| \le \frac{d}{2\pi}\oint_\Gamma \left|z\right|\,\mathbb{E}\left[X(z)^\ell\,\|\mathbf{R}_\mathbf{G}(z)\|_2\right]\,\left|dz\right| \le \frac{d}{\pi g}\oint_\Gamma \left|z\right|\,\mathbb{E}\left[X(z)^\ell\right]\,\left|dz\right|,\tag{39}\] where for the last inequality, we used the resolvent bound 12 \(\|\mathbf{R}_\mathbf{G}(z) \|_2\le \frac{2}{g}\).

6.13.0.2 Step 2: Bounding \(X(z)\) on \(\mathcal{G}\) and on \(\mathcal{G}^c\).

For every \(z\in\Gamma\), submultiplicativity and 12 gives \(X(z) \le \frac{2}{g}\|\boldsymbol{\Delta}\|_2\).

On \(\mathcal{G}\), the bound from Lemma 6 therefore yields \[X(z)\,\mathbf{1}_{\mathcal{G}} \le \frac{2}{g}\,C\left(\sigma_1\rho\sqrt d+\rho^2 n\right) = \alpha(\rho),\] with \(\alpha(\rho)\) as in 13 .

On \(\mathcal{G}^c\) the lemma no longer applies, but the entries of \(\mathbf{E}\) are bounded by \(\rho\), so the deterministic Frobenius bound \[\|\mathbf{E}\|_2\le \|\mathbf{E}\|_F \le \rho\sqrt{nd}\] still holds, and the triangle inequality on \(\boldsymbol{\Delta}=\mathbf{A}^\top\mathbf{E}+\mathbf{E}^\top\mathbf{A}+\mathbf{E}^\top\mathbf{E}\) gives \[\|\boldsymbol{\Delta}\|_2 \le 2\|\mathbf{A}\|_2\|\mathbf{E}\|_2+\|\mathbf{E}\|_2^2 \le 2\sigma_1\rho\sqrt{nd}+\rho^2 nd.\] Multiplying by \(2/g\), \[X(z)\,\mathbf{1}_{\mathcal{G}^c} \le \frac{2}{g}\left(2\sigma_1\rho\sqrt{nd}+\rho^2 nd\right) = \beta(\rho).\]

6.13.0.3 Step 3: Taking expectation and splitting on \(\mathcal{G}\).

Splitting the indicator \(1=\mathbf{1}_{\mathcal{G}}+\mathbf{1}_{\mathcal{G}^c}\) inside the expectation and using the bounds of Step 2, \[\begin{align} \mathbb{E}\left[X(z)^\ell\right] &= \mathbb{E}\left[X(z)^\ell\,\mathbf{1}_{\mathcal{G}}\right] + \mathbb{E}\left[X(z)^\ell\,\mathbf{1}_{\mathcal{G}^c}\right] \\ &\le \alpha(\rho)^\ell\,\mathbb{P}\left(\mathcal{G}\right) + \beta(\rho)^\ell\,\mathbb{P}\left(\mathcal{G}^c\right) \\ &\le \alpha(\rho)^\ell + c_1 e^{-c_2 d}\,\beta(\rho)^\ell, \end{align}\] where the last step uses \(\mathbb{P}\left(\mathcal{G}\right)\le 1\) and \(\mathbb{P}\left(\mathcal{G}^c\right)\le c_1 e^{-c_2 d}\) from Lemma 6. Substituting into 39 , \[\left|I_\ell\right| \le \frac{d}{\pi g}\left(\alpha(\rho)^\ell+c_1 e^{-c_2 d}\,\beta(\rho)^\ell\right)\oint_\Gamma\left|z\right|\,\left|dz\right|.\]

6.13.0.4 Step 4: Geometric parameters of the contour \(\Gamma\).

The contour \(\Gamma\) consists of two horizontal segments of length \(\lambda_1-\lambda_k\) joined by two half-circles of radius \(g/2\) (Figure 1), so \[\operatorname{len}(\Gamma) = 2(\lambda_1-\lambda_k)+\pi g.\] The point of \(\Gamma\) farthest from the origin is the rightmost point of the right half-circle, namely \(\lambda_1+g/2\), so \[\sup_{z\in\Gamma}\left|z\right| = \lambda_1+\frac{g}{2}.\] Bounding \(\oint_\Gamma\left|z\right|\left|dz\right|\le\operatorname{len}(\Gamma)\cdot\sup_{z\in\Gamma}\left|z\right|\) and substituting, \[\label{eq:Iell-step4} \left|I_\ell\right| \le \frac{d}{\pi g}\left(2(\lambda_1-\lambda_k)+\pi g\right)\left(\lambda_1+\frac{g}{2}\right)\left(\alpha(\rho)^\ell+c_1 e^{-c_2 d}\,\beta(\rho)^\ell\right).\tag{40}\]

6.13.0.5 Step 5: Geometric summation in \(\ell\).

Under the hypothesis \(\alpha(\rho),\beta(\rho)<1\), both geometric series converge and \[\sum_{\ell=3}^{\infty}\alpha(\rho)^\ell = \frac{\alpha(\rho)^3}{1-\alpha(\rho)}, \qquad \sum_{\ell=3}^{\infty}\beta(\rho)^\ell = \frac{\beta(\rho)^3}{1-\beta(\rho)}.\] Summing 40 over \(\ell\ge 3\), using \(\left|\sum_\ell I_\ell\right|\le\sum_\ell\left|I_\ell\right|\), and absorbing \(1/\pi\) into an absolute constant \(C\), \[\left|\sum_{\ell=3}^{\infty}I_\ell\right| \le \frac{Cd}{g}\left(2(\lambda_1-\lambda_k)+\pi g\right)\left(\lambda_1+\frac{g}{2}\right)\left(\frac{\alpha(\rho)^3}{1-\alpha(\rho)}+c_1 e^{-c_2 d}\,\frac{\beta(\rho)^3}{1-\beta(\rho)}\right).\]


6.14 Proof of Theorem 5↩︎

Let \(\mathcal{G}\) be the event of Lemma 6. Under the hypotheses \(\alpha(\rho)<1\) and \(\beta(\rho)<1\) with \(\alpha,\beta\) as in 13 , on \(\mathcal{G}\) \[\|\boldsymbol{\Delta}\|_2\le C\left(\sigma_1\rho\sqrt d+\rho^2 n\right) = \tfrac{g}{2}\,\alpha(\rho) < \tfrac{g}{2},\] while on \(\mathcal{G}^c\) the deterministic Frobenius bound \(\|\mathbf{E}\|_2\le\|\mathbf{E}\|_F\le\rho\sqrt{nd}\) together with \(\|\boldsymbol{\Delta}\|_2\le 2\|\mathbf{A}\|_2\|\mathbf{E}\|_2+\|\mathbf{E}\|_2^2\) gives \[\|\boldsymbol{\Delta}\|_2\le 2\sigma_1\rho\sqrt{nd}+\rho^2 nd = \tfrac{g}{2}\,\beta(\rho) < \tfrac{g}{2}.\] Hence \(\|\boldsymbol{\Delta}\|_2<g/2\) almost surely, so Lemma 2 applies and yields the convergent series ?? . Combining ?? with 11 , \[\mathcal{T}= \mathbb{E}\left[\|\mathbf{E}\|_F^2\right] - \sum_{\ell=1}^{\infty} I_\ell = \sum_{j=1}^{d}\nu_j - I_1 - I_2 - \sum_{\ell\ge 3} I_\ell,\] where \(\mathbb{E}\left[\|\mathbf{E}\|_F^2\right] = \mathop{\mathrm{tr}}(\mathbf{D}) = \sum_{j=1}^{d}\nu_j\). Bounding \(I_1\) from above by Lemma 3, \(I_2\) from above by Lemma 4, and \(\sum_{\ell\ge 3} I_\ell\) in absolute value by Lemma 7, and using \(\sum_{j=1}^{d}\nu_j - \sum_{j=1}^{k}\nu_j^\downarrow = \sum_{j=k+1}^{d}\nu_j^\downarrow\), \[\mathcal{T}\ge \sum_{j=k+1}^{d}\nu_j^\downarrow - \frac{C}{g}\left[\rho^2\sum_{i=1}^{k}\sum_{j=k+1}^{d}(\lambda_i+\lambda_j) + \sum_{j=1}^{k}(\nu_j^\downarrow)^2 + \rho^4 n\,k(d-k)\right] - \mathcal{R}(\rho).\]


6.15 Proof of Corollary 2↩︎

Let \(\bar c_1, \bar c_2 > 0\) denote the constants from Lemma 6.

Verifying the hypotheses of Theorem 5. Using \(\sigma_1 \le c_2\sigma_k\), \(\rho\sqrt{nd}\le c_1\sigma_k\), and \(n \ge d \ge 1\), \[\sigma_1\rho\sqrt d \le \frac{c_2\sigma_k\cdot c_1\sigma_k}{\sqrt n} \le c_1 c_2\,\sigma_k^2,\qquad \rho^2 n = \frac{(\rho\sqrt{nd})^2}{d} \le c_1^2\,\sigma_k^2,\] \[\sigma_1\rho\sqrt{nd} \le c_1 c_2\,\sigma_k^2,\qquad \rho^2 nd \le c_1^2\,\sigma_k^2.\] Substituting into 13 and using \(g \ge c_g\sigma_k^2\), \[\alpha(\rho) \le \frac{2C(c_1 c_2 + c_1^2)}{c_g},\qquad \beta(\rho) \le \frac{2(2 c_1 c_2 + c_1^2)}{c_g}.\] Choosing \(c_1\) sufficiently small, we may ensure \(\alpha(\rho), \beta(\rho) \le 1/2\) so that \(\tfrac{1}{1-\alpha(\rho)},\;\tfrac{1}{1-\beta(\rho)} \le 2\). In particular \(\alpha(\rho), \beta(\rho) < 1\), so the hypotheses of Theorem 5 are satisfied.

Bounding the leading terms. From \(\sigma_1 \le c_2\sigma_k\), \(\lambda_1 = \sigma_1^2 = O(\sigma_k^2)\) and hence \(\lambda_i+\lambda_j = O(\sigma_k^2)\) for all \(i\le k<j\). Therefore \[\frac{\rho^2}{g}\sum_{i=1}^{k}\sum_{j=k+1}^{d}(\lambda_i+\lambda_j) \le C\,k(d-k)\,\rho^2.\] Since each \(\nu_j \le n\rho^2/4\), \[\frac{1}{g}\sum_{j=1}^{k}(\nu_j^\downarrow)^2 \le C\frac{kn^2\rho^4}{\sigma_k^2}, \qquad \frac{\rho^4 n\,k(d-k)}{g} \le C\frac{kdn\,\rho^4}{\sigma_k^2}.\]

Bounding the remainder. Plugging the definitions 13 into Lemma 7 together with \(g=\Theta(\sigma_k^2)\) and \(\lambda_1=O(\sigma_k^2)\) gives \[\label{eq:R-cube-form} \mathcal{R}(\rho) \le \frac{Cd}{\sigma_k^4}\left[\left(\sigma_1\rho\sqrt d + \rho^2 n\right)^3 + \bar c_1 e^{-\bar c_2 d}\left(2\sigma_1\rho\sqrt{nd}+\rho^2 nd\right)^3\right].\tag{41}\] We claim the linear-in-\(\rho\) part dominates the quadratic-in-\(\rho\) part inside each cube. For the second cube, using \(\sigma_1\ge\sigma_k\) and \(\rho\sqrt{nd}\le c_1\sigma_k\), \[\frac{2\sigma_1\rho\sqrt{nd}}{\rho^2 nd} \;=\; \frac{2\sigma_1}{\rho\sqrt{nd}} \;\ge\; \frac{2\sigma_k}{c_1\sigma_k} \;=\; \frac{2}{c_1}.\] For the first cube, write \(\rho n = \rho\sqrt{nd}\cdot\sqrt{n/d}\le c_1\sigma_k\sqrt{n/d}\). Using \(\sigma_1\ge\sigma_k\) and the assumption \(d\ge c_3\sqrt n\), \[\frac{\sigma_1\rho\sqrt d}{\rho^2 n} \;=\; \frac{\sigma_1\sqrt d}{\rho n} \;\ge\; \frac{\sigma_k\sqrt d}{c_1\sigma_k\sqrt{n/d}} \;=\; \frac{1}{c_1}\cdot\frac{d}{\sqrt n} \;\ge\; \frac{c_3}{c_1}.\] Choosing \(c_1\le 1\) and \(c_3\ge 2c_1\), both ratios are \(\ge 2\), so \((a+b)^3\le 4 a^3\) for \(a\) the linear part and \(b\) the quadratic part. Thus \((\sigma_1\rho\sqrt d+\rho^2 n)^3\le 4(\sigma_1\rho\sqrt d)^3\) and \((2\sigma_1\rho\sqrt{nd}+\rho^2 nd)^3\le 4(2\sigma_1\rho\sqrt{nd})^3\). Substituting into 41 and using \(\sigma_1\le c_2\sigma_k\), \[\mathcal{R}(\rho) \le C\left[\frac{d^{5/2}\rho^3}{\sigma_k} + \bar c_1 e^{-\bar c_2 d}\,\frac{d^{5/2} n^{3/2}\rho^3}{\sigma_k}\right] \le C\,\frac{d^{5/2}\rho^3}{\sigma_k},\] where in the last step we absorb the exponentially small bad-event contribution. More precisely, by the assumption \(d\ge c_3 \sqrt{n}, \bar c_1 e^{-\bar c_2 d}\,n^{3/2}\) vanishes as \(n\rightarrow \infty\).

Substituting these bounds into Theorem 5 proves the claim.


6.16 Proof of Remark 6↩︎

Indeed, assume \(\sum_{j=k+1}^{d}\nu_j^\downarrow \ge (d-k) n \rho^2 \nu_0\) for some constant \(\nu_0\). Applying the assumption \(\rho \sqrt{n d}\le c_1 \sigma_k\) yields: \[\begin{align} \frac{k(d-k)\,\rho^2}{\sum_{j=k+1}^{d}\nu_j^\downarrow}\;&\le\; \frac{1}{\nu_0}\cdot\frac{k}{n}, \\ \frac{kn^2\rho^4/\sigma_k^2}{\sum_{j=k+1}^{d}\nu_j^\downarrow}\;&\le\; \frac{c_1^2}{\nu_0}\cdot\frac{k}{d\,(d-k)}, \\ \frac{d^{5/2}\rho^3/\sigma_k}{\sum_{j=k+1}^{d}\nu_j^\downarrow}\;&\le\; \frac{c_1}{\nu_0}\cdot\frac{d^{2}}{(d-k)\,n^{3/2}}. \end{align}\] Note that for the 3rd ratio to vanish, we need \(d = o (n^{3/4})\).

Recall the definition of \(\nu\) in 5 . For the special case \(k = d-1\), the leading term of our bound becomes: \[\begin{align} \sum_{j=k+1}^{d}\nu_j^\downarrow = \nu_d^\downarrow = \min_{1\le j\le d}\{\nu_j\} = n \rho^2 \nu = \left( \rho \sqrt{n \nu}\right)^2. \end{align}\]


6.17 Proof of Corollary 3↩︎

We specialize Corollary 2 to \(d=n\). The first assumption \(\rho\sqrt{nd}\le c_1\sigma_k\) reduces to \(\rho n\le c_1\sigma_k\), while the second is unchanged. Substituting \(d=n\) into the conclusion of Corollary 2, \[\mathcal{T}\ge \sum_{j=k+1}^{n}\nu_j^\downarrow - C\left(k(n-k)\,\rho^2 + \frac{kn^2\rho^4}{\sigma_k^2}+\frac{n^{5/2}\rho^3}{\sigma_k}\right).\] The middle term is absorbed into the last via \(\rho n\le c_1\sigma_k\): \[\frac{k n^2 \rho^4}{\sigma_k^2} = \frac{k n \rho^3}{\sigma_k^2} \cdot n \rho \le \frac{k n \rho^3}{\sigma_k^2} \cdot c_1 \sigma_k < c_1 \frac{n^2 \rho^3}{\sigma_k} < c_1 \frac{n^{5/2}\rho^3}{\sigma_k}.\] which gives the claimed bound.


6.18 Proof of Remark 7↩︎

Indeed, assume \(\sum_{j=k+1}^{d}\nu_j^\downarrow \ge (n-k) n \rho^2 \nu_0\) for some constant \(\nu_0\). Applying the assumption \(\rho n\le c_1 \sigma_k\) yields: \[\frac{k (n-k)\rho^2}{\sum_{j=k+1}^{n}\nu_j^\downarrow}\;\le\; \frac{1}{\nu_0}\cdot\frac{k}{n}, \qquad \frac{n^{5/2}\rho^3/\sigma_k}{\sum_{j=k+1}^{n}\nu_j^\downarrow}\;\le\; \frac{1}{\nu_0}\cdot\frac{n^{3/2}\rho}{(n-k)\,\sigma_k}\;\le\; \frac{c_1}{\nu_0}\cdot\frac{\sqrt n}{n-k}.\]


6.19 Alternative Bounds for Small Clusters↩︎

In this section we provide alternative bounds for the expected excess tail energy after rounding, \[\mathcal{T}:= \mathbb{E}\left[\sum_{i=k+1}^{d} \sigma_i(\tilde{\mathbf{A}})^2\right] - \sum_{i=k+1}^{d}\sigma_i(\mathbf{A})^2.\]

We start with an alternative upper bound for \(\|\boldsymbol{\Delta}\|_2\), and propagate the changes to Lemma 7, Theorem 5 and Corollary 2.

Lemma 10 (Alternative bound of \(\|\boldsymbol{\Delta}\|_2\)). Under the noise model of Section 2.2, with \(n \ge d\) and \(\mathbf{A}\) of full column rank, there exists a positive constant \(C\), such that: \[\mathbb{P}\left(\|\boldsymbol{\Delta}\|_2 \le C (\sigma_1(\mathbf{A}) \rho \sqrt{n} + \rho^2 n)\right) \ge 1 - 2 e^{-\left(n+d\right)}.\]

With triangle inequality and submultiplicity, we can bound \(\|\boldsymbol{\Delta}\|_2\) by \(2 \sigma_1(\mathbf{A}) \|\mathbf{E}\|_2+ \|\mathbf{E}\|_2^2\). Invoking Lemma 5 completes the proof.


Define another contraction parameter: \[\label{eq:alpha95tilde} \tilde{\alpha}(\rho):= \; \frac{2}{g}\,C\left(\sigma_1(\mathbf{A})\,\rho\sqrt{n} + \rho^2 n\right).\tag{42}\] On the high-probability event of Lemma 10, \(\|\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\|_2 \le \tilde{\alpha}(\rho)\) uniformly in \(z\in\Gamma\); on its complement, the deterministic Frobenius bound \(\|\mathbf{E}\|_2\le\|\mathbf{E}\|_F\le\rho\sqrt{nd}\) gives \(\|\mathbf{R}_\mathbf{G}(z)\boldsymbol{\Delta}\|_2\le\beta(\rho)\).

Lemma 11 (Higher-order remainder alternative). Assume \(\tilde{\alpha}(\rho) < 1\) and \(\beta(\rho) < 1\) with \(\tilde{\alpha},\beta\) as in 1342 . Then there is an absolute constant \(C > 0\) such that \[\left|\sum_{\ell=3}^{\infty} I_\ell\right| \le \frac{Cd}{g}\left(2(\lambda_1 - \lambda_k) + \pi g\right)\left(\lambda_1 + \frac{g}{2}\right)\left(\frac{\tilde{\alpha}(\rho)^3}{1 - \tilde{\alpha}(\rho)} + 2 e^{-(n+d)}\,\frac{\beta(\rho)^3}{1 - \beta(\rho)}\right).\]

Follow the same steps as in Appendix 6.13.


Theorem 9. Assume \(\tilde{\alpha}(\rho) < 1\) and \(\beta(\rho) < 1\) with \(\tilde{\alpha}, \beta\) as in 1342 . Then \[\mathcal{T}\ge \sum_{j=k+1}^{d}\nu_j^\downarrow - \frac{C}{g}\left[\rho^2 \sum_{i=1}^{k}\sum_{j=k+1}^{d}(\lambda_i + \lambda_j) + \sum_{j=1}^{k}\left(\nu_j^\downarrow\right)^2 + \rho^4 n\, k(d-k)\right] - \tilde{\mathcal{R}}(\rho),\] where the higher-order remainder \(\tilde{\mathcal{R}}(\rho)\) is the right-hand side of Lemma 11.

Follow the same steps as in Appendix 6.14.


Corollary 4. There exist absolute constants \(c_1, c_2, c_g > 0\) such that, if \[\rho\sqrt{nd} \le c_1\,\sigma_k, \qquad \sigma_1 \le c_2\,\sigma_k, \qquad g > c_g \lambda_k,\] then there is a constant \(C > 0\) such that \[\mathcal{T}\ge \sum_{j=k+1}^{d}\nu_j^\downarrow - C\left(k(d-k)\,\rho^2 + \frac{k\,n^2 \rho^4}{\sigma_k^2} + \frac{d n^{3/2}\rho^3}{\sigma_k}\right).\]

Verifying the hypotheses of Theorem 9. Using \(\sigma_1 \le c_2\sigma_k\), \(\rho\sqrt{nd} \le c_1\sigma_k\), and \(n \ge d \ge 1\), \[\sigma_1\rho\sqrt n \le \frac{c_2\sigma_k\cdot c_1\sigma_k}{\sqrt d} \le c_1 c_2\,\sigma_k^2,\qquad \rho^2 n = \frac{(\rho\sqrt{nd})^2}{d} \le c_1^2\,\sigma_k^2,\] \[\sigma_1\rho\sqrt{nd} \le c_1 c_2\,\sigma_k^2,\qquad \rho^2 nd \le c_1^2\,\sigma_k^2.\] Substituting into 42 and 13 , and using \(g \ge c_g\sigma_k^2\), \[\tilde{\alpha}(\rho) \le \frac{2C(c_1 c_2 + c_1^2)}{c_g},\qquad \beta(\rho) \le \frac{2(2 c_1 c_2 + c_1^2)}{c_g}.\] Choosing \(c_1\) sufficiently small, we may ensure \(\tilde{\alpha}(\rho), \beta(\rho) \le 1/2\) so that \(\tfrac{1}{1-\tilde{\alpha}(\rho)},\;\tfrac{1}{1-\beta(\rho)} \le 2\). The hypotheses of Theorem 9 are thus satisfied.

Bounding the leading terms. From \(\sigma_1 \le c_2\sigma_k\), \(\lambda_1 = \sigma_1^2 = O(\sigma_k^2)\) and hence \(\lambda_i + \lambda_j = O(\sigma_k^2)\) for all \(i \le k < j\). Therefore \[\frac{\rho^2}{g}\sum_{i=1}^{k}\sum_{j=k+1}^{d}(\lambda_i+\lambda_j) \le C\,k(d-k)\,\rho^2.\] Since each \(\nu_j \le n\rho^2/4\), \[\frac{1}{g}\sum_{j=1}^{k}(\nu_j^\downarrow)^2 \le C\frac{kn^2\rho^4}{\sigma_k^2}, \qquad \frac{\rho^4 n\,k(d-k)}{g} \le C\frac{\rho^4\,n\,k(d-k)}{\sigma_k^2} \le C\frac{kn^2\rho^4}{\sigma_k^2}.\]

Bounding the remainder. Plugging 42 and 13 into Lemma 11 together with \(g = \Theta(\sigma_k^2)\) and \(\lambda_1 = O(\sigma_k^2)\) gives \[\label{eq:R-cube-form-alter} \tilde{\mathcal{R}}(\rho) \le \frac{Cd}{\sigma_k^4}\left[\left(\sigma_1\rho\sqrt n + \rho^2 n\right)^3 + 2 e^{-(n+d)}\left(2\sigma_1\rho\sqrt{nd}+\rho^2 nd\right)^3\right].\tag{43}\] We claim the linear-in-\(\rho\) part dominates the quadratic-in-\(\rho\) part inside each cube. For the second cube, using \(\sigma_1 \ge \sigma_k\) and \(\rho\sqrt{nd} \le c_1\sigma_k\), \[\frac{2\sigma_1\rho\sqrt{nd}}{\rho^2 nd} \;=\; \frac{2\sigma_1}{\rho\sqrt{nd}} \;\ge\; \frac{2\sigma_k}{c_1\sigma_k} \;=\; \frac{2}{c_1}.\] For the first cube, write \(\rho\sqrt n = \rho\sqrt{nd}/\sqrt d \le c_1\sigma_k/\sqrt d\). Using \(\sigma_1 \ge \sigma_k\) and \(d \ge 1\), \[\frac{\sigma_1\rho\sqrt n}{\rho^2 n} \;=\; \frac{\sigma_1}{\rho\sqrt n} \;\ge\; \frac{\sigma_k}{c_1\sigma_k/\sqrt d} \;=\; \frac{\sqrt d}{c_1}.\] (In contrast to Corollary 2, no separate lower bound on \(d/\sqrt n\) is required: the alternative remainder uses \(\sqrt n\) rather than \(\sqrt d\) inside \(\tilde{\alpha}\), so dominance follows directly from \(\rho\sqrt{nd}\le c_1\sigma_k\).) Choosing \(c_1 \le 1/2\), both ratios are \(\ge 2\), so \((a+b)^3 \le 4 a^3\) for \(a\) the linear part and \(b\) the quadratic part. Thus \((\sigma_1\rho\sqrt n + \rho^2 n)^3 \le 4(\sigma_1\rho\sqrt n)^3\) and \((2\sigma_1\rho\sqrt{nd} + \rho^2 nd)^3 \le 4(2\sigma_1\rho\sqrt{nd})^3\). Substituting into 43 and using \(\sigma_1 \le c_2\sigma_k\), \[\tilde{\mathcal{R}}(\rho) \le C\left[\frac{d\, n^{3/2}\rho^3}{\sigma_k} + 2 e^{-(n+d)}\,\frac{d^{5/2} n^{3/2}\rho^3}{\sigma_k}\right] \le C\,\frac{d\, n^{3/2}\rho^3}{\sigma_k},\] where the last step absorbs the exponentially small bad-event contribution: since \(d \le n\), \(e^{-(n+d)} d^{3/2} \le e^{-n} n^{3/2} \to 0\) as \(n \to \infty\).

Substituting these bounds into Theorem 9 proves the claim.


References↩︎

[1]
P. Micikevicius et al., “Mixed precision training,” in International conference on learning representations, 2018.
[2]
T. Dettmers, M. Lewis, Y. Belkada, and L. Zettlemoyer, GPT3.int8(): 8-bit matrix multiplication for transformers at scale,” Advances in neural information processing systems, vol. 35, pp. 30318–30332, 2022.
[3]
A. Gholami, S. Kim, Z. Dong, Z. Yao, M. W. Mahoney, and K. Keutzer, “A survey of quantization methods for efficient neural network inference,” in Low-power computer vision, Chapman; Hall/CRC, 2022, pp. 291–326.
[4]
R. Jin et al., “A comprehensive evaluation of quantization strategies for large language models,” in Findings of the association for computational linguistics: ACL 2024, 2024, pp. 12186–12215.
[5]
A. Gupta et al., Effective Quantization of Muon Optimizer States,” arXiv preprint arXiv:2509.23106, 2025.
[6]
S. Gupta, A. Agrawal, K. Gopalakrishnan, and P. Narayanan, “Deep learning with limited numerical precision,” in International conference on machine learning, 2015, pp. 1737–1746.
[7]
L. Xia, S. Massei, M. E. Hochstenbach, and B. Koren, “On stochastic roundoff errors in gradient descent with low-precision computation,” Journal of Optimization Theory and Applications, vol. 200, no. 2, pp. 634–668, 2024.
[8]
K. Ozkara, T. Yu, and Y. Park, Stochastic Rounding for LLM Training: Theory and Practice,” in International conference on artificial intelligence and statistics, 2025, pp. 4402–4410.
[9]
G. E. Forsythe, “Reprint of a note on rounding-off errors,” SIAM Review, vol. 1, no. 1, pp. 66–67, 1959.
[10]
J. von Neumann and H. H. Goldstine, “Numerical inverting of matrices of high order,” Bulletin of the American Mathematical Society, 1947.
[11]
I. C. Ipsen and H. Zhou, “Probabilistic error analysis for inner products,” SIAM journal on matrix analysis and applications, vol. 41, no. 4, pp. 1726–1741, 2020.
[12]
M. P. Connolly, N. J. Higham, and T. Mary, “Stochastic rounding and its probabilistic backward error analysis,” SIAM Journal on Scientific Computing, vol. 43, no. 1, pp. A566–A585, 2021.
[13]
G. Dexter, C. Boutsikas, L. Ma, I. C. Ipsen, and P. Drineas, “Stochastic rounding implicitly regularizes tall-and-thin matrices,” SIAM Journal on Matrix Analysis and Applications, vol. 46, no. 1, pp. 341–369, 2025.
[14]
M. Mahoney and C. Martin, “Traditional and heavy tailed self regularization in neural network models,” in International conference on machine learning, 2019, pp. 4284–4293.
[15]
C. H. Martin and M. W. Mahoney, “Implicit self-regularization in deep neural networks: Evidence from random matrix theory and implications for learning,” Journal of Machine Learning Research, vol. 22, no. 165, pp. 1–73, 2021.
[16]
C. H. Martin, T. Peng, and M. W. Mahoney, “Predicting trends in the quality of state-of-the-art neural networks without access to training or testing data,” Nature Communications, vol. 12, no. 1, p. 4122, 2021.
[17]
D. He et al., Alphadecay: Module-wise weight decay for heavy-tailed balancing in LLMs,” International Journal of Artificial Intelligence and Robotics Research, 2026.
[18]
R. Saha, V. Srivastava, and M. Pilanci, “Matrix compression via randomized low rank and low precision factorization,” Advances in Neural Information Processing Systems, vol. 36, 2023.
[19]
M. Kwun, D. Morwani, H. Su, S. Gil, N. Anand, and S. M. Kakade, LOTION: Smoothing the Optimization Landscape for Quantized Training,” in OPT 2025: Optimization for machine learning.
[20]
K. Zhao, T. Tabaru, K. Kobayashi, T. Honda, M. Yamazaki, and Y. Tsuruoka, Direct Quantized Training of Language Models with Stochastic Rounding,” in The 17th asian conference on machine learning (conference track), 2025.
[21]
T. Brailovskaya and R. van Handel, “Universality and sharp matrix concentration inequalities,” Geometric and Functional Analysis, vol. 34, no. 6, pp. 1734–1838, 2024.
[22]
P. Tran and V. Vu, Davis–Kahan theorem under a moderate gap condition,” Communications in Contemporary Mathematics, vol. 28, no. 1, p. 2550035, 2026.
[23]
P. Tran and V. Vu, “Eigenvalue stability and new perturbation bounds for the extremal eigenvalues of a matrix.” 2026, [Online]. Available: https://arxiv.org/abs/2603.19758.
[24]
P. Tran and V. Vu, “New matrix perturbation bounds with relative norm: Perturbation of eigenspaces.” 2026, [Online]. Available: https://arxiv.org/abs/2409.20207.
[25]
P. Tran and V. Vu, “Matrix perturbation bounds via contour bootstrapping.” 2024, [Online]. Available: https://arxiv.org/abs/2407.05230.
[26]
P. Tran, V. Vu, and N. K. Vishnoi, “Spectral perturbation bounds for low-rank approximation with applications to privacy,” in The thirty-ninth annual conference on neural information processing systems, 2025.
[27]
P. Tran and V. Vu, “New perturbation bounds for low rank approximation of matrices via contour analysis,” arXiv preprint arXiv:2511.08875, 2025.
[28]
R. Vershynin, High-dimensional probability: An introduction with applications in data science, 2nd ed. Cambridge: Cambridge University Press, 2026.
[29]
C. Boutsikas, P. Drineas, and I. C. Ipsen, “Small singular values can increase in lower precision,” SIAM Journal on Matrix Analysis and Applications, vol. 45, no. 3, pp. 1518–1540, 2024.
[30]
S. B. Ali, S.-I. Filip, and O. Sentieys, “A stochastic rounding-enabled low-precision floating-point mac for dnn training,” in 2024 design, automation & test in europe conference & exhibition (DATE), 2024, pp. 1–6.
[31]
J. Green, The November meeting in Pasadena.” 1950.
[32]
E. Hallman and I. C. Ipsen, “Precision-aware deterministic and probabilistic error bounds for floating point summation,” Numerische Mathematik, vol. 155, no. 1, pp. 83–119, 2023.
[33]
P. de Oliveira Castro, E.-M. El Arar, E. Petit, and D. Sohier, “Error analysis of sum-product algorithms under stochastic rounding,” SIAM Journal on Scientific Computing, vol. 47, no. 6, pp. B1481–B1502, 2025.
[34]
E.-M. El Arar, D. Sohier, P. de Oliveira Castro, and E. Petit, “Stochastic rounding variance and probabilistic bounds: A new approach,” SIAM Journal on Scientific Computing, vol. 45, no. 5, pp. C255–C275, 2023.
[35]
M. Croci and M. B. Giles, “Effects of round-to-nearest and stochastic rounding in the numerical solution of the heat equation in low precision,” IMA Journal of Numerical Analysis, vol. 43, no. 3, pp. 1358–1390, 2023.
[36]
M. Croci, M. Fasi, N. J. Higham, T. Mary, and M. Mikaitis, “Stochastic rounding: Implementation, error analysis and applications,” Royal Society Open Science, vol. 9, no. 3, 2022.
[37]
P. Drineas and I. C. Ipsen, “Stochastic rounding 2.0, with a view towards complexity analysis,” SIAM News, vol. 57, no. 9, 2024.
[38]
D. Kalamkar et al., A study of BFLOAT16 for deep learning training,” arXiv preprint arXiv:1905.12322, 2019.
[39]
P. Micikevicius et al., FP8 formats for deep learning,” arXiv preprint arXiv:2209.05433, 2022.
[40]
A. Mishra, D. Stosic, S. Layton, and P. Micikevicius, Recipes for Pre-training LLMs with MXFP8,” arXiv preprint arXiv:2506.08027, 2025.
[41]
A. Tseng, T. Yu, and Y. Park, “Training LLMs with MXFP4,” arXiv preprint arXiv:2502.20586, 2025.
[42]
L. Xia, S. Massei, and M. E. Hochstenbach, On the convergence of the gradient descent method with stochastic fixed-point rounding errors under the Polyak–Łojasiewicz inequality,” Computational Optimization and Applications, vol. 90, no. 3, pp. 753–799, 2025.
[43]
S. Boucheron, G. Lugosi, and P. Massart, Concentration inequalities: A nonasymptotic theory of independence. Oxford University Press, 2013.
[44]
K. Fan, On a theorem of Weyl concerning eigenvalues of linear transformations I,” Proceedings of the National Academy of Sciences, vol. 35, no. 11, pp. 652–655, 1949.
[45]
K. B. Petersen and M. S. Pedersen, The Matrix Cookbook,” Technical University of Denmark, vol. 7, no. 15, p. 510, 2008.

  1. Equal contribution.↩︎

  2. Up to lower order terms. Here \(\rho\) is the resolution of the quantizer and \(\nu\) is a parameter measuring the column-wise variance of \(\mathbf{A}- \tilde{\mathbf{A}}\), which effectively measures the amount of available stochasticity when \(\mathbf{A}\) is rounded to form \(\tilde{\mathbf{A}}\). See Section 2.2 for the exact definition of \(\nu\).↩︎

  3. Up to lower order terms.↩︎

  4. Up to logarithmic terms.↩︎

  5. Here we refer to the first version of this paper on arXiv.↩︎

  6. The article can also be found here.↩︎

  7. Up to lower-order terms.↩︎

  8. Indeed, by Cauchy–Schwarz, \((a+b+c)^2=\bigl((a,b,c)\cdot(1,1,1)\bigr)^2\le (a^2+b^2+c^2)(1^2+1^2+1^2)=3(a^2+b^2+c^2).\)↩︎