arXiv CSopen access

Uniform Sampling from High-dimensional Spectral Norm Balls

Uniform Sampling from High-dimensional Spectral Norm Balls arXiv:2606.24134v1 [math.PR] 23 Jun 2026

Michael R. Metel∗ June 24, 2026

Abstract Motivated by an application in machine learning optimization, this paper focuses on the challenges of sampling a matrix uniformly from the unit spectral norm ball. It is proven that all singular values of sampled matrices converge to 1 almost surely as the matrix dimensions increase. This result provides the theoretical justification for a proposed simple sampling method applicable for large dimension sizes matching matrices found in modern large language models. Experimental results demonstrate both the convergence of the singular values, as well as the exact and proposed approximate sampling methods.

1

Introduction

This paper studies the uniform sampling of matrices ∆ ∈ Rn×m from the unit spectral norm ball, B := {∆ ∈ Rn×m : σmax (∆) ≤ 1}, where σmax (∆) denotes the largest singular value of ∆. A sample from this distribution will ˆ ∼ U(B). Given that a sample from the spectral norm ball of radius ρ > 0 be denoted as ∆ ˆ the focus on ρ = 1 is without loss of generality. After a possible can be generated as ρ∆, transpose, it will also be assumed that m ≥ n. A polynomial-time recipe to sample from this distribution can be found across two papers (Calafiore and Dabbene, 2002; Calafiore et al., 2000), which was motivated by probabilistic robustness techniques for control systems. The idea is to consider a normalized singular value decomposition of ∆, ∆ = U ΣV T , where U ∈ Rn,n and V ∈ Rm,n have orthonormal columns with V1,j > 0 (almost surely) for j = 1, ..., n, and Σ = diag(σ1 , .., σn ) for singular values σ1 ≥ σ2 ≥ ... ≥ σn . It is proven (Calafiore et al., 2000, Theorem 2) that to uniformly sample over B, U , Σ, and V ∗

[email protected]

1

