arXiv CSopen access

Beyond the $d^{2.5}$-mixing bound for Dikin walks on polytopes

Beyond the d2.5-mixing bound for Dikin walks on polytopes Yunbum Kook Georgia Tech [email protected]

arXiv:2607.13943v1 [cs.DS] 15 Jul 2026

Abstract Inspired by interior-point methods (IPM) for structured convex optimization, Kannan and Narayanan introduced the Dikin walk for sampling uniformly from polytopes in 2009. As in IPMs, the Dikin walk is affine-invariant, and its convergence is governed by the barrier geometry used to define its local proposal. They showed that the Dikin walk with the logarithmic barrier for a polytope in Rd with m linear inequalities mixes in md iterations. In 2017, Chen, Dwivedi, Wainwright, and Yu improved this to d2.5 using a Lewis-weight barrier, and conjectured that the correct mixing time should be d2 . We make progress toward this conjecture by improving the previous d2.5 -mixing bound. For exponential sampling over a polytope, we prove that the Dikin walk with a scaled Lee–Sidford metric mixes from a warm start in d2.25 iterations. This also yields an improved cold-start complexity via a known annealing framework. The main technical ingredient is improved average self-concordance of the Lee–Sidford metric, which gives high acceptance probability for the Metropolis filter along a random Dikin proposal. While previous analyses were effectively limited to second-order control due to technical difficulties, we develop a principled higher-order analysis. The proof combines a selective higher-order expansion of recursive bottleneck terms, a moving orthonormal-frame calculus for higher derivatives of the Lewis weights, and Wiener-chaos decompositions via multiple stochastic integrals to control the resulting Gaussian polynomials.

Contents 1 Introduction 1.1 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Technical overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1 4 5 9

2 Preliminaries 9 2.1 A mixing framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 The Lee–Sidford metric for linear inequalities . . . . . . . . . . . . . . . . . . . . . . 11 3 Average self-concordance for the LS metric 3.1 Proof sketch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Pathwise bound . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Concentration of base-point Gaussian polynomials . . . . . . . . . . . . . . . . . . .

12 12 14 18

A Self-concordance definitions

31

B Deferred computations B.1 Tensor computations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.2 Row-map calculus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.3 Elementary frame calculus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

32 32 33 33

Introduction

We study the complexity of sampling from a high-dimensional polytope. Throughout the paper, for A ∈ Rm×d and b ∈ Rm , let K := {x ∈ Rd : Ax ≥ b} be a bounded full-dimensional polytope with m linear inequalities. The goal of this problem is to output a sample whose law is close to the uniform distribution (more generally, an exponential distribution) over K. In theoretical terms, this problem has rich connections to convex optimization and convex-body sampling [DFK91, LS93, KLS97, LV06, CV18, JLLV26, KVZ26, KZ25], and on the applied side, it has found important applications in studying human metabolic networks in systems biology [Tea13, HCT+ 17, KLSV22, KLSV23]. This is an important and well-studied special case of the more general problem of log-concave sampling with only zeroth-order oracle access to the distribution [LV07, KV25b, KV25a, KV25c]. The best-known approaches currently have complexity that grows roughly quadratically in the dimension [KLS97, KVZ24, KVZ26] from a warm start (i.e., when already close to a target logconcave distribution) and sub-cubic [KV25a, KV25c] from a cold start, after an appropriate affine transformation called (near-)isotropic rounding [JLLV26]. Here, each step involves a zeroth-order query (i.e., checking whether a point x is in K), which for a polytope has arithmetic complexity O(md) (for checking all the inequalities)1 . A natural question is whether the structure of a polytope, and in particular its explicit inequality description, can be exploited in the design of the sampler, rather than only through oracle queries. General log-concave samplers treat the inequality description only as a black-box oracle. Hence, a more “structural” approach would leverage the additional information to accelerate mixing of polytope samplers. To motivate principled approaches to this question and highlight its deep connection to optimization theory, we first turn to the corresponding optimization problem. Sampling and optimization. Sampling and optimization can be related through the correspondence between (1) sampling from a Gibbs distribution dπ(x) ∝ e−V (x) dx and (2) minimizing the potential V : Rd → R ∪ {∞}. Indeed, one can use sampling to solve the optimization problem via simulated annealing [KGJV83, HJJ03, KV06]—sampling from the annealed distributions dπt ∝ e−V /t dx while gradually decreasing the temperature t > 0: as t becomes small, the mass of πt concentrates near the minimizers of V . Conversely, ideas from optimization have inspired the design and analysis of sampling algorithms. A prominent example is the connection between gradient flow and Langevin dynamics: the latter can be interpreted as a Wasserstein gradient flow over the space of probability measures [JKO98]. Under the above correspondence, the optimization counterpart of exponential sampling over a polytope K is linear programming (LP): minimizing a linear function over a set of linear inequalities. Since LP is a special case of convex optimization, one might first consider applying a general-purpose optimization method, such as gradient descent. However, such methods typically depend on the conditioning of the problem instance and may struggle on ill-conditioned problems. This brings us back to the structural question above: when solving LP, how can one exploit the explicit inequality description of the polytope, rather than treating it merely as a black-box convex body? 1

In the case of a polytope, by amortizing membership queries and using fast matrix multiplication, this can be reduced to mdc1 +0.5 for rounding and mdc2 for each sample thereafter, where c1 < 3.2 and c2 < 2.26 with current matrix multiplication bounds [JLLV26].

1

Interior-point methods (IPM). This question is precisely what led to the development of interior-point methods. In 1967, Dikin introduced an interior-point method based on the Dikin ellipsoid [Dik67], a local ellipsoidal approximation of the feasible region induced by the logarithmic barrier. Following the subsequent breakthroughs of Karmarkar [Kar84], Renegar [Ren88], and Nesterov and Nemirovski [NN94], the IPM literature has evolved into a general theory based on self-concordant barriers. In this theory, the iteration complexity is controlled by the barrier parameter. Roughly speaking, self-concordance imposes a regularity condition on the third-order derivative of the barrier, ensuring that its local Hessian geometry does not change too rapidly, while the barrier parameter quantifies the global quality of the barrier. For example, for a polytope K, P T the logarithmic barrier ϕ(x) = − m i=1 log(ai x − bi ) has barrier parameter ν = O(m). A crucial feature of IPM is affine-invariance. Unlike general-purpose optimizers, their performance is not affected by the conditioning of the problem instance, but by the geometry induced by the chosen barrier. Nesterov and Nemirovski showed that IPM based on a self-concordant barrier with parameter ν requires O(ν 1/2 ) iterations. Thus, using the logarithmic barrier gives an O(m1/2 )-iteration method for LP. A major direction in IPM theory has been to reduce this dependence on the number m of constraints, since m can be inflated by nearly redundant inequalities without changing the intrinsic geometry of the feasible region2 . This motivated the search for efficiently computable barriers with smaller √ parameters. In 1989, Vaidya’s breakthrough [Vai96] constructed such a barrier with parameter O( md), based on leverage scores of the constraint matrix. Finally in 2014, Lee and Sidford obtained an efficiently computable e barrier with parameter O(d) based on Lewis weights of the constraint matrix [LS14, LS19], which 1/2 e leads to O(d )-iteration IPM. Thus, in the optimization theory of LP, the right choice of a barrier removes the linear dependence on the number of inequalities. IPM for sampling: the Dikin walk. The connection above suggests that the IPM perspective can also be used to address the structural question for sampling. In the setting of uniform sampling from a polytope, which may be viewed as the “random” feasibility for linear inequalities, Kannan and Narayanan initiated an interior-point theory for sampling by proposing the Dikin walk [KN09, KN12] in 2009. This walk is based on the classical Dikin ellipsoid from convex optimization: given a local 2 metric g : int K → Sd++ and a radius r > 0, let Nrg (x) := N(x, rd g(x)−1 ) be the Gaussian proposal 2 with mean x and covariance matrix rd g(x)−1 . The Dikin ellipsoid of radius r at x is defined as Dg (x, r) := {y ∈ Rd : ∥y − x∥g(x) ≤ r} ,

∥v∥2g(x) = ∥v∥2x := v T g(x)v

for v ∈ Rd .

Starting from x0 ∼ π0 , the Dikin walk with target distribution π and local metric g repeats the following steps: for i = 0, 1, . . . , 1. Sample y ∼ Nrg (xi ). 2. If y ∈ / K, then xi+1 ← xi . If y ∈ K, set xi+1 ← y with probability π(y) det g(y) exp(− 2rd2 ∥y − xi ∥2g(y) ) o p

n

min 1,

π(xi ) det g(xi ) exp(− 2rd2 ∥y − xi ∥2g(xi ) ) p

,

(1.1)

and xi+1 ← xi otherwise. 2

In principle, this m-dependence can be avoided: the universal barrier [NN89] and the entropic barrier [BE15] are d-self-concordant [LY21, Che23], and hence yield O(d1/2 )-iteration IPM. However, these barriers are generally difficult to compute or implement efficiently.

2

Kannan and Narayanan used the Hessian of the logarithmic barrier for the local metric g for uniform sampling π ∝ 1K . This makes the walk affine-invariant, and therefore avoids the usual preprocessing step of putting the body into near-isotropic position. They showed that, from an M -warm start3 , the Dikin walk with the log-barrier mixes in e O(md log M ε ) iterations to total variation (TV) distance ε from the uniform distribution; see also [SV16] for a simplified proof. Moreover, each step of the log-barrier Dikin walk can be implee mented in amortized O(md) time [LLV20]. Subsequently, Narayanan extended this viewpoint to convex bodies equipped with self-concordant barriers, further developing the interpretation of the Dikin walk as a randomized IPM for sampling and optimization [Nar16]. e More recently, Kook and Vempala improved the iteration complexity to O(md log 1ε ) from a cold start (i.e., no warm start) [KV24]. A related recent direction is due to Jiang and Chen [JC25], who studied regularized Dikin walk for sampling log-concave distributions truncated to polytopes. Their e results include a soft-threshold Dikin walk with warm-start mixing O(md + κd) for log-concave and log-smooth targets with condition number κ. Toward dimension-square mixing. For polytopes with many constraints (i.e., m ≫ d), the Dikin walk may suffer from slow mixing. This mirrors the central question from the optimization IPM: can the linear dependence on m be improved by choosing a better barrier? There have been a series of works toward this goal, closely paralleling the development of improved barriers in optimization. Gustafson and Narayanan showed that by using a geometry related to the classical John ellipsoid, one can obtain an iteration complexity polynomial only in the dimension, though with a large exponent, namely d27 [GN23]. Subsequently, Chen, Dwivedi, Wainwright, and Yu showed that the Dikin walk with a variant of the Lewis-weights metric mixes e 2.5 log M ) steps from an M -warm start [CDWY18], where the Lewis-weights metric can be in O(d ε viewed as an efficiently implementable proxy for the John ellipsoid geometry. In particular, they left e 2 polylog M ). This is the main motivation of the the conjecture that the mixing time could be O(d ε present paper: can the iteration complexity of the Dikin walk be reduced from the log-barrier bound e e 2 polylog m)? O(md) to O(d There are several reasons why a dimension-square mixing bound is the natural target. First, in e the optimization IPM theory, the barrier parameter can be reduced from m to O(d) by replacing the log-barrier with the Lee–Sidford barrier. Second, for general-purpose samplers such as the Ball walk [KLS97] and In-and-Out [KVZ26, KV25b], a dimension-square rate arises naturally once the target distribution has been put in isotropic position. For example, the mixing bounds for the aforementioned samplers scale as d2 ∥cov π∥, where the largest eigenvalue ∥cov π∥ of the covariance matrix measures the geometric skewness of the target distribution. Hence, those samplers require (near-)isotropic rounding (i.e., ∥cov π∥ = O(1)), after which they mix in d2 steps from an O(1)-warm start. Since the Dikin walk is affine-invariant, it does not necessitate this global rounding step. Intuitively, the local metric used by the Dikin walk plays the role of implicitly rounding the body around the current point. From this perspective, the Dikin walk with a properly chosen metric should behave as though the target were locally isotropic, making d2 -mixing the fundamental target. Earlier attempts. Since then, there have been several attempts toward the conjectured d2 -mixing bound. The challenge is that we need delicate analysis of the change of local metric: in bounding the acceptance probability of the proposal y from the current point x, one must not only show that a 2 proposal is well-behaved in the forward ellipsoid centered at the current point (i.e., N(x, rd g(x)−1 )), 2 but also control the corresponding reverse proposal density at the proposal (i.e., N(y, rd g(y)−1 )). 3

An initial distribution µ is said to be M -warm with respect to π if µ(S) ≤ M π(S) for every measurable set S.

3

This is referred to as average self-concordance (ASC) in the terminology of [KV24]. This is where some previous claims of dimension-square mixing have gaps that cannot be readily fixed. In 2020, Laddha, Lee, and Vempala [LLV20] introduced extensions of the classical theory of self-concordance, including strong self-concordance and symmetry, as part of a general framework for analyzing the Dikin walk. Their work originally claimed a d2 -iteration bound, but this claim leaves open an important Metropolis-filter issue. Specifically, after drawing a proposal y from the forward Dikin ellipsoid centered at the current point x, the analysis also needs to verify that x lies in the Dikin ellipsoid centered at y with high probability. This step was initially missed, and a direct way of fixing this issue yields a d3 -bound [Vem26]. In 2024, [GKM+ 24] claimed d2 -mixing of the Dikin walk with Lewis weights. In the proof sketch of Lemma D.5, they assert that [SV16, Proposition 7] can be used similarly for the Lewis weights, but without the complete justification4 . This is exactly the ASC challenge underlying the Metropolis filter that prevents prior works from improving d2.5 -mixing of [CDWY18]. Despite several plausible approaches, a complete proof of d2 -mixing for a Lewis-weight-based Dikin walk has remained elusive.

1.1

Results

We make the first progress in nearly a decade toward the conjectured d2 -mixing bound for the weighted Dikin walk of Chen, Dwivedi, Wainwright, and Yu. We show that the Dikin walk with a e 2.25 ). smaller d1/4 scaling of the Lee–Sidford metric has the warm-start iteration complexity O(d Theorem 1.1 (Sampling from warm start). Let K = {x ∈ Rd : Ax ≥ b} be a full-dimensional bounded polytope, where A ∈ Rm×d has full column rank and no all-zero rows, and b ∈ Rm . Let g0 be the unscaled Lee–Sidford metric (2.2), L = Cd1/4 polylog(md) for a universal constant C > 0, and g := Lg0 . Consider the Dikin walk with local metric g, a sufficiently small Θ(1) radius r, an initial distribution π0 , and a target exponential distribution π over K. Let πn be the law of the n-th iterate of the 12 -lazy5 Dikin walk for n ∈ N. Then, for any ε > 0, we obtain χ2 (πn ∥ π) ≤ ε after e d2.25 polylog m log N =O

χ2 (π0 ∥ π)
ε

iterations.

This result assumes a warm start, and a natural next question is how efficiently such a warm start can be generated. One simple approach is to initialize from an easy distribution, such as the uniform distribution over a small ball contained in K. However, this typically introduces an extra factor of d in the iteration complexity. Thus, obtaining a sharp cold-start guarantee requires a more sophisticated approach. Such a framework was developed in [KV24]. Their approach constructs a sequence of annealing distributions and repeatedly uses a Dikin-walk sampler to move from one distribution to the next. In the polytope setting, the annealing path has the form proportional to e−ϕ/t · 1K , where ϕ is an appropriate self-concordant barrier and the temperature parameter t is gradually increased. The schedule is chosen so that consecutive distributions are sufficiently close, allowing the previous distribution to serve as an O(1)-warm start for the next phase. Using this framework with the 4 [GKM+ 24] uses closeness of Lewis weights at the base point x and a random proposal z, obtaining |∥z − x∥2z − ∥z − 2 x∥x | = O(1) ∥z − x∥2x ≲ r2 ∥h∥2 /d ≈ r2 . Since ASC needs the O(r2 /d)-scale, this closeness argument is not sufficient

for their claim. The log-barrier argument of [SV16] achieves this smaller scale through the Gaussian-polynomial concentration. For Lewis weights, however, the corresponding expansion contains wz and Taylor remainders along the random segment [x, x + ηh], so the similar polynomial argument does not apply since the estimated point wz has dependence on h. In fact, this is exactly the main challenge in showing ASC. 5 It stays at the current point (i.e., xi+1 ← xi ) with probability 1/2.

4

e 2.5 ) warm-start bound and obtained an O(d e 23/8 ) Lewis-weights metric, they recovered the known O(d cold-start iteration bound, giving the first sub-cubic iteration complexity for this problem. By plugging our improved warm-start guarantee into the annealing framework of [KV24], we obtain an improved cold-start bound of d41/16 ≈ d2.56 .

Corollary 1.2 (Sampling from cold start). Given ε > 0, there exists an algorithm that generates a e 41/16 polylog m log 1 ) iterations of the Dikin walk such that dTV (law X, π) ≤ ε. sample X using O(d ε Finally, we note that the above bounds are stated in terms of Dikin-walk iterations. The corresponding arithmetic complexities can be obtained by multiplying the iteration bounds by the per-step cost of computing Lewis weights, which contributes a factor of approximately mdω−1 per iteration, up to polylogarithmic factors, where ω < 2.37 is the matrix multiplication exponent.

1.2

Technical overview

We now summarize the technical ingredients of the main results. In this overview, we ignore any polylogarithmic factors in m. In the optimization IPM, the barrier parameter ν of a self-concordant barrier is sufficient to bound the iteration complexity of the IPM as O(ν 1/2 ). However, in the sampling counterpart, self-concordance and barrier parameters alone are not enough to provide a tight mixing guarantee of the Dikin walk. Mixing analysis through self-concordance. The framework of [KV24, Theorem 1] shows that if the local metric g : int K → Sd++ is strong self-concordant (SSC), lower trace self-concordant (LTSC), average self-concordant (ASC), and ν̄-symmetric (see §2.1), then the Dikin walk with Θ(1)2 radius, an initial distribution π0 , and a target exponential distribution π mixes in O(dν̄ log χ (πε0 ∥π) ) iterations to be ε-close to π in χ2 (Lemma 2.2). Recall the Lee–Sidford (LS) metric. For x ∈ int K, let Sx := Diag(Ax − b) ∈ Rm×m and Ax := Sx−1 A. For p ≍ polylog m, let Wx = Diag wx be the ℓp -Lewis weights of Ax , defined as 1/2−1/p wx = σ(Wx Ax ), where σ(B) denotes the leverage scores of B. The unscaled LS metric is 1−2/p e defined as g0 (x) := AT Ax . It was already shown in [LLV20] that g0 satisfies ν̄ = O(d), SSC, x Wx and LTSC. Hence, the only missing condition for applying the mixing framework is ASC. ASC of the Lee–Sidford metric. Recall that given a local metric g, the radius-r Dikin proposal 2 at x is z ∼ Nrg (x) := N(x, rd g(x)−1 ). The Metropolis filter then compares this forward proposal density with the reverse proposal density from z back to x, and a crucial part is to bound the acceptance probability (1.1) from below, by showing the closeness of the Gaussian proposal at x and z. This desired property is referred to as average self-concordance in [KV24]. Informally, ASC says that a random Dikin proposal does not significantly change its squared local length when the metric is evaluated at the current point x and the proposal z. Precisely, a metric g is ASC if for every ε > 0 there is a dimension-free radius rε > 0 such that for every x ∈ int K and every r ≤ rε , Pz∼Nrg (x) |∥z − x∥2g(z) − ∥z − x∥2g(x) | ≤

2εr2
≥ 1 − ε. d

If we can show ASC of the LS metric, then we would immediately obtain the d2 -mixing as desired. However, existing analyses [CDWY18, KV24, JC25] prove ASC of the LS metric only after scaling it by roughly d1/2 . Since scaling the metric by a factor L also scales the symmetry parameter e e 3/2 ), and then the general framework gives dν̄ = d2.5 by L, this changes ν̄ = O(d) into ν̄ = O(d 1/2 mixing bound. Namely, the extra d in the previous mixing bound comes from the amount of scaling needed to ensure constant-radius ASC. 5

In this paper, we improve this scaling from d1/2 to d1/4 , showing that the unscaled LS metric g0 satisfies the ASC condition at radius d−1/8 . If g = Lg0 , then a radius-r proposal for g is the same e 1/4 ) converts the radius-d−1/8 (in law) as a radius-r/L1/2 proposal for g0 . Hence, taking L = O(d ASC for g0 into constant-radius ASC for g. This scaling increases the symmetry parameter only e 5/4 ), and the mixing framework yields O(d e · d5/4 ) = O(d e 9/4 ) iterations. Thus, the proof of to O(d Theorem 1.1 is reduced to proving this sharper ASC estimate for the unscaled LS metric. Prior approaches to ASC. We now describe the ASC proof. By affine invariance, fix a base point x = 0 and normalize the unscaled LS metric so that g0 (x) = Id .√A proposal from the unscaled metric can be written as z = x + ηh, where h ∼ N(0, Id ) and η = r/ d. Define the path F (t) := hT g0 (x + th)h ,

0 ≤ t ≤ η.

Since F (0) = ∥z − x∥2g0 (x) /η 2 and F (η) = ∥z − x∥2g0 (z) /η 2 , ASC follows from a high-probability bound on |F (η) − F (0)| ≤ 2ε. A straightforward approach would be a direct Taylor expansion of F (η) around η = 0. For the Dikin walk with Lewis weights, one may use the first-order expansion, F (η) − F (0) = ηF ′ (t∗ ) for some t∗ ∈ (0, η). However, technical computations lead to supt≤η |F ′ (t)| ≲ d, and this requires a e 3 )-mixing of the Dikin walk). d-scaling of the metric (which yields O(d Chen et al. [CDWY18] expanded F up to the second-order terms, F (η) − F (0) = ηF ′ (0) + 2 η ′′ ′ 2 F (t∗ ). We note that at the base point x = 0, the first-order term F (0) is now a polynomial in Gaussian variables h. As already used in [SV16, CDWY18], the Gaussian polynomial P (h) concentrates around its L2 -norm (Eγ [P (h)2 ])1/2 (see Lemma 3.2), and this allows one to save an extra d1/2 -factor in F ′ (0). Indeed, after substantial technical computations, [CDWY18] managed to show that with high probability, |F ′ (0)| ≲ d1/2 and supt≤η |F ′′ (t)| ≲ d3/2 . In this case, the second-order contribution is η 2 d3/2 = r2 d1/2 . To make this O(ε), the radius r should be reduced to d−1/4 , which introduces the d1/2 -scaling of their metric for a constant-radius ASC. One could guess from this pattern that F (k) (0) ≲ dk/2 and supt≤η |F (k) (t)| ≲ d(k+1)/2 with high probability. To implement this idea, one could consider expanding further and attempt to obtain better bounds. However, a brute-force higher-order Taylor expansion quickly becomes unmanageable: F ′′ (0) already contains several complicated Gaussian polynomials involving derivatives of Lewis weights, while higher derivatives introduce many additional terms and require calculus beyond the known first-order Lewis-weight calculus in [LS19]. This suggests that we need a more organized approach and principled tools to handle this technical problem. Technical contributions in the improved ASC proof. The improvement is not obtained by simply differentiating the LS metric a few more times. The main difficulty is to make the higher-order expansion give the small number of bottleneck terms, while keeping both the pathwise Lewis-weight calculus and the base-point Gaussian-polynomial estimates tractable. Our proof introduces three organizing ideas. First, instead of Taylor expanding all derivatives of F , we isolate a recursive bottleneck chain Hk (t) and expand only this chain; all other differentiated terms are treated as controlled pathwise terms. Second, we formulate the higher-order LS calculus through a moving orthonormal frame for the column space of a half of the LS metric. This removes irrelevant rotations of the column space and converts derivatives of the Lewis-weight matrix into controlled frame and row-map estimates. Third, for the base-point terms, we replace direct high-moment Gaussian calculations by a Hermite-polynomial decomposition, implemented through multiple stochastic integrals. This gives a canonical orthogonal decomposition of the high-order Gaussian polynomials and reduces their L2 estimates to structured tensor-norm bounds.

6

(1) Isolation of bottleneck terms. The first idea is a selective expansion that avoids the combinatorial explosion of higher-order Taylor terms of F . Along the path xt = x + th, write At := Axt , Wt := Wxt , and let Nt denote the Lewis-weight derivative matrix from Lemma 2.5. Set α := 1 − p2 and β := α/2, and define (v ◦2 )i := vi2 for a vector v (and similarly for a matrix) and st := At h ,

ut := Wtβ st ,

1/2

vt := Wt

st ,

−1/2 ◦2 ut .

qt := Wt

With this notation, F (t) = ∥ut ∥2 . The Lewis-weight derivative formula gives a decomposition of the first derivative of the form F ′ (t) = −2βH1 (t) − 2E1 (t), where H1 (t) := qtT Nt vt and E1 (t) := ⟨st , Wtα s◦2 t ⟩. A simple but crucial observation is that not every term produced by differentiating F ′ is a (k−1) (j) bottleneck. To isolate the bottleneck terms, define Hk (t) := qtT Nt vt for k ≥ 1, where Nt denotes the j-th derivative of Nt along the path. Differentiating Hk gives Hk′ (t) = Gk+1 (t) + Hk+1 (t), where Gk+1 (t) collects the terms in which the derivative lands on qt or vt , while Hk+1 (t) is the (k−1) single term in which the derivative lands on Nt . The point of this decomposition is that the G-terms are controllable from the same pathwise estimates used at the previous level, whereas the new H-term contains the next derivative of Nt and is the only bottleneck term. Thus, we only expand H1 , H2 , H3 , H4 , rather than all terms appearing in the derivatives of F . The resulting bound is as follows: |F (η) − F (0)| ≲ η |E1 (0) + βH1 (0)| + η 2 sup |E1′ (t)| + sup |G2 (t)| + |H2 (0)|

t≤η

+η

3

t≤η

sup |G3 (t)| + |H3 (0)| + η sup |G4 (t)| + |H4 (t)| .

t≤η

4

(1.2)

t≤η

The proof then splits into two tasks: (1) pathwise bounds for the good terms E1′ , G2 , G3 , G4 , and H4 , and (2) Gaussian-polynomial estimates for the base-point terms E1 (0) and Hi (0) for i ∈ [3]. (2) Pathwise part. Recall preliminary calculus for Lewis weights in [LS19]: for cp := 1 − 2/p, −1/2 −1/2 ′ := DW [h] = − Diag(W 1/2 N W 1/2 s Λx := Wx −Px◦2 , Λ̄x := Wx Λ x Wx , they obtained Wx,h x x x x x,h ) (see Lemma 2.5 for details). For the pathwise part, we extend this computation to develop a higherorder calculus for Lewis weights in Lemma 3.3. The main difficulty in computing higher-order derivatives of Wx arises from Nx . Indeed, its derivative was computed together with several additional matrices in [LS19, Lemma 36 and 37], so a direct computation quickly becomes unwieldy. The second idea is to carry out this higher-order Lewis-weight calculus in a moving orthonormal bt := W β At and choose a smooth orthonormal frame Ut for col B bt satisfying frame. We work with B t T ′ the condition Ut Ut = 0 (see §B.3). This frame removes the arbitrary rotation of the basis inside the column space, streamlining computations for downstream calculus. 1/2 Letting ut,i and ϕt,i := (ut,i ⊗ ut,i )/wt,i be the i-th rows of Ut and Φt , respectively, we rewrite ′ ′′ ′′′ ′ ′′ ′′′ Λ̄t = I − Φt ΦT t . Then, bounds on Nt , Nt , and Nt reduce to bounds on Φt , Φt , and Φt , which ′ ′′ ′′′ ′ ′′ ′′′ in turn reduce to bound on Ut , Ut , Ut . Using the identities for Ut , Ut , Ut (Lemma B.2) together with row-wise estimates for the map u 7→ u ⊗ u/∥u∥ in §B.2, we prove the good-event bounds in Lemma 3.3. This is the step that turns the higher-order Lewis-weight calculus into a systematic set of matrix estimates. On this event, uniformly for all t ∈ [0, η], we control the relevant slack directions, Lewis-weight scores, and derivatives of Nt . Consequently, Lemma 3.4 gives |E1′ (t)|, |G2 (t)| ≲ d, |G3 (t)| ≲ d3/2 , |G4 (t)| ≲ d2 , and |H4 (t)| ≲ d5/2 . The recursive decomposition is important here: once Gk+1 is separated from Hk+1 in Hk′ , bounds on Gk+1 terms immediately follow from the same pathwise bounds on qt , qt′ , vt , vt′ , and 7

(k−1)

Nt . Thus, these terms are controlled inductively rather than by expanding all derivatives of F directly. The remaining term H4 is bounded by Cauchy–Schwarz as |H4 (t)| = |qtT Nt′′′ vt | ≤ ∥qt ∥ ∥Nt′′′ ∥ ∥vt ∥ ≲ d5/2 . This is the residual bottleneck of our argument. In (1.2), it contributes e −1/8 ). η 4 d5/2 = r4 d1/2 , which is exactly what limits our ASC radius to r = O(d (3) Base-point part. For the base-point terms, we exploit the fact that they are Gaussian polynomials in direction h, in order to save an extra d1/2 in their high-probability bound. Once their L2 -norms are bounded, the standard concentration inequality for Gaussian polynomials (Lemma 3.2) converts those L2 -estimates into high-probability bounds. The first base-point term is the cubic polynomial E1 (0)+βH1 (0). In previous analyses [CDWY18, KV24], controlling this cubic term was one of the most involved calculations. We streamline this computation by using the Gaussian Poincaré inequality, recovering the previous bound of ∥E1 (0) + βH1 (0)∥L2 ≲ d1/2 . Hence, the first term in (1.2) contributes ηd1/2 = O(r). The most technical part of this paper is to estimate the L2 -norms of the quartic and quintic terms H2 (0) and H3 (0). These are Gaussian polynomials of degrees four and five, so a direct L2 -norm computation would require expanding eighth and tenth Gaussian moments. Such an expansion produces too many tensor terms. It is not clear how to organize unwieldy computations by hand, and this is where we need a more principled tool to keep these computations under control. Our approach is to decompose them into orthogonal Hermite polynomials (equivalently, Wiener chaoses [Nua06, Nua19]) and then use orthogonality in estimating the L2 -norm. This replaces a non-canonical high-moment expansion by an orthogonal decomposition, so that downstream tasks are substantially streamlined. Then, the MSI (multiple stochastic integral) formalism provides the algebra needed for this decomposition. For a symmetric L2 -function f : Rn+ → R and Brownian motion (Bt )t≥0 , the MSI of f (see Definition 3.6) is defined as In (f ) := n!

Z ∞ Z tn 0

0

···

Z t2 0

f (t1 , . . . , tn ) dBt1 · · · dBtn .

Note that this is the usual Itô integral when n = 1. One can represent base Gaussian polynomials as MSIs, and then leverage known results on MSI (e.g., Theorem 3.8, Lemma 3.7, and (3.13)) to estimate the L2 -norm in a streamlined and organized way. Also, a known isometry between symmetric tensors and functions naturally induces the MSI of symmetric tensors (see (3.14)). Through the MSI isometry (Lemma 3.12), L2 -norms of those MSIs of tensors can be related to the Frobenius norm of those tensors. In summary, the L2 estimates of Gaussian polynomials are reduced to bounding the norm of induced tensors. We illustrate how this toolkit is used in bounding ∥H2 (0)∥2 . By simple algebra, one can P 1/2 rewrite it as H2 (0) = i wi (hT Ki h)(hT Si h) for some matrices Ki and Si . Using the identities hT Ki h = 1 + I2 (Ki ) and hT Si h = tr Si + I2 (Si ), we can represent H2 (0) in terms of tensor MSIs P 1/2 as H2 (0) = i wi {1 + I2 (Ki )}{tr Si + I2 (Si )}. With the product formula for MSIs (Lemma 3.7) and orthogonality (3.13), this decomposes into orthogonal terms as H2 (0) = I4 (T4 ) + I2 (T2 ) + T0 , where the Ti are induced tensors of degree i. By the MSI isometry (Lemma 3.12), the L2 -norm is reduced to bounding ∥T4 ∥F , ∥T2 ∥F , and |T0 |. The remaining estimates exploit the fact P ◦2 ⪯ W and Lewis-weight derivative bounds, yielding ∥H2 (0)∥L2 ≲ d. The quintic term H3 (0) is handled by the same strategy, with second-derivative estimates for Nt replacing the first-derivative estimates, and gives ∥H3 (0)∥L2 ≲ d3/2 . These two estimates contribute η 2 d = O(r2 ) and η 3 d3/2 = O(r3 ) in (1.2), respectively.

8

Combining the pathwise estimates with the base-point Gaussian-polynomial bounds gives, with high probability, e + r 2 + r 3 + r 4 d1/2 ). |F (η) − F (0)| = O(r Choosing r ≍ d−1/8 makes the RHS sufficiently small, proving ASC for the unscaled LS metric at radius d−1/8 . Scaling the metric by d1/4 then gives constant-radius ASC with symmetry parameter d5/4 , and applying the mixing theorem of [KV24] proves Theorem 1.1. Warm-start generation. We now sketch the proof of the cold-start guarantee in Corollary 1.2. No new mixing argument is needed for this step, and we simply use the better scaling of the LS metric and the annealing framework of [KV24] as a black box. Theorem 2 of [KV24] states that if the local metric g is SSC, LTSC, ASC, and ν̄-symmetric with a barrier parameter ν, then given ε > 0 and an exponential distribution π, their annealing e max(d, ν̄, ν 1/2 ν̄ 3/4 ) log 1 ) scheme can generate a sample X such that dTV (law X, π) ≤ ε using O(d ε iterations of the Dikin walk. Since we now know that the scaled LS metric d1/4 g0 (previously d1/2 g0 ) e 5/4 ), plugging this into the same annealing framework satisfies all these properties with ν, ν̄ = O(d 5/4 5/4 41/16 gives d (d ) = d .

1.3

Discussion

We close this section by suggesting a potential route toward the conjectured d2 -mixing. Recall (j) that the bottleneck term in the proof is Hj+1 (t) = qtT Nt vt . Our argument establishes pathwise control of this chain up to j = 3, giving the terminal bound |H4 (t)| ≲ d5/2 . Its contribution is e −1/8 ) and hence the d1/4 scaling of the metric. η 4 d5/2 = r4 d1/2 , which forces the ASC radius r = O(d (j) If one could extend the pathwise estimates to (1) |qtT Nt vt | ≲ d1+j/2 for 0 ≤ j ≤ k (namely, (j) ∥Nt ∥ ≲ dj/2 since ∥qt ∥, ∥vt ∥ ≲ d1/2 in Lemma 3.3), and (2) the base-point Gaussian estimates to (j) e 2+1/(k+1) ) ∥qxT Nx vx ∥L2 ≲ d(j+1)/2 for 0 ≤ j ≤ k − 1, then the mixing framework would yield O(d warm-start mixing. The present paper carries out this program for k = 3, giving d2+1/4 = d9/4 . Therefore, further progress toward the d2 -conjecture amounts to developing higher-order versions of the two estimates proved here. We view the present paper as isolating these two problems as the main remaining obstacles to the conjectured d2 -mixing for polytope sampling within our approach.

2

Preliminaries

We assume throughout that A ∈ Rm×d has full column rank and no all-zero rows. For x ∈ int K, let m×m and A := S −1 A ∈ Rm×d . ai be the i-th row of A and write Sx := Diag(aT x x i x − bi ) ∈ R

2.1

A mixing framework

Self-concordance for sampling. Recall the self-concordance terminology of [KV24, Definition 1.1]. The basic regularity condition is self-concordance (SC): for a local metric g, this means −2∥h∥g(x) g(x) ⪯ Dg(x)[h] ⪯ 2∥h∥g(x) g(x)

for x ∈ int K and h ∈ Rd .

The conditions SSC and LTSC are the stronger Frobenius and trace versions used in the mixing framework; for completeness their precise forms are recalled in §A. The Lee–Sidford metric is SSC and LTSC. The condition proved in this paper is ASC, so we spell it out explicitly.

9

Definition 2.1 (Average self-concordance, [KV24, Definition 1.1]). A matrix function g : int K → Sd++ is called average self-concordant (ASC) if for any ε > 0 there exists dimension-free rε > 0 such that, for every x ∈ int K, if r ≤ rε , then Pz∼Nrg (x) |∥z − x∥2g(z) − ∥z − x∥2g(x) | ≤

2εr2
≥ 1 − ε. d

In this condition, the event is understood to include z ∈ int K, so that g(z) is defined. Mixing analysis. We use the Dikin walk framework of [KV24]. In that framework, the local metric must satisfy SSC, LTSC, ASC, and the following symmetry condition. The symmetry parameter originates in [LLV20, Definition 2]; we call g ν̄-symmetric if √ (2.1) Dg (x, 1) ⊆ K ∩ (2x − K) ⊆ Dg (x, ν̄) for all x ∈ int K It holds that ν̄ = O(ν 2 ) in general. Lemma 2.2 (Mixing of Dikin walk [KV24, Theorem 1]). Suppose that a metric g is ν̄-symmetric, SSC, LTSC, and ASC. Let π be an exponential distribution on K, and initialize the Dikin walk at π0 . Then, for r = Θ(1), the walk outputs a sample X ∼ µ satisfying χ2 (µ ∥ π) ≤ ε in O dν̄ log

χ2 (π0 ∥ π)
ε

iterations.

Here we use the exponential-target specialization of [KV24, Theorem 1]. The one-step coupling √ and isoperimetry argument in [KV24] gives conductance Φ ≳ 1/ dν̄. Since the Dikin walk is lazy and reversible with respect to π, Cheeger’s inequality gives an L2 (π) spectral gap γ ≳ Φ2 . Thus, for f0 := dπ0 /dπ and the transition kernel P , there exists a universal constant c > 0 such that χ2 (π0 P t ∥ π) = ∥P t (f0 − 1)∥2L2 (π) ≤ (1 − γ)2t ∥f0 − 1∥2L2 (π) ≤ exp −

ct
2 χ (π0 ∥ π) . dν̄

2

Choosing t ≳ dν̄ log χ (πε0 ∥π) yields the displayed χ2 guarantee. We refer readers to [JC25] for a simpler mixing analysis. The final step of the reduction is the following simple scaling observation. Proposition 2.3. For ε > 0, let g0 be a metric satisfying the ASC condition with rε = r(ε) (possibly dimension-dependent). If g = Lg0 with L ≥ 1, then the metric g satisfies the ASC condition with rε′ := L1/2 rε . Moreover, if g0 is ν̄-symmetric, then g is Lν̄-symmetric. 2

Proof. Using the ASC of g0 , with respect to z ∼ N(x, rd g(x)−1 ), if r ≤ rε , |∥z − x∥2Lg0 (z) − ∥z − x∥2Lg0 (x) | = |L (∥z − x∥2g0 (z) − ∥z − x∥2g0 (x) )| ≤ L

2ε (r/L1/2 )2 2εr2 = d d

with probability at least 1 − ε. Hence, if r ≤ rε′ = L1/2 rε , then the ASC condition for g = Lg0 holds. The symmetry statement follows directly from ∥y − x∥2Lg0 (x) = L ∥y − x∥2g0 (x) .

10

2.2

The Lee–Sidford metric for linear inequalities

For p ≥ 2, the ℓp -Lewis weight of Ax [LS19] is the vector wx ∈ Rm satisfying wx = σ Diag(wx )1/2−1/p Ax ,

where σ(B) := Diag(B (B T B)−1 B T ) denotes the leverage scores of the matrix B. Equivalently, for Wx := Diag wx , 1

Wx = Diag Wx2

− p1

1− p2

Ax (AT x Wx

1

2 Ax )−1 AT x Wx

− p1

.

Throughout, we take p = Θ(polylog m) and suppress polylogarithmic factors in m, p, and ε−1 . We define the baseline Lee–Sidford (LS) metric [LLV20] as 1− p2

g0 (x) := AT x Wx

(2.2)

Ax .

The standard LS metric used in the mixing theorem is a fixed polylogarithmic multiple of g0 . This e scalar is suppressed in O(·) and does not affect the ASC calculation. It is known from [LLV20, Lemmas 4.2 and 4.3] that, with this standard polylogarithmic normalization, the LS metric is SSC and has convex log-determinant. The latter gives tr(g0 (x)−1 D2 g0 (x)[h, h]) ≥ tr(g0 (x)−1 Dg0 (x)[h]g0 (x)−1 Dg0 (x)[h]) ≥ 0 , which is the LTSC condition. Thus, the main task in this paper is to prove the ASC estimate for g0 at the largest possible radius before the final scaling. 2.2.1

Lewis-weight calculus

We use the following Lewis-weight identities and estimates throughout the ASC proof. The first lemma provides the basic size controls for slack directions and Lewis-weight derivatives. Since 0 ⪯ Wx ⪯ I and 1 − 2/p ∈ [0, 1], we have 1− p2

T AT x Wx Ax ⪯ Ax Wx

Ax = g0 (x) .

T Lemma 2.4 ([LS19]). Let Wx = Diag(wx (Ax )) ∈ Sm ++ be the ℓp -Lewis weights and g(x) = Ax Wx Ax the Lewis-weights metric, and h ∈ Rd . 1/2

2

x Ax )]i • (Lemma 26) maxi∈[m] [σ(W(w ≤ 2m p+2 . x )i

• (Lemma 33) ∥Ax h∥Wx = ∥h∥g(x) and ∥Ax h∥∞ ≤

√

1

2m p+2 ∥h∥g(x) .

′ ∥ • (Lemma 34) ∥Wx−1 wx,h Wx ≤ p ∥h∥g(x) .

The next identity is the main way derivatives of the Lewis weights enter the proof. Below, 2 ij = Aij for a matrix A.

A(2) (= A◦2 ) is the matrix defined as (A(2) )

Lemma 2.5 ([LS19], Lemma 24). The directional derivative of the ℓp -Lewis weight Wx in direction h ∈ Rd is 1 1 ′ Wx,h := DWx [h] = − Diag(Wx2 Nx Wx2 sx,h ) , where sx,h := Ax h, 1

Px := Wx2 (2)

− p1

1

2 Ax g0 (x)−1 AT x Wx

−1

−1

− p1

,

Px(2) := Px◦2 ,

cp := 1 − 2/p, Λx = Wx − Px , Λ̄x = Wx 2 Λx Wx 2 , and Nx = 2Λ̄x (I − cp Λ̄x )−1 . def

def

11

def

The matrices in Lemma 2.5 satisfy Px(2) , Λx ⪯ Wx ⪯ I .

(2.3)

We also use the following bounds on Nx and its first directional derivative. Lemma 2.6 ([LS19]). Let x ∈ int K and h ∈ Rd . For cp = 1 − 2/p with p > 2, let Λ̄x := −1

−1

−1

−1

(2)

Wx 2 Λx Wx 2 = I − Wx 2 Px Wx 2 , Nx = 2Λ̄x (I − cp Λ̄x )−1 and θx = AT x Wx Ax . def

• (Lemma 31) Nx is symmetric and 0 ⪯ Nx ⪯ pI. √ 1 ′ ∥ • (Lemma 34) ∥Wx−1 wx,h ∞ ≤ p( 2m p+2 + p/2) ∥h∥θx . 1

1

• (Lemma 37) ∥(I + Nx )− 2 DNx [h] (I + Nx )− 2 ∥ ≤ 4p5/2 ∥h∥θx . Finally, the closeness estimate below propagates coordinate bounds from the base point to nearby points on the proposal segment. Lemma 2.7 ([LS19], Lemma 35). In the same setting as above, let xt = x + t h, σt,i := aT i xt − bi , [wt,i ]α d m and wt = wxt . Define zt,α ∈ R by [zt,α ]i := dt log σt,i . Then, √ 1
∥zt,α ∥∞ ≤ 2(1 + |α|p)m p+2 + p |α| max(1, p/2) ∥h∥ATx Wxt Axt . t

3

Average self-concordance for the LS metric

The main estimate is the following ASC bound for the LS metric g0 . e −1/8 ) > 0 so that Proposition 3.1. Given ε > 0, choose Rε = Θ(d ε Rε + Rε2 + Rε4 d1/2 ≤ . polylog(m/ε)

Then, for any x ∈ int K, the following ASC estimate holds whenever 0 < r ≤ Rε : Pz∼Nrg (x) |∥z − x∥2g0 (z) − ∥z − x∥2g0 (x) | ≤ 0

3.1

2εr2
≥ 1 − ε. d

Proof sketch

The LS metric and ASC quantity are affine-invariant, so we may fix x = 0 and g0 (x) = Id . Let F (t) := hT g0 (x + th)h for 0 ≤ t ≤ η := r/d1/2 and h ∼ N(0, Id ) . Then, the ASC estimate follows from |F (η) − F (0)| ≤ 2ε with high probability. As in previous works, we expand F (t) at t = 0 and present high-probability bounds. To explain the proof idea, we describe approaches in prior works and then introduce notation used throughout in the computation. Prior approaches. Earlier works expand F (t) to second-order derivatives, and bound |F ′ (0)| ≲ d1/2 and |F ′′ (t)| ≲ d3/2 (w.h.p.), respectively. Since η = r/d1/2 , the second-order term η 2 F ′′ is of order d1/2 . This d1/2 -factor forces previous analyses to scale up their base Lewis metric by d1/2 , thereby incurring the previously best d2.5 -mixing of the Dikin walk. A natural idea is to expand the second-order term further and to attempt to bound |F ′′ (0)| ≲ d. The trouble is that F ′′ (0) already has 6 distinct terms, and bounding some of them requires intensive L2 -norm estimations of Gaussian polynomials. Even worse, additional expansion requires bounds on F ′′′ (t), but the number of distinct terms in F ′′′ dramatically increases, and Lewis-weight calculus was established up to second order in [LS19]. Taken together, these obstacles make it daunting to e 2.5 ) improve the mixing guarantee of the Dikin walk with Lewis weights, which has remained at O(d since [CDWY18, KV24, JC25]. 12

Refined approach. Define α := 1 − p2 and β := α2 = 12 − p1 . Along xt := x + th, define At := Axt , Wt := Wxt , Nt := Nxt , and st := At h ,

1/2

ut := Wtβ st ,

vt := Wt 1/2

−1/p

ρt := Wt

st ,

st ,

−1/2 ◦2 ut .

qt := Wt

1/2

By Lemma 2.5, Wt′ = − Diag(Wt Nt Wt st ). For ℓt := log wt (elementwise), we have ℓ′t = −1/2 −Wt Nt vt . We also define zt := βℓ′t − st and Zt := Diag zt . ′ With the notation above, we have F (t) = ∥ut ∥2 . Since s′t = −s◦2 t , we have ut = Zt ut . Therefore, ′ ◦2 T α ◦2 F ′ (t) = 2⟨ut , Zt ut ⟩ = 2β⟨u◦2 t , ℓt ⟩ − 2⟨ut , st ⟩ = −2βqt Nt vt − 2⟨st , Wt st ⟩ =: −2βH1 (t) − 2E1 (t) .

Unlike earlier works, we group terms in a better way, and isolate and expand only “bottleneck” terms Rη ′ ′ in F (t) that are not controlled directly on a high-probability event. Since F (η) − F (0) = 0 F (t) dt, we may want to bound H1 (t) and E1 (t) by O(η −1 ) = O(d1/2 )(w.h.p.), but, as shown below, the argument can only afford O(d) bounds, so in the first-order term, both E1 (t) and H1 (t) are problematic. This leads us to differentiate both terms further. To present our approach in an organized way, let us define (k−1)

Hk (t) := qtT Nt (j)

where Nt

vt

for k ≥ 1 ,

denotes the j-th derivative along the path xt = x + th. Then, we can write Hk′ (t) = Gk+1 (t) + Hk+1 (t) ,

(k−1)

Gk+1 (t) := (qt′ )T Nt

(k−1) ′ vt .

vt + qtT Nt

(3.1)

Additional expansion produces base-point Gaussian polynomials η (E1 (0) + H1 (0)) along with pathwise terms of order O(η 2 {E1′ (t) + H1′ (t)}). In contrast to E1 (t) and H1 (t), the base-point terms afford an additional d1/2 -reduction in its order through the following standard concentration estimate for Gaussian polynomials. Lemma 3.2 (Concentration of Gaussian polynomials). For d ≥ 1, let P : Rd → R be a polynomial of degree n. For any t ≥ (2e)n/2 , q

Ph∼N(0,Id ) |P (h)| ≥ t E[P (h)2 ] ≤ exp −

nt2/n
. 2e

Namely, E1 (0) and H1 (0) are bounded by their L2 -norm with high probability. As seen shortly, their L2 -norms are O(d1/2 ), so the first-order term η (E1 (0) + H1 (0)) is bounded as O(1). We next consider the second-order term η 2 (E1′ + H1′ ) = η 2 (E1′ + G2 + H2 ). Fortunately, E1′ (t) and G2 (t) are O(d) = O(η −2 ) in a good high-probability region, so these terms can be handled directly. However, the bottleneck in prior analyses is that H2 (t) was only bounded by O(d3/2 ). This extra d1/2 -factor forces prior analyses to scale up their base Lewis metric by d1/2 , thereby incurring the previously best d2.5 -mixing of the Dikin walk. We find that Hk (t) is indeed a bottleneck term in our argument, so expand only this term up to fourth-order. Similar to the first and second-order terms, the base-point Gaussian polynomial Hk+1 (0) is controlled in L2 and bounded by Lemma 3.2, and the pathwise term Gk+1 and the terminal term H4 are bounded deterministically on the good event in Lemma 3.4: |F (η) − F (0)| ≲ η |E1 (0) + βH1 (0)| + η 2 sup |E1′ (t)| + sup |G2 (t)| + |H2 (0)|

t≤η

13

t≤η

+ η 3 sup |G3 (t)| + |H3 (0)| + η 4 sup |G4 (t)| + |H4 (t)| .

t≤η

t≤η

(3.2)

The rest of the proof has two independent tasks: (1) pathwise bounds on the good event, and (2) L2 -norm estimates of the base-point terms E1 (0) + βH1 (0), H2 (0), and H3 (0). For the first task, we establish higher-order Lewis-weight calculus. As seen later in Lemma 3.4, isolation of Gk and ′ Hk in Hk−1 proves fruitful, since a high-probability bound for Gk follows inductively, rather than via ad hoc computation. Regarding the second task, L2 -norm estimation of Gaussian polynomials F ′ (0) was arguably the most technical part in prior works. We streamline this using the Gaussian Poincaré inequality. As for H2 (0) and H3 (0), these are Gaussian polynomials of degree 4 and 5, so we would need to compute the expectation of Gaussian polynomials of degree 8 and 10 for their L2 -norms. As far as we can tell, these are technically the most difficult parts and require principled approaches to go beyond the previous work. To this end, we work with a decomposition of H2 (0), H3 (0) into Hermite polynomials. Since Hermite polynomials are orthogonal in L2 (γ), their L2 -norm estimation can be simplified. Furthermore, this toolkit can be handled in a principled way via a multiple stochastic integral (MSI) (see §3.3.2), combined with Wiener chaos decomposition. Known results on MSI substantially streamline L2 -norm estimation, which were previously almost infeasible without it.

3.2

Pathwise bound

We establish higher-order calculus of Lewis weights as follows. Lemma 3.3 (Good event). In the affine normalization x = 0 and g0 (x) = Id , let h ∼ N(0, Id ), √ e −1/8 ), and η = r/ d. There exists a good event E of P(E) ≥ 1 − ε/20 xt = x + th, 0 < r ≤ Rε = Θ(d on which for all t ∈ [0, η], ∥ρt ∥∞ , ∥st ∥∞ ≲ 1 ∥ut ∥, ∥vt ∥, ∥Nt vt ∥, ∥qt ∥, ∥qt′ ∥, ∥vt′ ∥, ∥zt ∥∞ , ∥zt ∥Wt , ∥Nt′ ∥ ≲ d1/2 ∥zt′ ∥Wt , ∥vt′′ ∥, ∥Nt′′ ∥ ≲ d ∥zt′′ ∥Wt , ∥Nt′′′ ∥ ≲ d3/2 . −1/2

Proof. All estimates below are uniform for t ∈ [0, η]. Note that ρ0 = W0 U := W0β A0 . Then, U T U = Id and U U T = P with Pii = wi . Hence, −1/2

ρ0,i = wi

u0 , and define

(U h)i ∼ N(0, 1) .

Using the standard bound ∥h∥ ≲ d1/2 and a union bound, we have that with probability at least 1 − ε/20, ∥h∥ ≲ d1/2 , ∥ρ0 ∥∞ ≲ 1 after the hidden logarithmic factor in ≲ is chosen as a function of m and ε. Hereafter, we always work on this event. For L ≍ d1/4 polylog(md), we note that ∥xt − x∥g = ∥xt − x∥Lg0 ≤ rL1/2 · ∥h∥/d1/2 , so by reducing a constant scale of r, we can enforce ∥xt − x∥g ≤ 1 on this event. By the ν̄-symmetry (2.1), the Dikin ellipsoid of radius 1 with respect to g is contained K, so the path is contained in K. By direct computation, for σt,i := aT i xt − bi , d −1 log σt,i = −st,i , dt

d 1 log |ρt,i | = − ℓ′t,i − st,i . dt p 14

(3.3)

This derivative is bounded by ∥h∥AT Wt At (up to logarithmic factors) by Lemma 2.7, applied with t the exponents 0 and −1/p. Since ∥xt − x∥g0 (x) = t∥h∥ ≤ η∥h∥ ≲ r polylog(m/ε), the choice r ≤ Rε makes the last quantity a sufficiently small absolute constant. By the self-concordance of the standard LS metric (after absorbing some polylogarithmic scalar relating it to g0 into the suppressed logarithmic factors), uniformly for t ∈ [0, η], ∥h∥g0 (xt ) ≲ ∥h∥g0 (x) = ∥h∥ ≲ d1/2 . Since AT t Wt At ⪯ g0 (xt ), the RHS in Lemma 2.7 is a polylogarithmic multiple of ∥h∥g0 (xt ) , hence 1/2 e O(d ) on the whole interval [0, η]. Thus, uniformly over this interval, ∥ρt ∥∞ ≲ 1 and ∥st ∥∞ ≲ 1. Moreover, as zt = βℓ′t − st , the derivative of log ρt,i above implies that ∥ℓ′t ∥∞ , ∥zt ∥∞ ≲ d1/2 . Using these stability results, ∥ut ∥2 =

X

X

wt,i ρ2t,i ≲

i

wt,i = d ,

1/2

∥vt ∥ = ∥Wt

1/p

st ∥ = ∥Wt

ut ∥ ≤ ∥ut ∥ ≲ d1/2 .

i

By Lemma 2.6, ∥Nt ∥ ≲ 1 up to polylogarithmic factors, so ∥Nt vt ∥ ≲ d1/2 . Also, ∥qt ∥2 =

X

wt,i ρ4t,i ≤ ∥ρt ∥4∞ d ≲ d ,

1/2 1/2 ′ ℓt − Wt st ∥ = ∥−βNt vt − vt ∥ ≲ d1/2 .

∥zt ∥Wt = ∥βWt

i

Next, since

(3.4) 1 1 vt′ = ( ℓ′t − st ) ◦ vt = −st ◦ vt − (Nt vt ) ◦ st , 2 2

we have

(3.5)

∥vt′ ∥ ≲ ∥st ∥∞ (∥vt ∥ + ∥Nt vt ∥) ≲ d1/2 .

1/2 ′ 1/2 1/2 zt = − β2 Wt (ℓ′t )◦2 − βNt′ vt − βNt vt′ + Wt s◦2 t , the first term is bounded as

In Wt

1/2

∥Wt

1/2 ′ ℓt ∥ = ∥ℓ′t ∥∞ ∥Nt vt ∥ ≲ d .

(ℓ′t )◦2 ∥ ≤ ∥ℓ′t ∥∞ ∥Wt

(3.6)

The Nt′ vt term is O(d), since ∥Nt′ ∥ ≲ d1/2 due to Lemma 2.6-3 and self-concordance of g0 . The last two terms are O(d1/2 ) in norm. Thus, ∥zt′ ∥Wt ≲ d. −1/2 ◦2 Differentiating qt = Wt ut and using u′t = Zt ut = (βℓ′t − st ) ◦ ut 1 1 2
′ 1 −1/2 qt′ = − ℓ′t ◦ qt + 2Wt (ut ◦ u′t ) = − ℓ′t ◦ qt + 2β (ℓ′t ◦ qt ) − 2 (st ◦ qt ) = − ℓ ◦ qt − 2 (st ◦ qt ) . 2 2 2 p t Since ℓ′t ◦ qt = −(Nt vt ) ◦ ρ◦2 t , ∥qt′ ∥ ≲ ∥st ∥∞ ∥qt ∥ + ∥ρt ∥2∞ ∥Nt vt ∥ ≲ d1/2 . bt = W β At , choose a smooth orthonormal basis Ut of col B bt , and define Pt := Ut U T For Nt , let B t t (2) ⊥ ◦2 and Pt := I − Pt . For Pt := Pt , recall

Nt = 2Λ̄t (I − cp Λ̄t )−1 ,

−1/2

Λ̄t = I − Wt

(2)

−1/2

Pt Wt

. −1/2

Denote the i-th row of Ut by ut,i , and define the i-th row of Φt by ϕt,i := wt,i −1/2 (2) −1/2 Φt ΦT Pt Wt , so Λ̄t = I − Φt ΦT t = Wt t . By (2.3), ∥Φt ∥ ≤ 1.

15

(ut,i ⊗ ut,i ). Then

Writing Rt := (I − cp Λ̄t )−1 , we have Nt′ = 2Rt Λ̄′t Rt ,

Nt′′ = 2Rt Λ̄′′t Rt + 4cp Rt Λ̄′t Rt Λ̄′t Rt .

(3.7)

′ ′ Thus, ∥Nt′ ∥ ≲ ∥Λ̄′t ∥ and ∥Nt′′ ∥ ≲ ∥Λ̄′t ∥2 + ∥Λ̄′′t ∥. Also, since Λ̄t = I − Φt ΦT t , it holds that ∥Λ̄t ∥ ≲ ∥Φt ∥ ′′ ′′ ′ 2 and ∥Λ̄t ∥ ≲ ∥Φt ∥ + ∥Φt ∥ . Moreover, the chain rule yields

ϕ′i (t) = DR(ut,i )[u′t,i ] ,

ϕ′′i (t) = D2 R(ut,i )[u′t,i , u′t,i ] + DR(ut,i )[u′′t,i ] .

(3.8)

By Lemma B.1, applied row by row to ϕt,i = R(ut,i ) = R(ui (t)), ∥Φ′t ∥F ≲ ∥Ut′ ∥F ,

∥Φ′′t ∥2F ≲ ∥Ut′′ ∥2F +

X ∥u′t,i ∥4 i

wt,i

.

(3.9)

bt (due to u′ = Zt ut ) and bt = Ut Mt for invertible Mt , so equating B b ′ = Zt B We can write B t t −1 ′ ′ ′ ′ ′ b B = U Mt + Ut M leads to U = Zt Ut − Ut M M . Hence, t

t

t

t

t

t

Ut′ = Pt⊥ Ut′ + Pt Ut′ = Pt⊥ Zt Ut + Ut UtT Ut′ . bt such that U T U ′ = 0, By Lemma B.2, we can replace Ut by a “good” orthonormal basis for col B t t ′ ⊥ and this leads to Ut = Pt Zt Ut . Differentiating this identity gives

Ut′ = Pt⊥ Zt Ut ,

Ut′′ = −Pt′ Zt Ut + Pt⊥ Zt′ Ut + Pt⊥ Zt Ut′ .

(3.10)

Using Zt = Diag zt and xT (A ◦ B)y = tr(Diag(x)A Diag(y)B T ), ∥Ut′ ∥2F = ztT (Pt ◦ Pt⊥ )zt = ztT (Wt − Pt◦2 )zt ≤ ∥zt ∥2Wt ≲ d . For the second derivative in (3.10), note that by Pt′ = Ut′ UtT + Ut (Ut′ )T and ∥Ut ∥ ≤ 1, we have ∥Pt′ ∥ ≤ 2∥Ut′ ∥. Also, ∥Zt ∥ = ∥zt ∥∞ . Thus, ∥Pt′ Zt Ut ∥F ≤ ∥Pt′ ∥ ∥Zt Ut ∥F = ∥Pt′ ∥ ∥zt ∥Wt ≲ d1/2 · d1/2 = d , ∥Pt⊥ Zt′ Ut ∥F ≤ ∥Zt′ Ut ∥F = ∥zt′ ∥Wt ≲ d ,

(3.11)

∥Pt⊥ Zt Ut′ ∥F ≤ ∥Pt⊥ ∥ ∥Zt ∥ ∥Ut′ ∥F ≲ d . Hence, ∥Ut′′ ∥F ≲ d. T T 2 Finally, from u′t,i = eT i (I − Pt )Zt Ut = ei Zt Ut − ei Pt Zt Ut and ∥ut,i ∥ = wi , we obtain that 1/2

T 2 T 2 T ◦2 2 ∥eT i Zt Ut ∥ ≲ ∥zt ∥∞ wt,i . Also, ∥ei Pt Zt Ut ∥ = ei Pt Zt Pt Zt Pt ei ≤ tr(Eii Pt Zt Pt ) = ei Pt zt ≲ 2 ∥zt ∥∞ wt,i . Thus, 1/2 1/2 ∥u′t,i ∥ ≲ ∥zt ∥∞ wt,i ≲ d1/2 wt,i . (3.12)

Combining this with ∥Ut′ ∥2F ≲ d gives X ∥u′t,i ∥4 i

wt,i

≲d

∥u′t,i ∥2 = d ∥Ut′ ∥2F ≲ d2 .

X i

By (3.9), ∥Φ′t ∥F ≲ d1/2 and ∥Φ′′t ∥F ≲ d. The derivative identities (3.7) therefore give ∥Nt′ ∥ ≲ d1/2 and ∥Nt′′ ∥ ≲ d. 1/2 1/2 1/2 For the third-order estimates, differentiating Wt zt′ = − β2 Wt (ℓ′t )◦2 −βNt′ vt −βNt vt′ +Wt s◦2 t , β 1/2 1/2 ′′ 1/2 zt = −βWt (ℓ′t ◦ ℓ′′t ) − Wt (ℓ′t )◦3 − βNt′′ vt − 2βNt′ vt′ − βNt vt′′

Wt

4

16

+

1 ′ 1 ′ 1/2 1/2 1/2 ′ ℓt ◦ Wt s◦2 ℓ ◦ Wt zt′ . t + 2Wt (st ◦ st ) − 2 2 t

−1/2

1/2

1/2

Since ℓ′t = −Wt Nt vt , we have Wt ℓ′′t = 12 ℓ′t ◦ Nt vt − Nt′ vt − Nt vt′ , so ∥Wt ℓ′′t ∥ ≲ d. Hence, e 3/2 ). Using ∥ℓ′ ∥∞ , ∥vt ∥ ≲ d1/2 and ∥N ′′ ∥, ∥z ′ ∥W ≲ d with (3.6), the first term is bounded by O(d t t t t 1/2 1/2 1/2 we can bound the norms of Nt′′ vt , Wt (ℓ′t ◦ ℓ′′t ), Wt (ℓ′t )◦3 , and ℓ′t ◦ Wt zt′ by d3/2 . As for vt′′ , differentiating vt′ in (3.5) gives vt′′ = s◦2 t ◦ vt +

1 1 1 Nt vt − st ◦ vt′ + Nt′ vt + Nt vt′ . 2 2 2

Using ∥st ∥∞ , ∥Nt ∥ ≲ 1 and ∥vt ∥, ∥vt′ ∥, ∥Nt′ ∥ ≲ d1/2 , we have ∥vt′′ ∥ ≲ d, and ∥Nt vt′′ ∥ ≲ d. Lastly, by ′′ 3/2 . s′t = −s◦2 t and ∥st ∥∞ ≲ 1, the remaining terms can be bounded as O(d). Hence, ∥zt ∥Wt ≲ d ′′′ T ′ ′ T ′ T ′′ ′′ We now bound ∥Nt ∥. Recall that Pt = Ut Ut , Pt = Ut Ut + Ut (Ut ) , and Pt = Ut UtT + ′ 2Ut (Ut′ )T + Ut (Ut′′ )T , so ∥Pt′ ∥F ≲ d1/2 and ∥Pt′′ ∥F ≲ d. Then, from (3.10), Ut′′′ = −Pt′′ Zt Ut − 2Pt′ Zt′ Ut − 2Pt′ Zt Ut′ + Pt⊥ Zt′′ Ut + 2Pt⊥ Zt′ Ut′ + Pt⊥ Zt Ut′′ . Note the identities ∥Zt′ Ut ∥F = ∥zt′ ∥Wt ≲ d ,

∥Zt Ut ∥F = ∥zt ∥Wt ≲ d1/2 ,

∥Zt Ut′ ∥F ≤ ∥zt ∥∞ ∥Ut′ ∥F ≲ d ,

and recall ∥u′t,i ∥2 ≲ d wt,i . Using this, ∥Zt′ Ut′ ∥2F =

X

′ 2 (zt,i ) ∥u′t,i ∥2 ≲ d

X

i

′ 2 wt,i (zt,i ) = d ∥zt′ ∥2Wt ≲ d3 .

i

Using these bounds together with ∥Pt′ ∥F ≲ d1/2 , ∥Pt′′ ∥F ≲ d, ∥zt ∥∞ ≲ d1/2 , and ∥Ut′′ ∥F ≲ d, we obtain ∥Pt′′ Zt Ut ∥F ≤ ∥Pt′′ ∥ ∥Zt Ut ∥F ≲ d3/2 ,

∥Pt′ Zt′ Ut ∥F ≤ ∥Pt′ ∥ ∥Zt′ Ut ∥F ≲ d3/2 ,

∥Pt′ Zt Ut′ ∥F ≤ ∥Pt′ ∥ ∥Zt Ut′ ∥F ≲ d3/2 ,

∥Pt⊥ Zt′′ Ut ∥F ≤ ∥Zt′′ Ut ∥F = ∥zt′′ ∥Wt ≲ d3/2 ,

∥Pt⊥ Zt′ Ut′ ∥F ≤ ∥Zt′ Ut′ ∥F ≲ d3/2 ,

∥Pt⊥ Zt Ut′′ ∥F ≤ ∥zt ∥∞ ∥Ut′′ ∥F ≲ d3/2 .

Therefore, ∥Ut′′′ ∥F ≲ d3/2 . Differentiating (3.8), 3 ′ ′ ′ 2 ′′ ′ ′′′ ϕ′′′ i (t) = D R(ut,i )[ut,i , ut,i , ut,i ] + 3D R(ut,i )[ut,i , ut,i ] + DR(ut,i )[ut,i ] .

Using Lemma B.1,∥u′t,i ∥2 ≲ d wt,i ,

P ′ 2 ′ 2 ′′ 2 ′′ 2 2 i ∥ut,i ∥ = ∥Ut ∥F ≲ d, and i ∥ut,i ∥ = ∥Ut ∥F ≲ d ,

P

2 ′′′ 2 ∥Φ′′′ t ∥F ≲ ∥Ut ∥F +

X ∥u′t,i ∥6 i

2 wt,i

+

X ∥u′t,i ∥2 ∥u′′t,i ∥2 i

wt,i

≲ d3 .

′′′ ′′ ′ 3/2 . Differentiating (3.7), Finally, note that ∥Λ̄′′′ t ∥ ≲ ∥Φt ∥F + ∥Φt ∥F ∥Φt ∥F ≲ d ′ ′′ ′ 3 3/2 ∥Nt′′′ ∥ ≲ ∥Λ̄′′′ . t ∥ + ∥Λ̄t ∥ ∥Λ̄t ∥ + ∥Λ̄t ∥ ≲ d

This completes the proof. Using this lemma, we present pathwise bounds on the good event. 17

Lemma 3.4 (Pathwise bounds for good terms). On the good event E, sup |E1′ (t)|, sup |G2 (t)| ≲ d , t≤η

sup |G3 (t)| ≲ d3/2 ,

t≤η

t≤η

sup |G4 (t)| ≲ d2 , t≤η

sup |H4 (t)| ≲ d5/2 . t≤η

e 2 )-mixing The estimate of |H4 (t)| is the source of the new radius (but also a bottleneck for O(d 4 5/2 4 1/2 e −1/8 ). of the Dikin walk). In (3.2), it contributes η d = r d , which leads us to take Rε = Θ(d (k−1)

(k−1)

vt + qtT Nt vt′ , its pathwise estimates immediately follow from Proof. For Gk+1 (t) = (qt′ )T Nt Lemma 2.6, using the good-event estimates ∥qt ∥, ∥qt′ ∥, ∥vt ∥, ∥vt′ ∥ ≲ d1/2 , the deterministic bound ∥Nt ∥ ≲ 1 from Lemma 2.6, and ∥Nt′ ∥ ≲ d1/2 , ∥Nt′′ ∥ ≲ d, and ∥Nt′′′ ∥ ≲ d3/2 . As for H4 , by Cauchy–Schwarz, |H4 (t)| = |qtT Nt′′′ vt | ≤ ∥qt ∥ ∥Nt′′′ ∥ ∥vt ∥ ≲ d5/2 . ′ ◦2 Finally, differentiating E1 (t) = ⟨st , Wtα s◦2 t ⟩ and using st = −st , ◦2 α ◦2 α ◦3 E1′ (t) = α⟨st , (ℓ′t ◦ Wtα ) s◦2 t ⟩ − ⟨st , Wt st ⟩ − 2⟨st , Wt st ⟩ . 1/2 ′ ℓt = −Nt vt , we can bound the first term as follows:

Since ∥ut ∥ ≲ d1/2 and Wt X

α wt,i |ℓ′t,i | |st,i |3 ≤

i

X

α 2 wt,i st,i

X
1/2 X α ′ 2 4
1/2
1/2 wt,i (ℓt,i ) st,i = ∥ut ∥ wt,i (ℓ′t,i )2 ρ2t,i s2t,i

i

i

≲ d1/2

X

i

2
1/2

wt,i (ℓ′t,i )

= d1/2 ∥Nt vt ∥ ≲ d .

i α The last two terms are bounded by ∥st ∥2∞ sT t Wt st ≲ d.

3.3

Concentration of base-point Gaussian polynomials

At the base point, write W := Wx , w := wx , N := Nx , Λ̄ := Λ̄x , and note that U = Wxβ Ax ,

U T U = Id ,

P = UUT ,

B = W 1/p U = Wx1/2 Ax .

For uh := U h, we can write v0 = Bh = W 1/p uh , s0 = W −β U h, and q0 = W −1/2 (U h)◦2 . The three base-point terms in (3.2) are Gaussian polynomials of degrees three, four, and five: E1 (0) + βH1 (0) ,

H2 (0) = q0T N0′ [h]v0 ,

H3 (0) = q0T N0′′ [h, h]v0 .

The cubic term is handled by the original LS moment calculation and the Gaussian Poincaré inequality. The quartic and quintic terms will be handled via Wiener chaoses and some MSI calculations. 3.3.1

The cubic Gaussian polynomial

Lemma 3.5. ∥E1 (0) + βH1 (0)∥L2 ≲ d1/2 , so E1 (0) + βH1 (0) ≲ d1/2 with high probability. Proof. Define M := −W −β U − βW −1/2 N W 1/p U = −W −1/2 (I + βN )W 1/p U . Then, T −β E1 (0) + βH1 (0) = ⟨s0 , W0α s◦2 U h, (U h)◦2 ⟩ + β⟨W −1/2 (U h)◦2 , N W 1/p U h⟩ 0 ⟩ + βq0 N0 v0 = ⟨W

= ⟨(W −β U + βW −1/2 N W 1/p U ) h, (U h)◦2 ⟩ = −⟨M h, (U h)◦2 ⟩ := −T (h) .

18

Since T is odd, ET = 0. By the Gaussian Poincaré inequality, T 2 T ◦2 2 T 2 E[T 2 ] ≤ E[∥∇T ∥2 ] = E ∥M T u◦2 h + 2U (uh ◦ M h)∥ ≲ E[∥M uh ∥ ] + E[∥U (uh ◦ M h)∥ ] .





As for the second term, denote the i-th row of M by mi . Using E[(aT h)2 (bT h)2 ] = ∥a∥2 ∥b∥2 + 2 (aT b)2 for deterministic vectors a, b ∈ Rd and h ∼ N(0, Id ) (see [KV24, Proposition C.2]) with 2 ∥eT i U ∥ = wi , E[∥U T (uh ◦ M h)∥2 ] ≤ E[∥uh ◦ M h∥2 ] =

X

2 E[u2h,i (mT i h) ] ≤ 3

i

X

wi ∥mi ∥2 = 3 tr(M T W M ) .

i

As for the first term, define Cij := (M M T )ij and write 2 ◦2 T T ◦2 ∥M T u◦2 h ∥ = (uh ) M M uh =

X

Cij u2h,i u2h,j .

ij 2 Since uh = U h is a centered Gaussian with covariance E[uh uT h ] = P , E[uh,i ] = wi , and E[uh,i uh,j ] = Pij . Thus, 2 E[∥M T u◦2 h ∥ ]=

X

Cij E[u2h,i u2h,j ] =

ij

X

Cij wi wj + 2

ij

X

Cij Pij2 = wT M M T w + 2 tr(P ◦2 C)

ij

◦2

= ∥M w∥ + 2 tr(M P M ) ≤ ∥M w∥2 + 2 tr(M T W M ) . T

2

T

T

We now bound tr(M T W M ) and ∥M T w∥2 . As for the trace, note that W 1/2 M = −(I + βN )W 1/p U . Since N ⪯ pId and ∥W 1/p U ∥2F ≤ ∥U ∥2F = d, we have (ignoring any polylogarithmic factors) tr(M T W M ) ≲ d. As for ∥M T w∥2 , denoting C := (I + βN ) ≾ I, ∥M T w∥2 = wT W −1/2 CW 1/p U U T W 1/p CW −1/2 w ≲ wT W −1 w =

X

wi = d .

i

Putting all these together, E[T 2 ] ≲ d, and the concentration result in Lemma 3.2 proves the claim. 3.3.2

Multiple stochastic integral calculus

We now estimate the L2 -norm of H2 (0) and H3 (0), for which we would need to bound the expectation of Gaussian polynomials of degree 8 and 10 with numerous terms. This is practically infeasible to carry out. We need a more principled and organized approach to carry out this analysis, leveraging already well-developed theory. To this end, we use the finite-dimensional Wiener chaos decomposition of Gaussian polynomials. This is equivalently the Hermite expansion under Gaussian measure. Precisely, one may want to represent a polynomial with respect to some orthonormal polynomial basis with respect to Gaussian measure, and then its squared L2 -norm is the sum of the squared coefficients of the basis. Wiener chaos calculus. One can interpret Gaussian polynomials as multiple stochastic integrals (MSI), following Nualart [Nua06, §1.1]. We refer readers to [Nua19, §3] for a gentler introduction. Suppose that (Bt )t≥0 is a Brownian motion defined over (Ω, F, P), with F generated by (Bt )t≥0 . Let L2s (Rn+ ) be the space of symmetric L2 functions f : Rn+ → R. For a non-symmetric f : Rn+ → R, 1 P its symmetrization is defined as f˜(t1 , . . . , tn ) = n! σ f (tσ(1) , . . . , tσ(n) ) where the sum runs over all permutations of [n].

19

Definition 3.6 (Multiple stochastic integral (MSI)). For f ∈ L2s (Rn+ ), its order-n multiple stochastic integral is Z Z Z In (f ) := n!

∞

tn

t2

···

0

0

0

f (t1 , . . . , tn ) dBt1 · · · dBtn .

For non-symmetric f , we define In (f ) := In (f˜). When n = 1, this is the usual Itô integral I1 (f ) = R+ f (t) dBt . As seen later, one can think of In (f ) as base Gaussian polynomials of degree n. For n, m ∈ N, the standard isometry of multiple Itô integrals says that for f ∈ L2 (Rn+ ) and g ∈ L2 (Rm + ), R

E[In (f ) Im (g)] = 1n=m n! ⟨f˜, g̃⟩L2 (Rn+ ) .

(3.13)

n+m−2r 2 For f ∈ L2s (Rn+ ) and g ∈ L2s (Rm ) of order r ∈ {0, . . . , n∧m} + ), their contraction f ⊗r g ∈ L (R+ is defined as

(f ⊗r g)(t1 , . . . , tn−r , s1 , . . . , sm−r ) =

Z Rr+

f (t1 , . . . , tn−r , x1 , . . . , xr ) g(s1 , . . . , sm−r , x1 , . . . , xr ) dx1 · · · dxr .

The product formula is the basic rule that turns products of Gaussian polynomials into orthogonal chaos pieces. Lemma 3.7 (Product formula, [Nua06, Proposition 1.1.3]). For f ∈ L2s (Rn+ ) and g ∈ L2s (Rm + ), n∧m X

n In (f ) Im (g) = r! r r=0

!

!

m In+m−2r (f ⊗r g) . r

The next result, the Wiener chaos expansion, is the underlying principle behind our attempt. Theorem 3.8 (Wiener chaos expansion, [Nua06, Theorem 1.1.2]). Every F ∈ L2 (Ω) admits a P 2 n unique orthogonal expansion F = EF + ∞ n=1 In (fn ) for suitable fn ∈ Ls (R+ ). Let Hn be the L2 (Ω)-closure of the collection of all order-n MSI. The theorem says that 2 L (Ω) = ⊕∞ n=0 Hn . Thus, after decomposing a Gaussian polynomial into its Wiener-chaos pieces, its 2 L norm is the sum of the squared norms of those pieces. Tensor MSI. For the finite-dimensional Gaussian h ∼ N(0, Id ), choose the orthonormal family ei (t) = 1[i−1,i] (t) in L2 (R+ ). Then hi := I1 (ei ) = Bi − Bi−1 ∼ N(0, 1). A symmetric n-tensor T = (Ti1 ,...,in ) ∈ (Rd )⊗n can be identified with the function fT :=

d X

Ti1 ,...,in ei1 ⊗ · · · ⊗ ein ,

(3.14)

i1 ,...,in =1

where (ei1 ⊗ · · · ⊗ ein )(t1 , . . . , tn ) := nj=1 eij (tj ), and this isometry induces a tensor MSI as In (T ) := In (fT ). For a non-symmetric tensor T , we use the same notation for In (sym T ). We will use the following tensor forms of the MSI rules, deferring their proof to §B.1. Q

Lemma 3.9. For a symmetric matrix M , I2 (M ) = hT M h − tr M . Lemma 3.10. For a symmetric third-order tensor C ∈ (Rd )⊗3 , I3 (C) = C[h⊗3 ] − 3⟨tr C, h⟩, where P tr C ∈ Rd is defined as (tr C)a := db=1 Cabb . Lemma 3.11. For symmetric matrices M, N , I2 (M ) I2 (N ) = I4 (M ⊗ N ) + 4I2 sym(M N ) + 2 tr(M N ) .

Lemma 3.12. If M, N are tensors of orders n, m, then E[In (M ) Im (N )] = 0 for n ̸= m, and ∥In (M )∥2L2 (Ω) = n! ∥sym M ∥2F ≤ n! ∥M ∥2F . 20

3.3.3

The quartic Gaussian polynomial

For the quartic term, we introduce the notation Nk := DN (x)[ek ] for k ∈ [d]. Also, for each i ∈ [m], define matrices (Mi )kℓ := eT Si := sym Mi . i Nk Beℓ , Then,

eT i DN [h] Bh =

X

T T hk hℓ eT i Nk Beℓ = h Mi h = h Si h .

kℓ

Decomposition of H2 (0). We now decompose H2 (0) into its Wiener-chaos pieces. Define −1/2 −1/2 K := W −1/2 U = W −1/p Ax (so KK T = Wx Ux UxT Wx and (KK T )ii = 1), and denote each T row by ki (so ∥ki ∥ = 1). Let Ki := ki ki . By Lemma 3.9 and ∥ki ∥ = 1, I2 (Ki ) = hT Ki h − tr Ki = hT Ki h − 1 . −1/p

1/2

1/2

Recall that H2 (0) = q T N ′ [h]v for q = Wx (Wx Ax h)◦2 and v = Wx Ax h = Bx h. Hence, 1/2 1/2 ′ T qi = wi (kiT h)2 = wi (1 + I2 (Ki )) and eT i N [h]Bh = h Si h = tr Si + I2 (Si ). Using Lemma 3.11, H2 (0) =

X 1/2

=

X 1/2

=

X 1/2 

wi

1 + I2 (Ki ) tr Si + I2 (Si )

i

wi

I2 (Ki ) I2 (Si ) + tr Si I2 (Ki ) + I2 (Si ) + tr Si

i

wi

I4 (Ki ⊗ Si ) + 4I2 sym(Ki Si ) + 2 tr(Ki Si ) + tr Si I2 (Ki ) + I2 (Si ) + tr Si



i

=: F4 + F2 + F0 , where F4 = I4 (T4 )

for T4 :=

X 1/2

wi Ki ⊗ Si ,

i

F2 = I2 (C2 )

for C2 :=

X 1/2

wi

4 sym(Ki Si ) + tr(Si ) Ki + Si ,

i

F0 =

X 1/2

wi

2 tr(Ki Si ) + tr Si .

i

Since ∥H2 (0)∥2 = ∥F4 ∥2 + ∥F2 ∥2 + ∥F0 ∥2 due to the orthogonality of {In } (3.13), it suffices to estimate the L2 -norm of F4 , F2 , F0 . By Lemma 3.12, ∥F4 ∥2 = ∥I4 (T4 )∥2 ≤ 24 ∥T4 ∥2F and ∥F2 ∥2 = ∥I2 (C2 )∥2 = 2 ∥C2 ∥2F , so the remaining task is to bound the Frobenius norms of structured tensors such as T4 and C2 . Useful estimates.

We extract essential terms that repeatedly appear in the tensor analysis: E1 :=

d X

T1 :=

∥Nk B∥2F ,

X

(tr Si )2 .

i

k=1

We will show E1 ≲ d and T1 ≲ d2 . Lemma 3.13. It holds that Λ̄w1/2 = N w1/2 = 0 ,

Na w1/2 =

1 2 N Bea , 2 21

X 1/2

wi Si =

i

1 T 2 B N B. 2

Proof. Since P ◦2 1 = w, Λ̄w1/2 = (I − W −1/2 P ◦2 W −1/2 ) w1/2 = w1/2 − W −1/2 P ◦2 1 = 0 . Since (I − cp Λ̄) w1/2 = w1/2 , we also have (I − cp Λ̄)−1 w1/2 = w1/2 . Then, N w1/2 = 2Λ̄(1 − 1 cp Λ̄)−1 w1/2 = 2Λ̄w1/2 = 0. Moreover, using Dw[h] = −W 1/2 N W 2 Ax h (Lemma 2.5), ωa := Dw1/2 [ea ] =

1 −1/2 1 W Dw[ea ] = − N Bea . 2 2

(3.15)

Differentiating N w1/2 = 0 in direction ea gives Na w1/2 = −N ωa =

1 2 N Bea . 2

Using this identity, X 1/2

wi S i

 kℓ

i

= =

1 1 X 1/2 T wi (ei Nk Beℓ + eT (w1/2 )T Nk Beℓ + (w1/2 )T Nℓ Bek i Nℓ Bek ) = 2 i 2
1 1 T T 2 1 T T 2 1 1 T 2 ek B N Beℓ + eT ek B N Beℓ = (B T N 2 B)kℓ . ℓ B N Bek = 2 2 2 2 2

This completes the proof. Before proceeding, note that ∥B∥ ≤ 1 and ∥B∥2F ≤ d from B T B = U T W 2/p U ⪯ U T U = Id . Lemma 3.14. E1 =

Pd

2 k=1 ∥Nk B∥F ≲ d.

Proof. Fix a direction r ∈ Rd and consider a path xt = x + tr. Using the computations along the path as in Lemma 3.3, for z[r] = βDℓx [r] − Ax r, ∥DN [r]∥F ≲ ∥DΛ[r]∥F ≲ ∥DΦ[r]∥F ≲ ∥DU [r]∥F ≤ ∥z[r]∥W . Since ∥B∥ ≤ 1, ∥N ∥ ≲ 1, and W ⪯ I, ∥DN [r] B∥F ≤ ∥DN [r]∥F ∥B∥ ≲ ∥W 1/2 z[r]∥ = ∥(I + βN ) W 1/2 Ax r∥ ≲ ∥W

1/2

Ax r∥ ≤

q

α r T AT x Wx Ax r = ∥r∥ .

Substituting r = ek , we obtain ∥Nk B∥F ≲ 1, so squaring and summing over k gives d. Lemma 3.15. T1 =

2 2 i (tr Si ) ≲ d .

P

Proof. For Bk := Bek and v := Hence, T1 =

X i

(tr Si )2 =

X

k Nk Bk ∈ R

P

vi2 =

i

X

m , then the i-th coordinate of v is exactly v

2

Nk Bk

≤d

k

X k

where the last step is Lemma 3.14. Lemma 3.16.

2 k ∥Nk B∥F ≲

P

P T 2 2 i ∥Si ∥F ≲ d and i (ki Si ki ) ≲ d.

P

22

∥Nk Bk ∥2 ≤ d

X

∥Nk B∥2F ≲ d2 ,

k

i = tr Si .

Proof. For the Frobenius bound, by the definition of Si as the symmetrization of Mi , X

∥Si ∥2F ≤

∥Mi ∥2F =

2 (eT i Nk Beℓ ) =

X

X

i

i,k,ℓ

i

X

∥Nk B∥2F ≲ d

k

by Lemma 3.14. Next, using ∥ki ∥ = 1, kiT Si ki = kiT Mi ki =

X

T ′ eT i Nj Beℓ kij kiℓ = ei N [ki ]Bki

j,ℓ

=

X

kik eT i Nk Bki ≤

X

2 (eT i Nk Bki )

k

1/2

X

≤

2 ∥eT i Nk B∥

1/2

.

k

k

Squaring and summing over i gives X

kiT Si ki

2

≤

i

X

2 ∥eT i Nk B∥ =

i,k

X

∥Nk B∥2F ≲ d .

k

This completes the proof. Combining the computational lemmas above, we can readily bound ∥H2 (0)∥L2 .

Norm of H2 (0).

e Lemma 3.17. ∥H2 (0)∥L2 ≲ d, so H2 (0) = O(d) with high probability.

Proof. We estimate the three orthogonal chaos components F4 , F2 , F0 . By Lemma 3.12, ∥F2 ∥2L2 = 2∥C2 ∥2F .

∥F4 ∥2L2 ≤ 24∥T4 ∥2F , For F4 , ∥T4 ∥2F =

2

X 1/2

wi Ki ⊗ Si

=

F

i

X√

wi wj ⟨Ki , Kj ⟩⟨Si , Sj ⟩ .

i,j

Observe that √

wi wj ⟨Ki , Kj ⟩ =

√

wi wj (kiT kj )2 = √

Pij2
= W −1/2 P ◦2 W −1/2 ij . wi wj

(3.16)

√ Set G := W −1/2 P ◦2 W −1/2 . The preceding identity gives wi wj ⟨Ki , Kj ⟩ = Gij . Here we use the kernel inequality in the following form. If Hij := ⟨Ai , Aj ⟩ is any positive-semidefinite Gram matrix, then 0 ⪯ G ⪯ I implies X

Gij ⟨Ai , Aj ⟩ = tr(GH) ≤ tr H =

i,j

X

∥Ai ∥2 .

(3.17)

i

Applying this with Ai = Si and using the first estimate in Lemma 3.16 ∥T4 ∥2F ≤

X

∥Si ∥2F ≲ d .

i

Thus, ∥F4 ∥2L2 ≲ d. For F2 , recall C2 =

X 1/2

wi

i

1/2 (Si + 4 sym(Ki Si ) + tr(Si )Ki ). Note that Lemma 3.13 yields i wi

P

X

tr Si = tr

1/2



wi Si =

i

23

1 tr(B T N 2 B) ≲ tr(B T B) ≲ d . 2

(3.18)

For the first summand, we again use Lemma 3.13: X 1/2 1 wi Si = ∥B T N 2 B∥F ≲ tr(B T B) ≲ d . 2 F i For the second summand, by X 1/2

wi

sym(Ki Si )

F

i

≤

i wi ≲ d and the first estimate of Lemma 3.16,

P

X 1/2

wi ∥Ki Si ∥F ≤

X 1/2

wi ∥Si ∥F ≤

i

X

i

wi

1/2 X

i

∥Si ∥2F

1/2

≲ d.

i

For the trace summand, by (3.16) and Lemma 3.15, X 1/2

wi Ki tr Si

i

2 F

=

X 1/2 1/2

wi wj ⟨Ki , Kj ⟩ tr Si tr Sj ≤

i,j

(tr Si )2 ≲ d2 .

X i

Hence, ∥C2 ∥F ≲ d, and ∥F2 ∥2L2 ≲ d2 . P 1/2 For F0 , recall F0 = i wi (2 tr(Ki Si ) + tr Si ). Since tr(Ki Si ) = kiT Si ki , Lemma 3.16 leads to X 1/2

wi

tr(Ki Si ) ≤

i

X

wi

1/2 X

i

(kiT Si ki )2

1/2

≲ d.

i

The second term is bounded by O(d) as in (3.18). Hence, ∥F0 ∥2L2 = |F0 |2 ≲ d2 . Combining the orthogonal chaos pieces, ∥H2 (0)∥2L2 = ∥F4 ∥2L2 + ∥F2 ∥2L2 + ∥F0 ∥2L2 ≲ d2 . The high-probability statement follows from Lemma 3.2. 3.3.4

The quintic Gaussian polynomial

Define Nab := D2 Nx [ea , eb ] for a, b ∈ [d]. Also, for each row i, define the order-3 tensor Ci0 by Ci0 (a, b, ℓ) := eT i Nab Beℓ ,

Ci := sym Ci0 .

′′ With this notation, Ci [h⊗3 ] = eT i N0 [h, h]Bh. We also define τi := tr Ci and (τi )a :=

b Ci [a, b, b]. T T 2 As in the quartic term, note that (ki h) = h Ki h = 1 + I2 (Ki ). Also,

P

Decomposition of H3 (0). Ci [h⊗3 ] = I3 (Ci ) + 3I1 (τi ) by Lemma 3.10. Thus, H3 (0) =

X 1/2

wi

1 + I2 (Ki ) I3 (Ci ) + 3I1 (τi ) .

i

By Lemma 3.7, I2 (Ki )I3 (Ci ) = I5 (Ki ⊗Ci )+6I3 (Ki ⊗1 Ci )+6I1 (Ki ⊗2 Ci ) ,

I2 (Ki )I1 (τi ) = I3 (Ki ⊗τi )+2I1 (Ki ⊗1 τi ) .

This implies that 1 + I2 (Ki ) I3 (Ci ) + 3I1 (τi ) = I5 (Ki ⊗ Ci ) + I3 (Ci + 6Ki ⊗1 Ci + 3Ki ⊗ τi )+

+ I1 (3τi + 6Ki ⊗2 Ci + 6Ki ⊗1 τi ) . Hence, we can write H3 (0) = I5 (T5 ) + I3 (T3 ) + I1 (T1 ), where T5 :=

X 1/2

wi Ki ⊗ Ci ,

i

T3 :=

X 1/2

wi

(Ci + 6Ki ⊗1 Ci + 3Ki ⊗ τi ) ,

X 1/2

(3τi + 6Ki ⊗2 Ci + 6Ki ⊗1 τi ) .

i

T1 :=

wi

i

24

As in the analysis of H2 (0), we first focus on bounding useful terms

Useful estimates.

d X

E2 :=

T2 :=

∥Nab B∥2F ,

X

Sτ :=

∥tr Ci ∥2 ,

X 1/2

i

a,b=1

wi

tr Ci .

i

To this end, we need two computational lemmas. Lemma 3.18. Z2 :=

a,b ∥W

P

1/2 z ∥2 ≲ d2 and B := P ∥B ∥2 ≲ d2 . 2 ab a,b a,b

Proof. Let C0 := I + βN , Ba := Bea , and za := z[ea ]. Recall from (3.4) that W 1/2 za = −C0 W 1/2 Ax ea = −C0 Ba . Differentiating this identity in direction −eb , 1 − ℓb ◦ W 1/2 za − W 1/2 zab = βNb Ba + C0 Ba,b , 2 where Ba,b = DB[eb ]ea = Diag( 12 ℓb − sb ) Ba . Substituting this and W 1/2 za = −C0 Ba yields
1 1 −W 1/2 zab = βNb Ba + C0 ( ℓb − sb ) ◦ Ba − ℓb ◦ C0 Ba . 2 2

Note that ∥eb ∥AT W A ≤ ∥eb ∥g0 = 1. By Lemma 2.7 with (3.3), ∥ℓb ∥∞ , ∥sb ∥∞ ≲ 1, and recall ∥B∥2F ≲ d. Thus, X X 1 2 B2 = ∥Ba,b ∥2 = Diag( ℓb − sb ) B F ≲ d2 . 2 a,b b 2 k ∥Nk B∥F ≲ d (Lemma 3.14),

Since ∥Nb ∥ ≲ 1 (Lemma 2.6) and E1 =

P

Z2 =

X

X

X

a,b

b

b

∥W 1/2 zab ∥2 ≲

∥Nb B∥2F + B2 +

∥ℓb ∥2∞ ∥C0 B∥2F ≲ d2 .

This completes the proof. Lemma 3.19. It holds that X

X

∥Ua ∥2F ≲ d ,

a

XX

2 ∥eT i Ua ∥ ≲ d wi ,

a

i

2 ∥eT i Ua ∥ ≲ d ,

a

X

∥Uab ∥2F ≲ d2 .

a,b

Proof. By (3.10) and (3.11), Ua = P ⊥ Za U and ∥Ua ∥2F ≤ ∥za ∥2W . Using ∥sa ∥∞ , ∥ℓa ∥∞ ≲ 1 and ∥C0 ∥ ≲ 1, X X X ∥Ua ∥2F ≤ ∥za ∥2W = ∥C0 Ba ∥2 ≲ ∥B∥2F ≲ d . a

a

a

P P P 1/2 1/2 This implies that i a ∥eT U ∥2 = a ∥Ua ∥2F ≲ d. Next, as ∥eT ≲ wi from i Ua ∥ ≲ ∥za ∥∞ wi P iT a 2

(3.12), we obtain that a ∥ei Ua ∥ ≲ dwi . ⊥ := P ⊥ U . By Lemma B.3, As for the last item, define Uab ab

⊥ Uab = P ⊥ (P ⊥ Za U )b = −P ⊥ Pb Za U + P ⊥ Zab U + P ⊥ Za Ub .

Since ∥Pb ∥F ≲ ∥Ub ∥F from Pb = Ub U T + U UbT , by Lemma 3.18, ∥P ⊥ Pb Za U ∥2F ≲ d

X

X

a,b

b

X

∥P

a,b

∥Ub ∥2F ≲ d2 ,

⊥

Zab U ∥2F ≤

∥Zab U ∥2F =

X

X

a,b

a,b

25

∥W 1/2 zab ∥2 = Z2 ≲ d2 ,

∥P ⊥ Za Ub ∥2F ≲ d

Thus,

P

a,b ∥P

X

X

a,b

b

∥Ub ∥2F ≲ d2 .

⊥ U ∥2 ≲ d2 . By Lemma B.3 again, we can take a suitable U so that ab F

X

X

a,b

a,b

∥Uab ∥2F ≲

∥P ⊥ Uab ∥2F +

X

∥Ua ∥2F

2

≲ d2 .

a

This proves the lemma. We now bound the useful quantities. Lemma 3.20. E2 =

2 2 a,b ∥Nab B∥F ≲ d .

P

Proof. For ϕi = R(ui ) = ui ⊗ui /∥ui ∥, we can check that D2 ϕi [ea , eb ] = DR(ui )[ui,ab ]+D2 R(ui )[ui,a , ui,b ]. By Lemma B.1, X

∥D2 ϕi [ea , eb ]∥2 ≲

X

a,b

∥ui,ab ∥2 +

a,b

2 2 X 1 X 1 X T 2 ∥ui,a ∥2 = ∥eT ∥ei Ua ∥2 . i Uab ∥F + wi a wi a a,b

By Lemma 3.19, X

∥Φab ∥2F =

XX i

a,b

∥D2 ϕi [ea , eb ]∥2 ≲

∥Uab ∥2F +

X

a,b

X X a

i

a,b

2 2 ∥eT i Ua ∥ ≲ d .

d

T T Differentiating Λ̄ = I − ΦΦT gives Λ̄ab = −Φab ΦT − Φb ΦT a − Φa Φb − ΦΦab . Since ∥Φ∥ ≤ 1, the summation of the Frobenius norm of the second-derivative terms is bounded by O(d2 ). Due to (3.9), ∥Φa ∥ ≲ ∥Ua ∥. Hence,

X

X

a,b

a

2 ∥Φb ΦT a ∥F ≤

∥Φa ∥2F

2

≲

X

∥Ua ∥2F

2

≲ d2 .

a

Hence, a,b ∥Λ̄ab ∥2F ≲ d2 . Next, differentiating N = 2Λ̄R in ea and eb yields Nab = 2RΛ̄ab R+2cp RΛ̄a RΛ̄b R+2cp RΛ̄b RΛ̄a R. Then, X X X ∥RΛ̄ab RB∥2F ≤ ∥R∥2 ∥RB∥2 ∥Λ̄ab ∥2F ≲ ∥Λ̄ab ∥2F ≲ d2 . P

a,b

a,b

a,b

For the quadratic terms, using Na = 2RΛ̄a R and ∥R−1 ∥ ≤ 1, we have RΛ̄a RΛ̄b RB = 14 Na R−1 Nb B. Since ∥Na ∥ ≲ 1 (Lemma 2.6) and E1 ≲ d, X

∥Na R−1 Nb B∥2F ≤

a,b

Therefore, E2 =

∥Na R−1 ∥2

X

X

a

b

∥Nb B∥2F ≲ d2 .

2 2 a,b ∥Nab B∥F ≲ d .

P

Lemma 3.21. T2 =

i ∥tr Ci ∥

P

2 ≲ d3 .

T T Proof. Using Nab = Nba , we have Ci (a, b, c) = 13 (eT i Nab Bc + ei Nbc Ba + ei Nca Bb ). Thus, for τi = tr Ci ,

τi (a) =

X b

Ci (a, b, b) =

2 TX 1 X 2 1 ei Nab Bb + eT Nbb Ba =: eT ∆a ) . i i ( Γa + 3 3 3 3 b b 26

Hence, T2 =

i ∥τi ∥

P

2 ≲ P ∥Γ ∥2 + P ∥∆ ∥2 , where a a a a

X

∥Γa ∥2 ≤ d

X

a

X

X

a,b

∥∆a ∥2 =

∥Nab B∥2F ≤ dE2

a,b 2

X

a

Therefore, T2 =

∥Nab Bb ∥2 ≤ d

Nbb B

F

b

X

∥Nbb B∥2F ≤ dE2 .

≤d

b

2 3 i ∥τi ∥ ≲ d .

P

Lemma 3.22. ∥Sτ ∥2 ≲ d3 for Sτ =

1/2 tr Ci . i wi

P

Proof. Recall from the previous lemma that with Γa = b Nab Bb , ∆a = b Nbb Ba , and ω := w1/2 , 2 1
2X T 1X T (Sτ )a = ω T Γa + ∆a = Bb Nab ω + B Nbb ω . 3 3 3 b 3 b a P

P

Using Λ̄ω = 0, N ω = 0, and Na ω = −N ωa = 12 N 2 Ba (Lemma 3.13), and ωb = − 21 N Bb (3.15), 1 (Na N Bb + Nb N Ba + N Nb Ba + N 2 Ba,b ) . 2 Thus, ∥Sτ ∥2 is bounded, up to a universal constant, by the sum of the seven quantities Nab ω =

A1 :=

XX a

A3 := A5 :=

2

A2 :=

,

XX a

b

X

BaT

BbT N Nb Ba X

2

Nb N Bb

A4 :=

,

2

A6 :=

,

BaT

a

BbT Nb N Ba

X

N 2 Bb,b

2

2

,

b

XX a

b

X

BaT

BbT N 2 Ba,b

2

X

2

a

b

X

XX a

b

a

A7 :=

BbT Na N Bb

N Nb Bb

, ,

b

.

b

We now estimate these terms. We use b ∥Bb ∥2 = ∥B∥2F ≲ d, ∥N ∥, ∥Na ∥, ∥B T ∥, ∥B T N ∥, ∥B T N 2 ∥ ≲ P P 1, b ∥Nb B∥2F = E1 ≲ d, and B2 = a,b ∥Ba,b ∥2 ≲ d2 . Then, P

A1 ≤

X X

A4 ≤

X X

2

∥Na Bb ∥ ∥N Bb ∥

a

X

X

a

b

∥Na ∥2

≲

b

∥N Bb ∥ ∥N Ba,b ∥

a

2

≲

X X

A7 = B T

N 2 Bb,b

2

X

≲

b

∥Bb ∥2

a 2

b

X

Bb,b

b

≤d

2

∥Bb ∥2

≲ d3 ,

X



∥Ba,b ∥2 ≲ dB2 ≲ d3 ,

b

X

b

∥Bb,b ∥2 ≤ dB2 ≲ d3 .

b

For the remaining terms, A2 = ∥B T N

X

X

Nb Bb ∥2 ≲ ∥

b

A3 ≲

XX a

∥Bb ∥2

b

X

b

A5 ≲

X

A6 ≲

X

Nb N Bb

b

∥Nb Ba ∥2 ≲ d

X

∥Nb Bb ∥2 ≲ d

∥Nb B∥2F ≲ d2 ,

b

≤d

b

b

≤d

X

X

b

b

∥Nb N Bb ∥2 ≲ d ∥N Nb Bb ∥2 ≲ d

Combining these bounds, we prove the claim. 27

∥Bb ∥2 = d2 ,

X b

X

b 2

X b



b 2

N Nb Bb

Nb Bb ∥2 ≤ d

X

∥Bb ∥2 ≲ d2 , ∥Bb ∥2 ≲ d2 .

Norm of H3 (0).

We are now ready to bound ∥H3 (0)∥L2 .

e 3/2 ) with high probability. Lemma 3.23. ∥H3 (0)∥2L2 ≲ d3 , so H3 (0) = O(d

Proof. Recall that ∥H3 (0)∥2L2 ≲ ∥T5 ∥2F + ∥T3 ∥2F + ∥T1 ∥2 , where T5 = T3 =

X 1/2

(Ci + 6Ki ⊗1 Ci + 3Ki ⊗ τi ) ,

wi

T1 =

X 1/2

wi

i

P

1/2 i wi Ki ⊗ Ci and

(3τi + 6Ki ⊗2 Ci + 6Ki ⊗1 τi ) .

i

First of all, since Ci = sym Ci0 , we have i ∥Ci ∥2F ≤ i ∥Ci0 ∥2F = a,b ∥Nab B∥2F = E2 ≲ d2 . Following (3.17), X X ∥T5 ∥2F = Gij ⟨Ci , Cj ⟩ ≤ ∥Ci ∥2F ≤ d2 . P

P

ij

P

i

Recall that for a matrix K, a third-order tensor C, and a vector τ , (K⊗1 C)abc =

X

(K⊗2 C)a =

Kar Crbc ,

r

X

Krs Crsa ,

(K⊗τ )abc = Kab τc ,

(K⊗1 τ )a =

X

r,s

Kar τr .

r

Since Ki = ki kiT and ∥ki ∥ = 1, we have ∥Ki ⊗1 Ci ∥2F =

X

∥Ki ⊗2 Ci ∥ =

X

∥Ki (Ci ):bc ∥2 ≤

b,c 2

∥(Ci ):bc ∥2 = ∥Ci ∥2F ,

X b,c

tr Ki (Ci )::c ≤ 2

c

∥Ki ∥2F ∥(Ci )::c ∥2F ≤

∥(Ci )::c ∥2F = ∥Ci ∥2F .

X

X

c

c

For T3 , following (3.17) for the last line, 2

X 1/2

wi Ci

i

X 1/2

wi Ki ⊗1 Ci

i

X 1/2

wi Ki ⊗ τi

i

For T1 , note that ∥

F 2 F 2 F

≤

X

wi ·

X

X

i

i

i

≤

X

wi ·

X

=

X√

∥Ci ∥2F ≤ d

∥Ci ∥2F ≤ d3 ,

∥Ki ⊗1 Ci ∥2F ≤ d

X

wi wj ⟨Ki , Kj ⟩⟨τi , τj ⟩ ≤

X

i

∥Ci ∥2F ≤ d3 ,

i

∥τi ∥2 = T2 ≲ d3 .

ij

1/2 2 2 3 i wi τi ∥ = ∥Sτ ∥ ≲ d . For the remaining terms,

P

X 1/2

wi Ki ⊗2 Ci

2

≤

X

=

X 1/2

i

wi ·

X

∥Ki ⊗2 Ci ∥2 ≤ d

X i

i

X 1/2

wi Ki ⊗1 τi

i

2

∥Ci ∥2F ≲ d3 ,

wi Ki τi

i

2

=

X√

wi wj ⟨ki , kj ⟩⟨ki , τi ⟩⟨kj , τj ⟩ ≤

ij

X i

⟨ki , τi ⟩2 ≤

X

∥τi ∥2 ≲ d3 .

i

Putting all these together, ∥H3 (0)∥2L2 ≲ d2 + d3 + d3 ≲ d3 . Finally, the high-probability bound follows from Lemma 3.2. Proof of the main claims. Proof of Proposition 3.1. Combining (3.2) with Lemmas 3.4, 3.5, 3.17, and 3.23, and choosing Rε suitably so that the union of the exceptional events has probability at most ε, we get e 1/2 ) + η 2 O(d) e e 3/2 ) + η 4 O(d e 5/2 ) = O(r e + r 2 + r 3 + r 4 d1/2 ) . |F (η) − F (0)| ≲ η O(d + η 3 O(d

By the defining condition on Rε , with hidden constants small enough, the RHS is at most 2ε whenever r ≤ Rε . Since ∥y − x∥2g0 (y) − ∥y − x∥2g0 (x) = η 2 (F (η) − F (0)), this proves the claim. 28

e −1/8 ). By Proof of Theorem 1.1. By Proposition 3.1, the baseline metric satisfies ASC for Rε = Θ(d e 1/4 ) g0 while ensuring ASC. Proposition 2.3, we can take r = O(1) and g⋆ = O(d Recall that the unscaled LS metric is SSC, has convex log determinant, and is O(d log3 m)e 1/4 ) multiplies ν̄ by O(d e 1/4 ) symmetric [LLV20, Lemmas 4.2 and 4.3]. Scaling the metric by O(d e 1+1/4 ) and satisfies SSC while preserving SSC and LTSC. Thus, g⋆ has symmetry parameter O(d and LTSC. Lemma 2.2 gives the stated mixing time.

Acknowledgement. We are deeply grateful to Santosh Vempala for discussions, encouragement, and comments on an earlier draft of this paper, and to Yuansi Chen and Minhui Jiang for their feedback and discussions about future directions. This work was supported in part by NSF Award CCF-2106444. We also acknowledge the use of ChatGPT to assist with proofreading an earlier draft. In particular, it helped with notation choices, computational simplifications in Lemma 3.3, and the orthonormal-frame calculus in Appendix B.

References [AMS08]

Pierre-Antoine Absil, Robert Mahony, and Rodolphe Sepulchre. Optimization Algorithms on Matrix Manifolds. Princeton University Press, Princeton, NJ, 2008.

[BE15]

Sébastien Bubeck and Ronen Eldan. The entropic barrier: a simple and optimal universal self-concordant barrier. In Conference on Learning Theory, volume 40, pages 279–279, Paris, France, 03–06 Jul 2015. PMLR.

[CDWY18]

Yuansi Chen, Raaz Dwivedi, Martin J. Wainwright, and Bin Yu. Fast MCMC sampling algorithms on polytopes. Journal of Machine Learning Research, 19(55):1–86, 2018.

[Che23]

Sinho Chewi. The entropic barrier is n-self-concordant, pages 209–222. Springer International Publishing, 2023.

[CV18]

Ben Cousins and Santosh S. Vempala. Gaussian cooling and O∗ (n3 ) algorithms for volume and Gaussian volume. SIAM Journal on Computing, 47(3):1237–1273, 2018.

[DFK91]

Martin Dyer, Alan Frieze, and Ravi Kannan. A random polynomial-time algorithm for approximating the volume of convex bodies. Journal of the ACM, 38(1):1–17, 1991.

[Dik67]

Iliya Iosiphovich Dikin. Iterative solution of problems of linear and quadratic programming. In Soviet Math. Dokl., volume 8, pages 674–675, 1967.

[GKM+ 24]

Yuzhou Gu, Nikki Lijing Kuang, Yi-An Ma, Zhao Song, and Lichen Zhang. Log-concave sampling from a convex body with a barrier: a robust and unified Dikin walk. In Advances in Neural Information Processing Systems, volume 37, pages 69230–69298. Curran Associates, Inc., 2024.

[GN23]

Adam Gustafson and Hariharan Narayanan. John’s walk. Advances in Applied Probability, 55(2):473–491, 2023.

[HCT+ 17]

Hulda S Haraldsdóttir, Ben Cousins, Ines Thiele, Ronan M T Fleming, and Santosh S. Vempala. CHRR: coordinate hit-and-run with rounding for uniform sampling of constraint-based models. Bioinformatics, 33(11):1741–1743, 01 2017.

[HJJ03]

Darrall Henderson, Sheldon H. Jacobson, and Alan W. Johnson. The theory and practice of simulated annealing, pages 287–319. Kluwer Academic Publishers, 2003.

[JC25]

Minhui Jiang and Yuansi Chen. Regularized Dikin walks for sampling truncated logconcave measures, mixed isoperimetry and beyond worst-case analysis. In Conference on Learning Theory, volume 291 of Proceedings of Machine Learning Research, pages 3017–3078. PMLR, 30 Jun–04 Jul 2025.

[JKO98]

Richard Jordan, David Kinderlehrer, and Felix Otto. The variational formulation of the Fokker-Planck equation. SIAM Journal on Mathematical Analysis, 29(1):1–17, 1998.

29

[JLLV26]

He Jia, Aditi Laddha, Yin Tat Lee, and Santosh Vempala. Reducing isotropy and volume to kls: Faster rounding and volume algorithms. J. ACM, 73(2), April 2026.

[Kar84]

Narendra Karmarkar. A new polynomial-time algorithm for linear programming. In Symposium on Theory of Computing, STOC ’84, pages 302–311. ACM Press, 1984.

[KGJV83]

Scott Kirkpatrick, C. Daniel Gelatt Jr, and Mario P. Vecchi. Optimization by simulated annealing. Science, 220(4598):671–680, 1983.

[KLS97]

Ravi Kannan, László Lovász, and Miklós Simonovits. Random walks and an O∗ (n5 ) volume algorithm for convex bodies. Random Structures & Algorithms, 11(1):1–50, 1997.

[KLSV22]

Yunbum Kook, Yin Tat Lee, Ruoqi Shen, and Santosh S. Vempala. Sampling with Riemannian Hamiltonian Monte Carlo in a constrained space. In Advances in Neural Information Processing Systems, volume 35, pages 31684–31696. Curran Associates Inc., 2022.

[KLSV23]

Yunbum Kook, Yin Tat Lee, Ruoqi Shen, and Santosh S. Vempala. Condition-numberindependent convergence rate of Riemannian Hamiltonian Monte Carlo with numerical integrators. In Conference on Learning Theory, volume 195, pages 4504–4569. PMLR, 2023.

[KN09]

Ravi Kannan and Hariharan Narayanan. Random walks on polytopes and an affine interior point method for linear programming. In Symposium on Theory of Computing, STOC ’09, pages 561–570, New York, NY, USA, 2009. ACM.

[KN12]

Ravi Kannan and Hariharan Narayanan. Random walks on polytopes and an affine interior point method for linear programming. Mathematics of Operations Research, 37(1):1–20, 2012.

[KV06]

Adam Tauman Kalai and Santosh S. Vempala. Simulated annealing for convex optimization. Mathematics of Operations Research, 31(2):253–266, 2006.

[KV24]

Yunbum Kook and Santosh S. Vempala. Gaussian cooling and Dikin walks: the interiorpoint method for logconcave sampling. In Conference on Learning Theory, volume 247, pages 3137–3240. PMLR, 2024.

[KV25a]

Yunbum Kook and Santosh S. Vempala. Faster logconcave sampling from a cold start in high dimension. In Symposium on Foundations of Computer Science, pages 997–1006. IEEE, 2025.

[KV25b]

Yunbum Kook and Santosh S. Vempala. Sampling and integration of logconcave functions by algorithmic diffusion. In Symposium on Theory of Computing, pages 924–932. ACM, 2025.

[KV25c]

Yunbum Kook and Santosh S. Vempala. Zeroth-order logconcave sampling. arXiv preprint arXiv:2507.18021, 2025.

[KVZ24]

Yunbum Kook, Santosh S. Vempala, and Matthew S. Zhang. In-and-Out: algorithmic diffusion for sampling convex bodies. In Advances in Neural Information Processing Systems, volume 37, pages 108354–108388. Curran Associates, Inc., 2024.

[KVZ26]

Yunbum Kook, Santosh S. Vempala, and Matthew S. Zhang. In-and-Out: algorithmic diffusion for sampling convex bodies. Random Structures & Algorithms, 68(3):e70061, 2026.

[KZ25]

Yunbum Kook and Matthew S. Zhang. Rényi-infinity constrained sampling with d3 membership queries. In Symposium on Discrete Algorithms, pages 5278–5306. SIAM, 2025.

[LLV20]

Aditi Laddha, Yin Tat Lee, and Santosh S. Vempala. Strong self-concordance and sampling. In Symposium on Theory of Computing, pages 1212–1222. ACM, June 2020.

[LS93]

László Lovász and Miklós Simonovits. Random walks in a convex body and an improved volume algorithm. Random Structures & Algorithms, 4(4):359–412, 1993.

[LS14]

Yin Tat Lee and√Aaron Sidford. Path finding methods for linear programming: solving linear programs in Õ( rank) iterations and faster algorithms for maximum flow. In Symposium on Foundations of Computer Science, pages 424–433. IEEE, 2014. √ Yin Tat Lee and Aaron Sidford. Solving linear programs with rank linear system solves. arXiv preprint arXiv:1910.08033, 2019.

[LS19]

30

[LV06]

László Lovász and Santosh S. Vempala. Hit-and-run from a corner. SIAM Journal on Computing, 35(4):985–1005, 2006.

[LV07]

László Lovász and Santosh S. Vempala. The geometry of logconcave functions and sampling algorithms. Random Structures & Algorithms, 30(3):307–358, 2007.

[LY21]

Yin Tat Lee and Man-Chung Yue. Universal barrier is n-self-concordant. Mathematics of Operations Research, 46(3):1129–1148, August 2021.

[Nar16]

Hariharan Narayanan. Randomized interior point methods for sampling and optimization. The Annals of Applied Probability, 26(1):597–641, 2016.

[NN89]

Yurii Nesterov and Arkadii Nemirovskii. Self-concordant functions and polynomial time methods in convex programming. preprint, central economic & mathematical institute, ussr acad. Sci. Moscow, USSR, 3(3.2):3–1, 1989.

[NN94]

Yurii Nesterov and Arkadii Nemirovskii. Interior-point polynomial algorithms in convex programming. Society for Industrial and Applied Mathematics (SIAM), January 1994.

[Nua06]

David Nualart. The Malliavin calculus and related topics. Probability and its Applications (New York). Springer-Verlag, Berlin, second edition, 2006.

[Nua19]

David Nualart. Malliavin calculus and normal approximations. XXIII Brazilian School of Probability, 2019.

[Ren88]

James Renegar. A polynomial-time algorithm, based on Newton’s method, for linear programming. Mathematical Programming, 40(1-3):59–93, January 1988.

[SV16]

Sushant Sachdeva and Nisheeth K. Vishnoi. The mixing time of the Dikin walk in a polytope—a simple proof. Operations Research Letters, 44(5):630–634, 2016.

[Tea13]

Ines Thiele et al. A community-driven global reconstruction of human metabolism. Nature Biotechnology, 31(5):419–425, March 2013.

[Vai96]

Pravin M. Vaidya. A new algorithm for minimizing convex functions over convex sets. Mathematical Programming, 73(3, Ser. A):291–341, 1996.

[Vem26]

Santosh S. Vempala. Personal communication, 2026.

A

Self-concordance definitions

We recall several self-concordance used in Lemma 2.2 from [KV24, Definition 1.1]. Definition A.1 (SC, SSC, and LTSC). Let g : int K → Sd++ be a C 2 local metric. We say that g is self-concordant (SC) if for every x ∈ int K and h ∈ Rd , −2∥h∥g(x) g(x) ⪯ Dg(x)[h] ⪯ 2∥h∥g(x) g(x) . We say that g is strongly self-concordant (SSC) if ∥g(x)−1/2 Dg(x)[h] g(x)−1/2 ∥F ≤ 2∥h∥g(x) for every x ∈ int K and h ∈ Rd . We say that g is lower trace self-concordant (LSTC) if tr g(x)−1 D2 g(x)[h, h] ≥ −∥h∥2g(x)

for every x ∈ int K and h ∈ Rd . Note that SSC implies SC since the operator norm is bounded by the Frobenius norm.

31

B

Deferred computations

B.1

Tensor computations

Proof of Lemma 3.9. By Lemma 3.7, hi hj = I1 (ei )I1 (ej ) = I2 (ei ⊗ ej ) + ei ⊗1 ej = I2 (ei ⊗ ej ) + δij . Thus I2 (ei ⊗ ej ) = hi hj − δij . Writing fM = I2 (M ) = I2 (fM ) =

X

ij Mij ei ⊗ ej gives

P

Mij (hi hj − δij ) = hT M h − tr M .

ij

This proves the claim. Proof of Lemma 3.10. For indices i, j, k, Lemma 3.9 and 3.7 give hi hj hk = I2 (ei ⊗ ej ) + δij I1 (ek ) = I3 (ei ⊗ ej ⊗ ek ) + δik hj + δjk hi + δij hk .

Hence, I3 (ei ⊗ ej ⊗ ek ) = hi hj hk − δij hk − δik hj − δjk hi . Multiplying by Cijk and summing over i, j, k gives X X X I3 (C) = C[h⊗3 ] − Cijk δij hk − Cijk δik hj − Cijk δjk hi . i,j,k

i,j,k

i,j,k

Since C is symmetric, each of the sums equals ⟨tr C, h⟩, where (tr C)a = I3 (C) = C[h⊗3 ] − 3⟨tr C, h⟩ = C[h⊗3 ] − 3I1 (tr C). Proof of Lemma 3.11. Write fM =

P

ij Mij ei ⊗ ej and fN =

P

b Cabb .

Therefore

ij Nij ei ⊗ ej . By Lemma 3.7,

P

I2 (M )I2 (N ) = I4 (fM ⊗0 fN ) + 4I2 (fM ⊗1 fN ) + 2I0 (fM ⊗2 fN ) . The zeroth contraction gives I4 (fM ⊗0 fN ) = I4 (M ⊗ N ). For the first contraction, (fM ⊗1 fN )(t, s) =

Z X

Mij ei (t)ej (x)

ij

=

X

 X



Nkl ek (s)el (x) dx

kl

Mij Njk ei (t)ek (s) = fM N (t, s) .

ijk

Since I2 only sees the symmetric part of its matrix coefficient, this contributes I2 (sym(M N )). Finally, Z fM ⊗2 fN =

X ij

 X

Mij ei (x)ej (y)



Nkl ek (x)el (y) dxdy = tr(M N ) .

kl

Combining the three contractions proves the identity. Proof of Lemma 3.12. This is the MSI isometry transported through the tensor-kernel identification. Indeed, ∥In (M )∥2L2 (Ω) = n! ⟨f˜M , f˜M ⟩ = n! ⟨fsym M , fsym M ⟩ = n! ∥sym M ∥2F . The identity ∥fT ∥L2 (Rn+ ) = ∥T ∥F follows from the orthonormality of the functions ei .

32

B.2

Row-map calculus

Lemma B.1 (Derivatives of the row map). For the map R(a) := (a ⊗ a)/∥a∥ and vectors a ̸= 0, b, and c, ∥DR(a)[b]∥ ≤ 3 ∥b∥ ,

∥D2 R(a)[b, c]∥ ≲

∥b∥∥c∥ , ∥a∥

∥D3 R(a)[b, b, b]∥ ≲

∥b∥3 . ∥a∥2

Proof. By direct differentiation, DR(a)[b] = Since ∥v ⊗ w∥ = ∥v∥∥w∥,

b ⊗ a + a ⊗ b ⟨a, b⟩ a ⊗ a. − ∥a∥ ∥a∥3

∥DR(a)[b]∥ ≤ 2∥b∥ + ∥b∥ = 3∥b∥ .

For the second derivative, differentiating DR(a)[b] in direction c, D2 R(a)[b, c] =

b ⊗ c + c ⊗ b ⟨a, c⟩ − (b ⊗ a + a ⊗ b) ∥a∥ ∥a∥3 ⟨b, c⟩ ⟨a, b⟩ ⟨a, b⟩⟨a, c⟩ − a⊗a− (c ⊗ a + a ⊗ c) + 3 a ⊗ a. 3 3 ∥a∥ ∥a∥ ∥a∥5

Each term is bounded by a universal multiple of ∥b∥∥c∥/∥a∥, so ∥D2 R(a)[b, c]∥ ≲ ∥b∥∥c∥/∥a∥. For the third derivative, write R = θΦ, where θ(a) := ∥a∥−1 and Φ(a) := a ⊗ a. Denoting α := ⟨a, b⟩ and β := ∥b∥2 , and differentiating θ and Φ in direction b, β 3α2 α 2 , D θ(a)[b, b] = − + , ∥a∥3 ∥a∥3 ∥a∥5 DΦ(a)[b] = b ⊗ a + a ⊗ b , D2 Φ(a)[b, b] = 2b ⊗ b , Dθ(a)[b] = −

9αβ 15α3 − . ∥a∥5 ∥a∥7 D3 Φ(a)[b, b, b] = 0 . D3 θ(a)[b, b, b] =

Thus the product rule yields D3 R(a)[b, b, b] =

9αβ 15α3
3α2 β
6α − a ⊗ a + 3 − (b ⊗ a + a ⊗ b) − b ⊗ b. 5 7 5 3 ∥a∥ ∥a∥ ∥a∥ ∥a∥ ∥a∥3

Using |α| ≤ ∥a∥∥b∥ and ∥v ⊗ w∥ = ∥v∥∥w∥, we obtain ∥D3 R(a)[b, b, b]∥ ≲ ∥b∥3 /∥a∥2 .

B.3

Elementary frame calculus

The short calculation below is a standard linear-algebra fact about choosing a “nice” orthonormal basis for a smoothly varying subspace. In differential-geometric language, it is part of the standard calculus of orthonormal frames on Stiefel/Grassmann manifolds; see [AMS08] for references. For self-containment, we rederive exactly what we need using only orthogonal changes of basis and differentiation of U T U = I. The next two lemmas use the same elementary idea in two settings. Lemma B.2 treats a one-dimensional path and ensures UtT Ut′ = 0 for every t. Lemma B.3 treats a multi-parameter family and records only the derivatives needed at one base point (in our setting, when z = x). Lemma B.2. For t ≥ 0, consider a path of smoothly varying subspaces Et and an orthogonal matrix Ut for Et (i.e., col Ut = Et ). There exists an orthogonal matrix path Ot such that Ūt := Ut Ot satisfies ŪtT Ūt′ = 0 and Ūt′ = Pt⊥ Ūt′ for every t. 33

Proof. Write Ūt = Ut Ot for an orthogonal matrix Ot . Since Ūt′ = Ut′ Ot + Ut Ot′ , ŪtT Ūt′ = OtT UtT (Ut′ Ot + Ut Ot′ ) = OtT (UtT Ut′ Ot + Ot′ ) . Hence, to ensure ŪtT Ūt′ = 0, it suffices to enforce Ot′ = −UtT Ut′ Ot with O0 = I. Since UtT Ut = I, we have (Ut′ )T Ut = −UtT Ut′ . For St := UtT Ut′ , we have St = −StT . Then, this ODE is written as Ot′ = −St Ot with O0 = I. Since St is smooth on [0, η], a standard result on the matrix ODE theory guarantees the existence of such Ot . Moreover, Ot remains orthogonal: ∂t (OtT Ot ) = (Ot′ )T Ot + OtT Ot′ = −OtT StT Ot − OtT St Ot = 0 . Lastly, as Pt = Ūt ŪtT , Ūt′ = Pt⊥ Ūt′ + Pt Ūt′ = Pt⊥ Ūt′ + Ūt ŪtT Ūt′ = Pt⊥ Ūt′ , which proves the claim. e (v) be Lemma B.3. Let P (v) be a smooth orthogonal projector, defined near v = 0, and let U e smooth orthogonal matrix with col U (v) = im P (v). There exists an orthogonal matrix O(v) such e O satisfies at v = 0 that for every coordinate direction a, b, that U := U

Ua = P ⊥ Ua = Pa U ,

U T Ua = 0 ,

P ⊥ Uab = P ⊥ (Pa U )b .

This also satisfies U T Uab = − 21 (UbT Ua + UaT Ub ) and X

X

X

a,b

a,b

a

∥Uab ∥2F ≲

∥P ⊥ Uab ∥2F +

∥Ua ∥2F

2

.

Proof. All quantities in this proof are evaluated at v = 0 unless stated otherwise. Proof for the e TU e = I, differentiating first-order condition is similar to the previous lemma in overall. Since U P T T T e e e e e e gives U Ua + Ua U = 0, and Ωa := U Ua is skew-symmetric. Choose O(v) := exp(− c vc Ωc (0)). e O, we can check that U T Ua = 0 and Ua = P ⊥ Ua Then, O(0) = I and Oa (0) = −Ωa (0). For U := U in a similar fashion to Lemma B.2. Lastly, since P U = U , we have Pa U + P Ua = Ua . This implies Pa U = P ⊥ Ua . For second derivatives, note that Ua = P ⊥ Ua + P Ua = Pa U + P Ua = Pa U + U U T Ua . Defining Ωa := U T Ua and differentiating Ua = Pa U + U Ωa in direction b at v = 0, Uab = (Pa U )b + U (Ωa )b + Ub Ωa = (Pa U )b + U (Ωa )b . This implies that P ⊥ Uab = P ⊥ (Pa U )b . For Γab := U T Uab , write Uab = P ⊥ Uab + P Uab = P ⊥ Uab + U Γab .

(B.1)

Just as we enforced U T Ua = 0, we will enforce some condition on U so that Γab is symmetric, b (v) := U (v)O(v) with while preserving Ua and ensuring U T Ua = 0. We will replace U (v) by U 1P O(v) = exp( 2 a,b va vb Kab ) where Kab = Kba and each Kab is skew-symmetric. Clearly, O(0) = I, b = U O and estimating it at v = 0 yields Oa (0) = 0, and Oab (0) = Kab . Differentiating U ba = Ua O + U Oa = Ua , U bab = Uab O + Ua Ob + Ub Oa + U Oab = Uab + U Kab . U

34

b TU ba = 0 and U ba = P ⊥ U ba remain valid. Next, from above, P ⊥ U bab = P ⊥ Uab and Hence, U b )b = P ⊥ (Pa U O)b = P ⊥ (Pa U )b , P ⊥ (Pa U b )b . bab = P ⊥ (Pa U so P ⊥ U b b TU bab = U T (Uab + U Kab ) = Γab + Kab . Since Γab = Γba , choosing Kab = Note that Γab = U b symmetric. From now on, we relabel U b − skew(Γab ) is compatible with Kab = Kba and makes Γ ab b as Γ . Differentiating U T U = I in directions a and b and estimating it at v = 0, as U and Γ ab ab T UbT Ua + U T Uab + Uab U + UaT Ub = 0 . T T Since Γab = U T Uab , we have Γab + ΓT ab = −Ub Ua − Ua Ub . Then, since Γab is symmetric,

1 Γab = − (UbT Ua + UaT Ub ) . 2 Therefore, X

∥Γab ∥2F ≲

a,b

The final item follows from Uab

= P ⊥U

X

X

a,b

a

∥Ua ∥2F ∥Ub ∥2F =

ab + U Γab (B.1).

35

∥Ua ∥2F

2

.