can be sampled independently, with U and V sampled uniformly over the orthogonal and normalized semi-orthogonal matrices, respectively, which can be done using samples from the standard Normal distribution and QR decomposition as shown in Algorithm 1. Algorithm 1 Uniformly sampling U and V Sample Û ∈ Rn,n and V̂ ∈ Rm,m : Ûj,k ∼ N (0, 1) and V̂j,k ∼ N (0, 1) ∀j, k [Q(Û ), R(Û )] = QR(Û ) [Q(V̂ ), R(V̂ )] = QR(V̂ ) Û = Q(Û ) diag(sign(R(Û )11 , ..., sign(R(Û )nn ) V̂ = Q(V̂ )[:, 1 : n] V̂ = V̂ diag(sign(V̂11 ), ..., sign(V̂1n )) Output: Û ,V̂ In order to sample Σ requires significantly more effort. Singular values are proposed to be sampled iteratively following the conditional distribution method (Devroye, 1986, Chapter 11). The algorithm presented in (Calafiore and Dabbene, 2002, page 1269) samples singular values using exact conditional densities which contain both constants and exponents which result in over and underflows for modest problem sizes. As the motivation is to sample matrices at the scale of modern large language models, e.g., the largest linear layer weight matrix of the 32 billion-parameter Qwen 3 model (Qwen Team, 2025) is n × m = 5, 120 × 25, 600, the first step, detailed in Section 2, was to simplify the algorithm by using functions which are only proportional to the exact conditional densities. After failing to scale up the algorithm directly, the distribution of the singular values as the dimensions of ∆ increase is studied in Section 3. The main contribution presented as Theorem 1 is that all singular values converge to 1 almost surely as n → ∞, and in Theorem 2, it is proven that all singular values converge to 1 almost surely as only m → ∞ for a fixed value of n. These results present a viable approximate sampling algorithm consisting of setting Σ = I, and simply using Algorithm 1 to sample from U(B). Section 4 contains further implementation details of the exact sampling algorithm, and empirical results demonstrating the convergence of the exact and proposed approximate sampling methods. There appears to have been no attempts to develop a practical approach to sample from U(B) at the scale of real-world problems. The software package RACT (Tremba et al., 2008) supports sampling from U(B), but a warning stating “This can take a LONG time...” will be given when min(m, n) > 6. The experiments presented in (Calafiore and Dabbene, 2002) are also restricted to problem sizes upper bounded by max(n, m) ≤ 4. This section ends with an application that motivated this work, requiring an efficient sampling approach for high-dimensional matrices.

1.1

Neural Network Training with Random Weight Perturbation

Gaussian, uniform l2 , and l∞ -norm ball perturbations have been applied to loss function decision variables when computing stochastic gradients for first-order optimization methods (Duchi et al., 2012, Appendix E). This has a smoothing effect, resulting in non-differentiable 2

Lipchitz continuous functions becoming smooth in expectation, enabling provably convergent optimization algorithms. From a neural network perspective, random weight perturbation pushes convergence to flatter local minima, improving a trained model’s ability to generalize (Bisla et al., 2022). These forms of perturbation implicitly view a model’s weights as a vector w ∈ RN . Considering that neural networks contain several layers, each of which consisting of matrices, a line of research (Bernstein and Newhouse, 2025) has focused on incorporating this structure into the design of optimization methods. Letting each weight matrix W i live in its own vector space V i with norm ∥ · ∥i , the neural network is contained in the product space Πi V i , with norm max ai ∥W i ∥i for per-matrix weights ai > 0. Assuming i

that the input xi ∈ Rm and output y i ∈ Rn of each linear layer lie in Euclidean space, ∥ · ∥i is set to the induced matrix norm of W i : the spectral norm. Given an efficient method to ˆ i ∼ U(B), W i can then be perturbed by ρi ∆ ˆ i for per-matrix radii ρi , enabling the sample ∆ extension of random weight perturbation to architecture-aware optimization algorithms.

2

Singular Value Sampling Algorithm

This section describes an implementation of the algorithm described in (Calafiore and Dabbene, 2002, page 1269) to sample iteratively the singular values σ1 ≥ σ2 ≥ ... ≥ σn . While the differences between the described implementation and the published method are described, for a full derivation of this algorithm, it is recommended to consult the original publication.

2.1

Definitions 

 0 1 := , and 1n := [1, 1, ..., 1]T ∈ Rn . Let Let v := −1 0 i−j , let F ∈ Rn be a vector S ∈ Rn×n be a skew-symmetric matrix, where Sij := (i+v)(j+v)(i+j+2v) 1 defined as Fj := j+v , and let X : R → Rn be defined as Xj (x) := xj−1 . When sampling σi , the algorithm depends on the parity of n and i, so the resulting four cases will need to be considered. The functions defined below are used for the recursive definitions of the marginal densities of the singular values. 1 (m − n − 1), β := n2 (n − 1), J 2

Even number of rows n: Let β β−1 Z0 (x) := diag(xβ , xβ−1 , ..., xβ−n−1 )Q−1 , ..., xβ−n−1 ), 0 diag(x , x

or element-wise, (Z0 (x))ij = x2β−i−j+2 (Q−1 0 )ij , where Q0 := S, and let he (x) := X(x). The functions Z0 and he are used to generate even-indexed samples σ̂i . Let β β−1 Z1 (x) := diag(xβ , xβ−1 , ..., xβ−n−1 , xβ , xβ )Q−1 , ..., xβ−n−1 , xβ , xβ ), 1 diag(x , x

or  2β−i−j+2 −1 x (Q1 )ij    x2β−i+1 (Q−1 ) 1 ij (Z1 (x))ij := 2β−j+1 x (Q−1  1 )ij   2β −1 x (Q1 )ij

if i, j ∈ [1, 2, ..., n] if i ∈ [1, 2, ..., n] and j ∈ [n + 1, n + 2] if j ∈ [1, 2, ..., n] and i ∈ [n + 1, n + 2] if i, j ∈ [n + 1, n + 2], 3

(1)

where



 S F 1n Q1 := −F T 0 0  , −1Tn 0 0

and let ho (x) := [X(x); 0; 0], which are used to generate odd-indexed samples σ̂i . Odd number of rows n: Let β β−1 , ..., xβ−n−1 , xβ ), Z0 (x) := diag(xβ , xβ−1 , ..., xβ−n−1 , xβ )Q−1 0 diag(x , x

where



 S F Q0 := , −F T 0 β β−1 , ..., xβ−n−1 , xβ ), Z1 (x) := diag(xβ , xβ−1 , ..., xβ−n−1 , xβ )Q−1 1 diag(x , x

where



 S 1n Q1 := , −1Tn 0

and he (x) := ho (x) := [X(x); 0]. The exponent of x in (Z0 (x))ij and (Z1 (x))ij for all 1 ≤ i, j ≤ n + 1 is identical to (1). In the original algorithm, Z0 and Z1 contain the matrices Adj(Q0 ) and Adj(Q1 ) instead −1 of Q−1 0 and Q1 , with our implementation not including the constants det(Q0 ) and det(Q1 ), which both underflow when n ≥ 28 using 64-bit floating point arithmetic1 . Let the scalar function p(x) := xβ p β det(Q0 ) and p1 (x) = also be defined, which will replace the scalar functions p (x) = x 0 p β x det(Q1 ) used in the original algorithm. The next section verifies that these simplifications preserve the proportionality to the conditional probability distribution of all σi .

2.2

Algorithm Implementation

With the given definitions, the basic algorithm is shown in Algorithm 2. 2(v+1)(n−1)

K σ1 p1 (σ12 )σ1m−n , Line 1: The probability density function of σ1 is equal to f1 (σ1 ) = 2n−1 where K is a normalization constant (Calafiore and Dabbene, 2002, Eq. 5). Removing constants and simplifying the exponents results in the given function proportional to f1 .

Lines 3, 5, & 11: These are exponents that σi will be taken to in its conditional density. Simplifying, e1 (i) = (m − n)(n − i + 1) + n − i, and by considering the case m = n, e1 (i) ≥ n − i. When i is even, e2 (i) = −(i−1)n(n−1), and when odd, e2 (i) = −(i−2)n(n−1). The exponent e1 (i) + e2 (i) is decreasing in i: For even i, e1 (i + 1) + e2 (i + 1) − (e1 (i) + e2 (i)) = −(m − n) − 1, and for odd i, e1 (i + 1) + e2 (i + 1) − (e1 (i) + e2 (i)) = −(m − n) − 1 − 2n(n − 1). 1

This was tested using numpy.linalg.det in the same Python environment described in Section 4.

4

Line 16: This is the recursive calculation of Zi (x) as given in (Calafiore and Dabbene, 2002, page 1269). Defining Z̃0 (x) := det(Q0 )Z0 (x) and Z̃1 (x) := det(Q1 )Z1 (x) as the original defi

initions of Z0 and Z1 , and using the same recursion to compute Z̃i , Z̃i (x) = det(Q0 )2 2 Zi (x) i−1

when i is even, and Z̃i (x) = det(Q1 )2 2 Zi (x) when i is odd. Lines 7 & 13: The marginal density of {σj }ij=1 is equal to K 2(v+1)(n−i)+m−n i−1 m−n Πj=1 σj , pi (σi2 )σi n−i 2 2 )T Z̃i−2 (x)he (x)| when i is even and where pi (x) = p0 (x)−(i−1) |he (σi−1 2 2 )| when i is odd. Given that {σj = σ̂j }i−1 )T Z̃i−2 (x)ho (σi−1 pi (x) = p1 (x)−(i−2) |ho (σi−2 j=1 are known, the required conditional density of σi is proportional to this marginal density. By m−n K removing the constant 2n−i Πi−1 , and recalling that p0 (x) & p1 (x) ∝ xβ and Zi−2 ∝ j=1 σj Z̃i−2 , the functions on Lines 7 & 13 are proportional to fi (σi |σ1 = σ̂1 , ..., σi−1 = σ̂i−1 ).

Algorithm 2 Sampling singular values σ1 ≥ σ2 ≥ ... ≥ σn 1: Generate sample σ̂1 from density f1 (σ1 ) ∝ σ1mn−1 2: for i = 2, ..., n do 3: e1 (i) = 2(v + 1)(n − i) + m − n 4: if i mod 2 = 0 then 5: e2 (i) = −2(i − 1)β 2 6: di (x) = he (σ̂i−1 )T Zi−2 (x)he (x) e (i)+e2 (i) 7: Generate sample σ̂i from density fi (σi |σ1 = σ̂1 , ..., σi−1 = σ̂i−1 ) ∝ σi 1 |di (σi2 )| 2 T 2 8: Update di (x) = he (σ̂i−1 ) Zi−2 (x)he (σ̂i ) 2 2 ), he (σ̂i2 )]T ), he (σ̂i2 )]J [he (σ̂i−1 9: H(i) = [he (σ̂i−1 10: else 11: e2 (i) = −2(i − 2)β 2 2 12: di (x) = ho (σ̂i−2 )T Zi−2 (x)ho (σ̂i−1 ) e (i)+e2 (i) 13: Generate sample σ̂i from density fi (σi |σ1 = σ̂1 , ..., σi−1 = σ̂i−1 ) ∝ σi 1 |di (σi2 )| 2 2 2 2 14: H(i) = [ho (σ̂i−2 ), ho (σ̂i−1 )]J [ho (σ̂i−2 ), ho (σ̂i−1 )]T 15: end if 16: Zi (x) = di (x)Zi−2 (x) + Zi−2 (x)H(i)Zi−2 (x) 17: end for

3

Convergence of {σi} in High Dimension

The main result of this work is that as n increases, all singular values converge to 1 almost surely, where we recall that σ1 ≥ σ2 ≥ ... ≥ σn by definition. Theorem 1. For matrices uniformly distributed in the unit spectral norm ball, ∆ ∼ U(B), lim σn = 1 almost surely. n→∞

5

Proof. From the definitions of Z0 and Z1 , the recursions defining Zi for i > 1, and the definitions of di , all elements of Zi and the functions di are polynomials in x, with all σi sampled from univariate polynomial functions. Let the unnormalized conditional density e (i)+e2 (i) function of σi , f˜i (σi |σ1 , ..., σi−1 ) := σi 1 |di (σi2 )|, with domain [0, σi−1 ], be written as X t aij σi ij , j

where all coefficients aij > 0. The polynomial changes with n, but as will be seen, the result is independent of the number of polynomial terms or the values of their coefficients. By normalization, fi (σi |σ1 , ..., σi−1 ) is equal to P tij j aij σi P aij tij +1 . j tij +1 σi−1 Viewing fi (σi |σ1 , ..., σi−1 ) as a weighted average of the terms t

aij σi ij

t +1 = aij σ ij tij +1 i−1

tij + 1 σi tij ( ) σi−1 σi−1

with weights t +1 aij σ ij tij +1 i−1 P aik tik +1 , k tik +1 σi−1

fi (σi |σ1 , ..., σi−1 ) can be bounded as inf j

tij + 1 σi tij tij + 1 σi tij ( ) ≤ fi (σi |σ1 , ..., σi−1 ) ≤ sup ( ) . σi−1 σi−1 σi−1 σi−1 j

(2)

The task now is to lower bound the exponents {tij }j for each i. By showing that these exponents are growing sufficiently fast with n, the theorem can be proven. To first compute the entries of Zi (x) = di (x)Zi−2 (x) + Zi−2 (x)H(i)Zi−2 (x), the values of di are derived: When i and n are even, 2(n−1)

2 2 di (x) = X(σi−1 )T Zi−2 (x)X(x) = [1, σi−1 , ..., σi−1 n n X X 2(p−1) = σi−1 xq−1 (Zi−2 (x))p,q ,

]Zi−2 (x)[1, x, ..., xn−1 ]T

p=1 q=1

and when n is odd, 2(n−1)

2 2 di (x) = [X(σi−1 )T , 0]Zi−2 (x)[X(x); 0] = [1, σi−1 , ..., σi−1 n X n X 2(p−1) = σi−1 xq−1 (Zi−2 (x))p,q . p=1 q=1

6

, 0]Zi−2 (x)[1, x, ..., xn−1 , 0]T

When i is odd and n is even, 2 2 ); 0; 0] = )T , 0, 0]Zi−2 (x)[X(σi−1 di (x) = [X(σi−2

n X n X

2(p−1) 2(q−1) σi−1 (Zi−2 (x))p,q ,

σi−2

p=1 q=1

and when n is odd, 2 2 ); 0] = )T , 0]Zi−2 (x)[X(σi−1 di (x) =[X(σi−2

n X n X

2(p−1) 2(q−1) σi−1 (Zi−2 (x))p,q .

σi−2

p=1 q=1

For the entries of Zi , when i is even, (Zi (x))jk = (

n X n X

2(p−1) q−1

σi−1

x

(Zi−2 (x))p,q )(Zi−2 (x))jk

(3)

p=1 q=1

+

n X n X

(Zi−2 (x))jp H(i)pq (Zi−2 (x))qk ,

p=1 q=1

and when i is odd, (Zi (x))jk = ( +

n X n X

2(p−1) 2(q−1) σi−1 (Zi−2 (x))p,q )(Zi−2 (x))jk

σi−2

(4)

p=1 q=1 n X n X

(Zi−2 (x))jp H(i)pq (Zi−2 (x))qk ,

p=1 q=1

where the summations involving H(i) only go up to n given that H(i)pq = 0 when p or −1 q > n. The matrices Q−1 0 and Q1 in the definitions of Z0 (x) and Z1 (X), and the matrices H(i) are skew-symmetric, hence zero-diagonal. For simplicity, this will be ignored by lower bounding all exponents {tij }j of each σi , even those whose corresponding coefficient aij = 0. Let E(i)jk be the infimum of the exponents of x in the polynomial (Zi (x))jk , and E(i) := inf E(i)jk . From Section 2.1, it can be observed that E(0) = E(0)nn and E(1) = E(1)nn j,k

for even or odd n, with the equalities holding uniquely for (j, k) = (n, n) when n ≥ 2. From the above recursions, for i > 1, E(i) = E(i)nn and E(i) = 2E(i − 2) for both even and odd i, by considering the term (Zi−2 (x))nn H(i)nn (Zi−2 (x))nn in (3) and (4), and also 2(n−1) 2(n−1) σi−2 σi−1 (Zi−2 (x))n,n (Zi−2 (x))nn in (4). Let D(i) be a lower bound on the exponents of x in the polynomial di (x). By considering odd i in particular, where x is only an argument of Zi−2 , setting D(i) = E(i − 2) for all i is i−2 valid, where D(2) = D(3) = 2β −2n+2 = n2 −3n+2, and in general, D(i) = 2 2 (n2 −3n+2) i−3 and D(i) = 2 2 (n2 − 3n + 2) for all even and odd i > 1, respectively. These exponents are with respect to the argument x, whereas di is evaluated at σi2 , so by adding 2D(i) to the lower bound of e1 (i) + e2 (i), by considering m = n from Section 2.2, let ef˜(i) be a lower bound on the exponents of σi contained in f˜i for i > 1, where for even i, i−2

ef˜(i) = n − i − (i − 1)n(n − 1) + 2 ∗ 2 2 (n2 − 3n + 2) i 2

i 2

= (2 − i + 1)n2 + (i − 3 ∗ 2 )n + 2 7

i+2 2

− i,

(5)

and for odd i, i−3

ef˜(i) = n − i − (i − 2)n(n − 1) + 2 ∗ 2 2 (n2 − 3n + 2) i−1

i−1

(6)

i+1

= (2 2 − i + 2)n2 + (i − 3 ∗ 2 2 − 1)n + 2 2 − i. i−1

i

The function ef˜ can be further lower bounded by c2 2 n2 and c2 2 n2 for even and odd i using a constant 0 < c < 14 for n sufficiently large: For even i, ef˜(i) i

2 2 n2

= =

i

i

i

2 2 n2 i (i − 3 ∗ 2 2 )

i+2

(2 2 − i + 1)n2 + (i − 3 ∗ 2 2 )n + 2 2 − i i

(2 2 − i + 1)

+

i

22 1 3 ≥ − + 0, 4 n

i

22 n

i+2

+

2 2 −i i

2 2 n2

where the first term equals 41 when i = 4. For odd i, i−1

ef˜(i) i−1

2 2 n2

= ≥

i−1

i+1

(2 2 − i + 2)n2 + (i − 3 ∗ 2 2 − 1)n + 2 2 − i i−1

2 2 n2 1 3 − + 0, 4 n

i−1

i−1

2

where the bound c2 2 n2 < 14 2 2 n2 is also valid for i = 1: mn − 1 ≥ n2 − 1 > n4 when i−1 n ≥ 2. For both even and odd i, let nc ≥ 2 be sufficiently large such that ef˜(i) ≥ c2 2 n2 holds for i ≥ 1 and n ≥ nc for a given c > 0. Given that a lower bound on {tij }j for each i has been established, the next step is to upper bound P(1 − σn > ϵ) for any ϵ > 0. Defining σ0 := 1,

P(1 − σn > ϵ) = P(

n X

σi−1 − σi > ϵ)

i=1

≤ P(∪ni=1 (σi−1 − σi > ≤

n X

P(σi−1 − σi >

i=1

ϵ P(σi−1 − σi > ) = n

(7)

n

··· 0

fi (σ1 , ..., σi )dσi dσi−1 ...dσ1 0

0

Z σi−2 Z σi−1 − ϵ

Z 1

n

··· 0

ϵ ), where n

Z σi−2 Z σi−1 − ϵ

Z 1

=

ϵ )) n

0

fi (σi |σ1 , ..., σi−1 )dσi fi−1 (σ1 , ..., σi−1 )dσi−1 ...dσ1 , and

0

(8)

8

Z σi−1 − ϵ

n

Z σi−1 − ϵ

n

fi (σi |σ1 , ..., σi−1 )dσi ≤

sup

0

j

0

Z 1− =

ϵ σi−1 n

tij + 1 σi tij ( ) dσi σi−1 σi−1

sup(tij + 1)y tij dy j

0

Z 1− ϵ

n

≤

sup(tij + 1)y tij dy,

(9)

j

0

i where the first inequality uses (2), and the equality uses the change of variable y = σσi−1 . i−1 Given that the upper bound (9) is independent of {σj }j=1 , applying it to (8),

ϵ P(σi−1 − σi > ) ≤ n

Z 1− ϵ

n

···

sup(tij + 1)y dy j

0

Z 1− ϵ

n

=

Z σi−2

Z 1

tij

0

fi−1 (σ1 , ..., σi−1 )dσi−1 ...dσ1 0

sup(tij + 1)y tij dy.

(10)

j

0 i−1

For all j, tij ≥ ef˜(i) ≥ c2 2 n2 when n ≥ nc . It will now be shown that for sufficiently large i−1

i−1

2

n, sup(tij + 1)y tij ≤ (c2 2 n2 + 1)y c2 2 n for all y ∈ [0, 1 − nϵ ]. For a positive exponent t > 0 j

and y ∈ (0, 1 − nϵ ], d (t + 1)y t = y t (1 + (t + 1) ln(y)) ≤ 0 dt

(11)

1 1 when t ≥ −( ln(y) + 1), (11) holds for all y ∈ (0, 1 − nϵ ]. + 1). By enforcing t ≥ −( ln(1− ϵ ) 2

Considering t ≥ cn , i.e., i = 1 in c2

i−1 2

n

n , and n ≥ n1c := max{ ϵc1 , nc }, 2

ϵcn ≥ 1 ϵ ⇒ − (cn2 + 1) ≤ −1 n ϵ 2 ⇒ e− n (cn +1) ≤ e−1 ϵ 2 ⇒ (1 − )cn +1 ≤ e−1 n ϵ 2 ⇒(cn + 1) ln(1 − ) ≤ −1 n  2 ⇒ t ≥ cn ≥ −

 1 +1 , ln(1 − nϵ )

(12)

where the fourth inequality uses the bound (1 + x)k ≤ exk ,

(13)

which holds for k ≥ 0 and x ≥ −1. For n ≥ n1c , (11) holds for y ∈ (0, 1 − nϵ ] when t ≥ cn2 . i−1

i−1

2

2

It follows that sup(tij + 1)y tij ≤ (c2 2 n2 + 1)y c2 2 n ≤ (cn2 + 1)y cn for y ∈ [0, 1 − nϵ ] when j i−1 1 n ≥ nc , given that tij ≥ c2 2 n2 ≥ cn2 for all j and i ≥ 1.

9

Continuing from (10) assuming that n ≥ n1c , Z 1− ϵ i−1 n i−1 ϵ 2 P(σi−1 − σi > ) ≤ (c2 2 n2 + 1)y c2 2 n dy n 0 i−1  ϵ c2 2 n2 +1 = 1− n i−1

ϵ

2

≤ e− n (c2 2 n +1) √

ϵ

2

≤ e− n (c(1+(i−1)( 2−1))n +1) √

1

= e−ϵ(cn+ n ) e−ϵcn( 2−1)(i−1) ,

(14)

where the second inequality uses (13) and the third uses Bernoulli’s inequality, (1 + x)k ≥ √ i−1 1 + kx for nonnegative integer k and x > −1, to bound 2 2 ≥ 1 + (i − 1)( 2 − 1) for i ≥ 1. Plugging (14) into (7) for n ≥ n1c ,

P(1 − σn > ϵ) ≤ e

1 −ϵ(cn+ n )

n−1 X

√

e−ϵcn( 2−1)i

i=0 2

1 −ϵ(cn+ n )

=e ≤

√

1 − e−ϵcn ( 2−1) √ 1 − e−ϵcn( 2−1)

e−ϵcn √ . 1 − e−ϵcn( 2−1)

The series ∞ X

P(1 − σn > ϵ) ≤

n=n1c

∞ X

e−ϵcn √ −ϵcn( 2−1) 1 − e 1 n=n c

is convergent for all ϵ > 0 using the ratio test: √

e−ϵc(n+1) 1 − e−ϵcn( 2−1) √ ≤ e−ϵc < 1, e−ϵcn 1 − e−ϵc(n+1)( 2−1) hence σn → 1 almost surely using the first Borel-Cantelli lemma (Williams, 1991, Appendix A13). The previous theorem assumed that n → ∞, implying that m → ∞ as well. The following theorem considers the case where m → ∞ for a fixed value of n. Theorem 2. For matrices uniformly distributed in the unit spectral norm ball, ∆ ∼ U(B), with a fixed number of rows n, lim σn = 1 almost surely. m→∞

Proof. The proof begins by redefining ef˜(i) for i > 1 in equations (5) and (6) by not assuming that m = n. Using the exact value of e1 (i) found in Section 2.2, for even i, i

ef˜(i) = (m − n)(n − i + 1) + n − i − (i − 1)n(n − 1) + 2 2 (n2 − 3n + 2) =: m(n − i + 1) + z(i, n), 10

and for odd i, i−1

ef˜(i) = (m − n)(n − i + 1) + n − i − (i − 2)n(n − 1) + 2 2 (n2 − 3n + 2) =: m(n − i + 1) + z(i, n), where the terms in ef˜(i) not including m have been collected into the function z, and ef˜ is extended to the case i = 1 by defining z(1, n) := −1. Given that min z(i, n) ∈ Z is a fixed i

value independent of m, the function ef˜(i) ≥ m + min z(i, n) can be further lower bounded i by cm for a constant 0 < c < 1 when m ≥ mc for a sufficiently large mc ∈ Z. Wanting to again upper bound, in this case only over a finite number of terms, max(tij + j

n 1)y tij in (10), by considering m ≥ mnc := max{ ϵc , mc },

ϵ cm ≥ 1 n ϵ ⇒ − (cm + 1) ≤ −1 n ϵ ⇒ e− n (cm+1) ≤ e−1 ϵ ⇒ (1 − )cm+1 ≤ e−1 n ϵ ⇒(cm + 1) ln(1 − ) ≤ −1 n  ⇒ cm ≥ −

 1 +1 , ln(1 − nϵ )

using the same steps to show (12). It follows that (11) holds for y ∈ (0, 1 − nϵ ] when t ≥ cmnc , hence max(tij + 1)y tij ≤ (cm + 1)y cm for y ∈ [0, 1 − nϵ ] when m ≥ mnc , given that tij ≥ cm j

for all j and i ≥ 1. Continuing from (10) assuming that m ≥ mnc , Z 1− ϵ n ϵ P(σi−1 − σi > ) ≤ (cm + 1)y cm dy n 0  ϵ cm+1 = 1− n − nϵ (cm+1) ≤e . ϵ

Plugging (15) into (7), P(1 − σn > ϵ) ≤ ne− n (cm+1) . By applying the ratio test, ϵ

ϵ ne− n (c(m+1)+1) = e− n c < 1, − nϵ (cm+1) ne

it follows that the series ∞ X

P(1 − σn > ϵ) ≤

m=mn c

∞ X

ϵ

ne− n (cm+1)

m=mn c

is convergent for all ϵ > 0, proving that σn → 1 almost surely. 11

(15)

4

Computational Experiments

This section empirically verifies the convergence of σn → 1, demonstrates the inability to scale up Algorithm 2, and shows the viability of only using Algorithm 1 to approximately sample from U(B). Before presenting these results, implementation details are first described.

4.1

Implementation Details

All computation was done using 64-bit floating-point arithmetic using Numpy 2.4.4 in Python 3.13.13. All function evaluations when sampling σi are performed over a grid of M := 1000 evenly-spaced points X i in the interval [X1 , σ̂i−1 ], where the minimum value was fixed to X1 := 0.01: The need to take large negative exponents of σi results in overflows for small values, while at the same time, the likelihood of σi < X1 decreases rapidly as n increases. When sampling σi , which will be described next, the point X0 = 0 is also included, which still allows, through linear interpolation, sampled values of σi ∈ [0, X1 ]. Let the extended set of points be denoted as X̃ i := X i ∪ {X0 }. Sampling σi : Recalling that f˜i (·|σ1 , ..., σi−1 ) is the unnormalized distribution proportional to fi (·|σ1 , ..., σi−1 ), σi is sampled based on the inversion method (Devroye, 1986, Chapter 2): An unnormalized conditional CDF is first computed given the realized samples {σ̂j }i−1 j=1 , F̃i (Xji |σ̂1 , ..., σ̂i−1 ) :=

j X

f˜i (Xki |σ̂1 , ..., σ̂i−1 ),

k=1

with the conditional CDF of σi estimated as Fi (Xji |σ̂1 , ..., σ̂i−1 ) :=

F̃i (Xji |σ̂1 , ..., σ̂i−1 ) i F̃i (XM |σ̂1 , ..., σ̂i−1 )

over the points X i . The function Fi := Fi (·|σ̂1 , ..., σ̂i−1 ) is then extended to include X0 by setting Fi (X0 ) = 0. By sampling û ∼ U[0, 1], setting j ∗ = max(min{j : û ≤ Fi (X̃ji )}, 1), and using linear interpolation, σ̂i =

((X̃ji∗ − û)Fi (X̃ji∗ −1 ) + (û − X̃ji∗ −1 )Fi (X̃ji∗ )) (X̃ji∗ − X̃ji∗ −1 )

.

Computing Zi : Versus Algorithm 2, where Zi is updated in iteration i on line 16, given that Zi must be evaluated at the points X i+2 , which are not known until iteration i + 1 has completed, Zi is instead computed at the beginning of iteration i + 2. In order to delay underflows and overflows, if √ 0 < Zimax := maxi+2 |(Zi (x))jk | < 1 or Zimax > max value, j,k,x∈X

12

√ Zi is divided by Zimax or Zimax / max value, respectively, where max value√ is the maximum representable number in the 64-bit floating-point representation, with max value = 1.341 × 10154 . As a final detail, there were only two Zi tensors stored, Ze and Zo for even and odd indexed singular values, respectively, with one being updated at the beginning of each iteration. Exponent Canceling: The unnormalized distribution f˜i (·|σ1 , ..., σi−1 ) consists of the product e (i)+e2 (i) of σi 1 and |di (σi2 )|. As given by 2D(i) in the proof of Theorem 1, the exponents of σi 2 in |di (σi )| are non-negative and increasing with n and i. The exponent e1 (i) + e2 (i), found in Section 2.2, is always decreasing in i, and for values of m < 2n, e1 (i) + e2 (i) is negative: e1 (i) + e2 (i) ≤ (m − n)(n − 1) + n − 2 − n(n − 1) ≤ (n − 1)2 + n − 2 − n(n − 1) = −1, where the first inequality considers i = 2, and the second inequality uses m = 2n − 1. For values of σi < 1, as n increases, large negative exponents e1 (i) + e2 (i) cause overflows, whereas large positive exponents within |di (σi2 )| cause underflows, which can be remedied by shifting negative exponents from e1 (i) + e2 (i) to |di (σi2 )|. For i ∈ {2, 3}, e1 (i) + e2 (i) is initially increased by a quantity e3 (2) = e3 (3) > 0. To correct for this, within Z0 and Z1 , β is replaced by β − e34(2) , resulting in all of the polynomial exponents within Z0 (σ22 ) and Z1 (σ32 ), and subsequently within d2 (σ22 ) and d3 (σ32 ), to decrease by e3 (2). This initial injection of −e3 (2) decreases exponentially given the recursive updates of Zi on line 16 of Algorithm 2. At iteration i, the exponent of σi in di (σi2 ) has decreased by i i −e3 (2)2⌊ 2 ⌋−1 . To correct for this, for i > 3 e1 (i) + e2 (i) is increased by e3 (i) := e3 (2)2⌊ 2 ⌋−1 . The initial quantity e3 (2) is set to e3 (2) : =

min(max(−(e1 (n) + e2 (n)), 0), 2D(n)) . n 2⌊ 2 ⌋−1

(16)

When e1 (n) + e2 (n) < 0 and −(e1 (n) + e2 (n)) ≤ 2D(n), e1 (n) + e2 (n) + e3 (n) = 0, meaning that by the nth iteration the negative exponent e1 (n) + e2 (n) has been cancelled out with the positive exponents of σn within |dn (σn2 )|, and likewise, when −(e1 (n) + e2 (n)) > 2D(n), the positive exponents of σn within |dn (σn2 )| are cancelled out with the negative exponent e1 (n) + e2 (n). In addition to using exponent canceling, the convention that 0 ∗ ∞ = 0 is applied when computing f˜i (·|σ̂1 , ..., σ̂i−1 ) over X i by setting f˜i (x|σ̂1 , ..., σ̂i−1 ) = 0 when di (x2 ) = 0, regardless of xe1 (i)+e2 (i) overflowing or not.

4.2

Results

With the described implementation, for each choice of (n, m), n singular values were generated 1000 times. Beginning with square matrices, it was only possible to generate up to n = 19 singular values. The experiment with n = m = 20 was stopped after the computed conditional CDF underflowed for all x ∈ X 20 in the 99th trial. Fixing n = 19, singular values 13

were then generated for m ∈ {50, 100, 250, 500, 1000}. The mean and standard deviation of all sampled singular values are given in Table 1. In each trial, samples of U and V were also generated using Algorithm 1, with the computed errors n

m

1 XX |(Û Σ̂V̂ T − Û V̂ T )jk | and ave error := nm j=1 k=1

(17)

max error := max |(Û Σ̂V̂ T − Û V̂ T )jk | j,k

from replacing Σ̂ := diag(σ̂1 , .., σ̂n ) with an identity matrix, bypassing the use of Algorithm n

m

1

1

σ̂1 σ̂2 σ̂3 σ̂4 σ̂5 σ̂6 σ̂7 σ̂8 σ̂9 σ̂10 σ̂11 σ̂12 σ̂13 σ̂14 σ̂15 σ̂16 σ̂17 σ̂18 σ̂19 µ 0.509 s 0.287 µ 0.803 0.299 2 2 s 0.159 0.203 µ 0.902 0.604 0.210 3 3 s 0.087 0.168 0.146 µ 0.942 0.759 0.478 0.164 4 4 s 0.053 0.115 0.143 0.113 µ 0.962 0.840 0.644 0.394 0.052 5 5 s 0.035 0.081 0.111 0.125 0.070 µ 0.974 0.887 0.745 0.556 0.325 0.309 6 6 s 0.025 0.059 0.085 0.105 0.125 0.118 µ 0.980 0.916 0.809 0.664 0.489 0.472 0.452 7 7 s 0.019 0.044 0.066 0.085 0.101 0.097 0.092 µ 0.985 0.935 0.852 0.737 0.597 0.582 0.564 0.562 8 8 s 0.014 0.035 0.052 0.070 0.086 0.083 0.081 0.080 µ 0.988 0.948 0.882 0.789 0.675 0.662 0.649 0.647 0.645 9 9 s 0.011 0.028 0.042 0.057 0.072 0.070 0.071 0.071 0.071 µ 0.990 0.958 0.903 0.828 0.735 0.725 0.718 0.717 0.715 0.714 10 10 s 0.009 0.023 0.035 0.049 0.065 0.065 0.065 0.065 0.065 0.065 µ 0.992 0.965 0.920 0.884 0.839 0.834 0.828 0.827 0.825 0.824 0.823 11 11 s 0.008 0.019 0.029 0.050 0.078 0.078 0.078 0.078 0.077 0.077 0.077 µ 0.993 0.970 0.932 0.919 0.900 0.896 0.892 0.890 0.889 0.888 0.887 0.886 12 12 s 0.006 0.016 0.025 0.027 0.034 0.034 0.034 0.034 0.034 0.034 0.034 0.034 µ 0.994 0.975 0.942 0.931 0.917 0.913 0.910 0.908 0.907 0.906 0.905 0.904 0.904 13 13 s 0.006 0.014 0.021 0.023 0.026 0.026 0.026 0.026 0.026 0.026 0.026 0.026 0.026 µ 0.995 0.980 0.952 0.942 0.930 0.927 0.924 0.923 0.922 0.921 0.920 0.919 0.919 0.918 14 14 s 0.005 0.011 0.018 0.020 0.022 0.022 0.023 0.022 0.023 0.022 0.022 0.023 0.023 0.023 µ 0.996 0.982 0.960 0.951 0.942 0.939 0.936 0.935 0.934 0.933 0.933 0.932 0.931 0.931 0.930 15 15 s 0.004 0.010 0.015 0.017 0.019 0.019 0.019 0.019 0.019 0.019 0.019 0.019 0.019 0.019 0.019 µ 0.996 0.984 0.964 0.957 0.949 0.947 0.945 0.944 0.943 0.942 0.941 0.941 0.940 0.940 0.939 0.939 16 16 s 0.004 0.009 0.015 0.016 0.018 0.018 0.018 0.018 0.018 0.018 0.018 0.018 0.018 0.018 0.018 0.018 µ 0.997 0.987 0.972 0.966 0.959 0.957 0.954 0.953 0.952 0.952 0.951 0.95 0.95 0.949 0.949 0.948 0.948 17 17 s 0.003 0.008 0.011 0.012 0.014 0.014 0.014 0.014 0.014 0.014 0.014 0.014 0.014 0.014 0.014 0.014 0.014 µ 0.997 0.988 0.975 0.970 0.964 0.962 0.960 0.959 0.958 0.957 0.957 0.956 0.956 0.955 0.955 0.954 0.954 0.953 18 18 s 0.003 0.007 0.012 0.012 0.014 0.014 0.014 0.014 0.014 0.014 0.014 0.014 0.014 0.014 0.014 0.014 0.014 0.014 µ 0.997 0.989 0.976 0.972 0.966 0.965 0.963 0.962 0.961 0.961 0.960 0.959 0.959 0.958 0.958 0.958 0.957 0.957 0.956 19 19 s 0.003 0.006 0.009 0.01 0.011 0.011 0.011 0.011 0.011 0.011 0.011 0.011 0.011 0.011 0.011 0.011 0.011 0.011 0.011 µ 0.999 0.997 0.996 0.994 0.992 0.990 0.989 0.988 0.987 0.987 0.986 0.985 0.985 0.984 0.984 0.983 0.983 0.982 0.982 19 50 s 0.001 0.002 0.002 0.003 0.003 0.004 0.004 0.004 0.004 0.004 0.004 0.004 0.004 0.004 0.004 0.004 0.004 0.004 0.004 µ 0.999 0.998 0.998 0.997 0.996 0.995 0.994 0.993 0.992 0.991 0.991 0.990 0.989 0.989 0.988 0.988 0.987 0.987 0.986 19 100 s 0.001 0.001 0.001 0.001 0.002 0.002 0.002 0.002 0.002 0.002 0.002 0.002 0.002 0.003 0.003 0.003 0.003 0.003 0.003 µ 1.000 0.999 0.998 0.998 0.997 0.997 0.996 0.995 0.995 0.994 0.993 0.993 0.992 0.992 0.991 0.991 0.990 0.990 0.989 19 250 s 0.000 0.000 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.002 0.002 0.002 0.002 0.002 0.002 0.002 µ 1.000 0.999 0.999 0.998 0.998 0.997 0.996 0.996 0.995 0.995 0.994 0.994 0.993 0.993 0.992 0.992 0.991 0.991 0.990 19 500 s 0.000 0.000 0.000 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 µ 1.000 0.999 0.999 0.998 0.998 0.997 0.997 0.996 0.996 0.995 0.995 0.994 0.994 0.993 0.993 0.992 0.992 0.991 0.991 19 1000 s 0.000 0.000 0.000 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001

Table 1: The mean µ and standard deviation s from 1000 samples of the singular values {σi } of uniformly sampled n × m matrices from the unit spectral norm ball B using Algorithm 2. 14

2. The mean and standard deviation of these errors over the 1000 trials are given in Table 2. The values of e3 (2), e1 (n), e2 (n), and 2D(n) are also given, which were used for exponent canceling following (16). Examining Table 1, starting from n = m = 6, µ(σ̂n ) is monotonically increasing with µ(σ̂19 ) = 0.956 when n = m = 19. This is interpreted as verification that σn → 1 following Theorem 1. From the starting point of n = m = 19, by further increasing m, µ(σ̂n ) continues to monotonically increase to a final value of µ(σ̂n ) = 0.991 when m = 1000, verifying Theorem 2. n m e3 (2) e1 (n) e2 (n) 2D(n) 1

1

-

-

2

2

0

0

3

3

4

0

4

4

12

0

5

5

24

0

6

6

37.5

0

7

7

52.5

0

8

8

49

0

9

9

63

0

10 10 50.625

0

11 11 61.875

0

12 12 45.375

0

µ s µ -2 0 s µ -6 4 s µ -36 24 s µ -60 48 s µ -150 160 s µ -210 240 s µ -392 672 s µ -504 896 s µ -810 2,304 s µ -990 2,880 s µ -1452 7,040 s -

-

ave error max error 0.491 0.491 0.287 0.287 0.316 0.574 0.100 0.179 0.246 0.573 0.054 0.136 0.207 0.552 0.034 0.108 0.191 0.569 0.019 0.092 0.148 0.437 0.025 0.093 0.116 0.346 0.021 0.075 0.091 0.280 0.017 0.063 0.073 0.227 0.015 0.050 0.058 0.182 0.013 0.045 0.035 0.112 0.015 0.051 0.022 0.072 0.007 0.023

n

m

e3 (2) e1 (n) e2 (n) 2D(n)

µ s µ 14 14 36.969 0 -2366 19,968 s µ 15 15 42.656 0 -2730 23,296 s µ 16 16 28.125 0 -3600 53,760 s µ 17 17 31.875 0 -4080 61,440 s µ 18 18 20.320 0 -5202 139,264 s µ 19 19 22.711 0 -5814 156,672 s µ 19 50 22.590 31 -5814 156,672 s µ 19 100 22.395 81 -5814 156,672 s µ 19 250 21.809 231 -5814 156,672 s µ 19 500 20.832 481 -5814 156,672 s µ 19 1000 18.879 981 -5814 156,672 s 13

13

53.625

0

-1716

8,448

ave error max error 0.018 0.060 0.005 0.018 0.015 0.050 0.004 0.015 0.012 0.042 0.004 0.012 0.011 0.036 0.003 0.011 0.009 0.030 0.002 0.009 0.007 0.027 0.002 0.009 0.007 0.024 0.002 0.007 0.001 0.006 0.000 0.002 0.001 0.003 0.000 0.001 0.000 0.002 0.000 0.000 0.000 0.001 0.000 0.000 0.000 0.001 0.000 0.000

Table 2: For each matrix size (n > 1) × m, the value of e3 (2) used for exponent canceling following (16) is shown given the values of e1 (n) and e2 (n) in Algorithm 2 and 2D(n) from the proof of Theorem 1. The columns ave and max error contain the mean µ and standard deviation s from 1000 samples of computing these errors following (17). Examining Table 2, beginning again from n = m = 6, the means of ave and max error are monotonically decreasing with µ(ave error) = 0.007 and µ(max error) = 0.024 when n = m = 19. By further increasing m, these mean errors continue to decrease to 0.001 and 0.006 when m = 50, and finally to 0.000 and 0.001 when m = 1000, demonstrating ˆ The the viability in terms of accuracy of only using Algorithm 1 to generate samples ∆. remaining columns show the use of (16) in choosing e3 (2), and perhaps more interestingly, the scale of the exponents required when only sampling up to 19 singular values. The n lower bound on the exponents contained in |dn (σn2 )|, 2D(n) = O(n2 2 2 ), hence its value will continue to grow exponentially, demonstrating the futility of trying to sample {σi } for larger values of n. In addition, |e1 (n) + e2 (n)| = O(n3 ), limiting the benefit of exponent canceling as n increases, which can already be seen when n = m = 19, where e1 (n) + e2 (n) = −5, 814 is already much smaller than 2D(n) = 156, 672. As a final experiment, Algorithm 1 was tested on the example matrix size of 5, 120 × ˆ = Û V̂ T can be 25, 600 given in Section 1, and it was verified that approximate samples of ∆ 15

generated using numpy.linalg.qr (NumPy Developers, 2026) with 64 and 32-bit floating-point arithmetic.

5

Conclusion

This paper studied uniform sampling from the unit spectral norm ball in high dimension. Motivated by numerical challenges of directly sampling from this distribution, it was proven that as the size of the sampled matrix increases, all singular values converge to 1 almost surely. This finding led to a much simpler approximate sampling algorithm which is scalable and shown empirically to have monotonically decreasing error as the size of the matrix increases.

References Jeremy Bernstein and Laker Newhouse. Modular Duality in Deep Learning. In ICML, pages 3920–3930, 2025. Devansh Bisla, Jing Wang, and Anna Choromanska. Low-Pass Filtering SGD for Recovering Flat Optima in the Deep Learning Optimization Landscape. In AISTATS, pages 8299– 8339, 2022. Giuseppe Calafiore and Fabrizio Dabbene. A probabilistic framework for problems with real structured uncertainty in systems and control. Automatica, 38(8):1265–1276, 2002. Giuseppe C. Calafiore, Fabrizio Dabbene, and Roberto Tempo. Randomized Algorithms for Probabilistic Robustness with Real and Complex Structured Uncertainty. IEEE Transactions on Automatic Control, 45(12):2218–2235, 2000. Luc Devroye. Non-Uniform Random Variate Generation. Springer, 1986. John C. Duchi, Peter L. Bartlett, and Martin J. Wainwright. Randomized Smoothing for Stochastic Optimization. SIAM Journal on Optimization, 22(2):674–701, 2012. NumPy Developers. numpy.linalg.qr. generated/numpy.linalg.qr.html, 2026.

https://numpy.org/doc/2.4/reference/

Qwen Team. Qwen3 Technical Report, 2025. Andrey Tremba, Giuseppe Calafiore, Fabrizio Dabbene, Elena Gryazina, Boris Polyak, Pavel Shcherbakov, and Roberto Tempo. RACT: Randomized Algorithms Control Toolbox for MATLAB. In IFAC World Congress, pages 390–395. Elsevier, 2008. David Williams. Probability with Martingales. Cambridge University Press, 1991.

16